Model-agnostic linear-memory online learning in spiking neural networks

Domain: neuroscience
Published: 2026-01-19
Editorial outcome: major_revision

Abstract
Spiking neural networks (SNNs) offer a promising paradigm for modeling brain dynamics and developing neuromorphic intelligence, yet an online learning system capable of training rich spiking dynamics over long horizons with low memory footprints has been missing. Existing online approaches either incur quadratic memory growth, sacrifice biological fidelity through oversimplified models, or lack end-to-end automated tooling. Here, we introduce BrainTrace, a model-agnostic, linear-memory, and automated online learning system for spiking neural networks. BrainTrace standardizes model specification to encompass diverse neuronal and synaptic dynamics; implements a linear-memory online learning rule by exploiting intrinsic properties of spiking dynamics; and provides a compiler that automatically generates optimized online-learning code for arbitrary user-defined models. Across diverse dynamics and tasks, BrainTrace achieves strong learning performance with a low memory footprint and high computational throughput. Critically, these properties enable online fitting of a whole-brain-scale Drosophila SNN that recapitulates region-level functional activity. By reconciling generality, efficiency, and usability, BrainTrace establishes a foundation for spiking network modeling at scale.

Introduction
Spiking neural networks (SNNs) offer a principled route to modeling neural computation at millisecond resolution while retaining the energy efficiency and temporal precision that characterize biological circuits1,2. Yet learning in SNNs remains fundamentally constrained by the challenge of long-horizon temporal credit assignment: synaptic updates must depend on state trajectories that unfold over thousands to millions of discrete time steps. The dominant solution, backpropagation through time (BPTT) with surrogate gradients3–5, computes exact sequence-level gradients but does so offline, requiring storage of all intermediate states and thus incurring memory costs that scale linearly with sequence length T. For biological timescales, where one second may correspond to thousands of solver steps, this memory footprint quickly becomes prohibitive, limiting both the duration of trainable sequences and the size and complexity of SNN models6.

Online learning algorithms for SNNs have emerged as a compelling alternative. By updating parameters in a forward-only fashion, they restrict memory usage to a small set of per-time-step variables and avoid storing entire trajectories7. However, despite encouraging progress, existing methods have not achieved both model-agnostic generality and linear memory scaling simultaneously. State-of-the-art approaches that approach BPTT-level accuracy, including e-prop8 and OSTL9, instantiate large collections of eligibility variables whose memory and computation scale quadratically with network size, rendering them impractical for large-scale spiking networks. Methods that achieve linear memory complexity10–12 typically rely on restrictive neuron abstractions (e.g., simplified leaky integrate-and-fire (LIF) dynamics) that do not generalize to the rich neuronal dynamics needed for biological fidelity. Most critically for practice, there is no analogue of modern automatic differentiation (AD) ecosystems for online learning in SNNs: users cannot write arbitrary spiking dynamics and expect the system to synthesize numerically stable, efficient, and accurate online updates automatically. The absence of such an automated, model-agnostic framework has hampered adoption outside specialized groups and has limited real-world impact.

Here, we introduce BrainTrace, a model-agnostic, linear-memory, and automated online learning system for spiking neural networks that closes this gap. This system provides a theoretically guaranteed, linear-memory online gradient algorithm for training general spiking networks, together with a compiler that automatically generates efficient online-learning code from user-defined SNN models. For efficiency, BrainTrace achieves linear memory scaling relative to the number of neurons while preserving gradient fidelity needed for effective training. For generality, this efficiency extends to a broad class of spiking network dynamics whose synaptic states can be organized along either presynaptic or postsynaptic axes. For usability, BrainTrace’s compiler automatically generates efficient online learning code for any user-defined SNNs, facilitating implementations across various hardware platforms, including central processing units (CPUs), graphics processing units (GPUs), and tensor processing units (TPUs).

We validate BrainTrace across diverse SNN architectures and tasks. Relative to offline surrogate-gradient BPTT and online baselines with quadratic memory growth, BrainTrace delivers orders of magnitude lower memory consumption and high throughput, enabling long-horizon training on commodity accelerators. It attains strong training performance on standard neuromorphic classification benchmarks and brain-simulation settings, where it trains SNNs that reproduce neural dynamics underlying canonical cognitive tasks. Notably, BrainTrace scales to a whole-brain SNN of the Drosophila connectome15,16 and successfully fits whole-brain functional activity patterns17,18. In sum, BrainTrace provides a unified solution to general, efficient, and automated SNN online learning, paving the way toward large-scale, biologically grounded spiking network modeling.

Results
BrainTrace supports a wide range of SNNs with the following discretized dynamics (Supplementary Note A):

|It = f(ht-1, It-1, It) and It = θ ⊛ xt, (1)

where ht ∈ RHd represent the hidden states of H neurons, with each neuron quantified by d state variables; It ∈ RH represent the synaptic currents; and t represent the timestamp; xt ∈ RI represents synaptic inputs (e.g., presynaptic spikes); and θ ∈ RP denotes connection weights between pre- and postsynaptic neurons (Fig. 1).

SNNs described by Eq. (1) are implemented using the AlignPre and AlignPost frameworks19,20 (Fig. 1; see examples in Supplementary Note B). Specifically, in AlignPre, synapses with identical parameters that share the same presynaptic neuron exhibit identical dynamical trajectories (Supplementary Fig. 8), enabling synaptic variables to be aligned with the presynaptic neuronal dimension (Fig. 1A). In AlignPost, when synapses feature exponential dynamics with identical time constants, multiple synapses from presynaptic neurons converging onto a single postsynaptic neuron are effectively represented by a single set of synaptic variables. This is achieved by combining the input to the postsynaptic neuron (Supplementary Fig. 8), allowing synaptic variables to be aligned along the postsynaptic dimension (Fig. 1B).

AlignPre and AlignPost decompose an SNN into two parts (Eq. (1)): (1) an “intrinsic dynamics” part that dictates the temporal evolution of neuronal states through the element-wise activation function f, and (2) a “neuronal interaction” part that converts presynaptic inputs into postsynaptic currents via the interaction operator ⊛. The decomposition can largely reduce computational complexity and memory consumption, enabling a simplified algorithm for SNN online learning (described in the next section).

BrainTrace achieves linear-memory online learning algorithm
The objective of SNN training is to find a group of synaptic weights θ that minimizes a loss function L that quantifies the discrepancy between target values and network outputs. When target values are provided at discrete time points t ∈ T, we compute the weight gradient by:

∇θ L = ∑t∈T ∂Lt/∂ht ϵt, (2)

where Lt represent the instantaneous loss at time t, ∂Lt/∂ht denotes the learning signal that arrival at the hidden state ht, and ϵt is the eligibility trace indicating how parameters influence the hidden state, which is computed as13:

ϵt = ∂ht/∂ht-1 ϵt-1 + ∂ht/∂θt, (3)

where the hidden Jacobian ∂ht/∂ht-1 = Dt + Jt consists of a diagonal matrix Dt derived from intrinsic dynamics f, and a square matrix Jt derived from neuronal interactions ⊛ (see Methods). The hidden-to-weight Jacobian ∂ht/∂θt captures the effect of weight changes on hidden states. Importantly, Eq. (3) imposes a memory overhead of O(H3) (Fig. 2A).

By leveraging SNN’s inherent properties (see Methods), we can reduce the complexity to O(H). Without loss of generality, we consider each neuron has only 1 state variable with d = 1 (but see Methods for generalized results for d > 1). First, we found that a diagonal approximation of the hidden Jacobian captures over 99% of the cosine similarity with the full Jacobian (∂ht/∂ht-1 ≈ Dt, Fig. 2B). This fact aligns with earlier practice8, which suggests that ignoring the gradient flow through recurrent connections in ALIF neurons does not significantly hinder learning performance. Second, we proved (Supplementary Note C) that the hidden-to-weight Jacobian can be decomposed as ∂ht/∂θt = Dt ⊗ xt, where Dt is a diagonal matrix with entries Dt,ii = ∂ht,i/∂It,i, and ⊗ represents the Kronecker product. These observations enable us to simplify the eligibility trace in Eq. (3) to:

ϵt ≈ ϵtD = Dtϵt-1 + diag(Dt) ⊗ xt, (4)

which we referred to as the diagonal-approximated RTRL (D-RTRL) (Fig. 2C). This diagonal approximation reduced the memory complexity from O(H3) to O(H2).

We further observed that ϵtD can be expressed as a summation of rank-one matrices: ϵtD = ∑t k=1 bt ⊗ xk, where bt = ∑i=k+1 diag(Di) ∘ diag(Dk) (Supplementary Fig. 9). Since inputs xt maintain consistent signs across time steps (whether as binary spikes in AlignPost or synaptic conductance in AlignPre (Fig. 1)), we proved that the sum of multiple rank-one matrices can be approximated as a single rank-one matrix between summed vectors (Supplementary Note D). Consequently, the eligibility trace can be simplified to:

ϵt ≈ ϵtD ≈ ∑t k=1 bt ⊗ xk, (5)

which still preserves a high cosine similarity with the full eligibility trace across LIF (Fig. 2D), ALIF (Fig. 2E), and other SNN architectures (Supplementary Fig. 15).

In summary, this approximation yields the pre-post eligibility trace propagation (pp-prop) algorithm (Fig. 2F):

ϵt ≈ ϵt f ∘ ϵt x, (6)

ϵt x = αϵt-1 x + xt, (7)

ϵt f = αdiag(Dt) ∘ ϵt-1 f + (1 - α)diag(Dt), (8)

where ∘ denotes element-wise multiplication, and α (0 < α < 1) is the smoothing factor corresponding to time constant τ = -Δt/ln(α). ϵt x ∈ RI and ϵt f ∈ RH represent eligibility traces tracking pre- and post-synaptic activity, respectively. The eligibility trace in pp-prop operates temporally, locally, and in an unsupervised way, effectively achieving memory complexity of O(H).

BrainTrace facilitates automated online learning compilation
The widespread success of ANNs primarily stems from the automation achieved in frameworks such as PyTorch21 and TensorFlow22 for their training algorithms, namely backpropagation and BPTT algorithms. To provide similar training convenience for SNNs, we introduced a compiler (Fig. 3A–E) that automates the translation of SNNs into efficient online learning implementations of D-RTRL and pp-prop, and subsequently deploys both models and generated algorithms across various hardware backends, including CPUs, GPUs, and TPUs (Fig. 3).

At the Primitive layer, the system encapsulates three primitives that capture the essential components of general SNNs (as formalized in Eq. (1), Fig. 3A). Specifically, Hidden State primitives represent internal neuronal and synaptic variables, such as membrane potentials and synaptic conductances; Parameter primitives denote trainable parameters, such as synaptic weights; while Connection primitives describe network connectivities, including dense, convolutional, and element-wise operations. By flexibly composing these primitives, the Complex Network layer supports the construction of complex SNNs with sophisticated dynamics and topologies (Fig. 3B).

The core innovation of the compiler lies in the ETraceGraph and ETraceAlgorithm layers. At the ETraceGraph stage, user-defined Python-based models are transformed into symbolic intermediate representations (IRs), which explicitly encode data dependencies between model parameters and neuronal states. The compiler analyzes IRs to extract higher-level computational graphs that describe how each model parameter, θ, influences the hidden state through specific connections ⊛. The resulting graph captures all critical information necessary for computing eligibility traces, including the hidden-to-hidden Jacobian (diag(Dt)) and hidden-to-weight Jacobians (diag(Dt) and xt) (Fig. 3C).

Subsequently, at the ETraceAlgorithm layer, the compiler automatically generates optimized online learning implementations based on the user-selected algorithm (Fig. 3D). The generated code includes both the logic for updating eligibility traces during state evolution (see Eq. (4) for D-RTRL and Eqs. (6)–(8) for pp-prop) and the gradient computation triggered by learning signals (see Eq. (2)). During this code generation process, the compiler also performs advanced optimizations, such as reusing shared eligibility traces within the same neuron population to minimize memory footprint, and vectorizing hidden-group Jacobian computations to maximize computational efficiency (Supplementary Note I).

Finally, leveraging just-in-time (JIT) compilation via JAX and XLA23, both user-defined models (Fig. 3B) and compiler-generated algorithms (Fig. 3D) are automatically deployed onto the target hardware platform, including CPUs, GPUs, and TPUs (Fig. 3E).

Through this hierarchical online learning compilation, BrainTrace provides a unified and automated solution for computing forward gradients (see comparison of Fig. 3F for BPTT and Fig. 3G for online learning) without introducing additional programming difficulties (see comparison of Fig. 3H for BPTT and Fig. 3I for online learning; Supplementary Note N). This end-to-end online learning system will greatly lower the barrier to entry for using online learning in SNNs, making advanced algorithms more accessible and deployable in real-world applications.

Algorithmic evaluations on approximation, long-term dependency, memory, and speed
We first assessed the effectiveness of BrainTrace compiled algorithms in approximating the gradient computation. Evaluations were performed over different spiking datasets (including IBM DVS Gesture24, Spiking Heidelberg Digits (SHD)25, and Neuromorphic-MNIST (N-MNIST)26 datasets), neuronal dynamics (ranging from LIF to ALIF neurons), and synapses (including Delta, exponential, short-term depression (STD), and short-term plasticity (STP)) (Supplementary Note B). We observed consistently high cosine similarity between the single-step hidden Jacobian (1/T ∑T ∂hi+1/∂hi; Fig. 4A) and long-range hidden Jacobian (∂hT/∂h1; Fig. 4B) compared to BPTT-computed Jacobian. Notably, high similarity values were especially prominent in networks with complex intrinsic dynamics, such as ALIF-based architectures (Fig. 4B). We further found that both D-RTRL and pp-prop achieved high gradient approximation accuracy relative to BPTT-computed gradients (Fig. 4C-D), validating their effectiveness for the gradient approximation. Moreover, both algorithms exhibited lower gradient similarities on N-MNIST (Fig. 4C-D), indicating some dataset dependence. Further investigation revealed several key insights about factors affecting the approximation accuracy (Supplementary Note G). Network size showed no significant impact on approximation accuracy (Supplementary Fig. 11). However, approximation declined greatly when either feedforward (Supplementary Fig. 12) or recurrent (Supplementary Fig. 13) weight strength increased to levels that produced unrealistically high firing rates. Additionally, deeper recurrent layers led to lower approximation scores (Supplementary Fig. 14), suggesting a preference for shallow network architectures in our algorithms.

We then assessed whether BrainTrace algorithms can capture long-term dependencies in datasets, a capability that is critical for effective SNN modeling. We trained recurrent spiking networks (Supplementary Eq. 73) on the delayed match-to-sample (DMTS) task27, in which the network is required to maintain directional information from the sample phase over a long delay period before making a decision based on the test period stimulus (see Methods). We found that while D-RTRL, pp-prop, and BPTT all effectively solve the task with surrogate gradients3 (Fig. 4E), only D-RTRL and pp-prop can learn successfully without surrogate gradients (Fig. 4F). We then investigated how our algorithm achieves long-term dependency learning. We found that both pp-prop and D-RTRL maintain temporal information through eligibility traces, but in distinct ways. In D-RTRL, ϵg responds to incoming spikes, increasing with each presynaptic event; ϵa captures slow Jacobian dynamics, enabling information to persist across delays; and ϵv combines information from ϵg and ϵa, remaining sensitive to input changes like ϵg while decaying gradually like ϵa (Fig. 4G). In pp-prop, the presynaptic trace ϵpre is driven solely by presynaptic spikes, while ϵpost,g tracks postsynaptic currents and remains constant; ϵpost,a and ϵpost,v evolves slowly, propagating long-term postsynaptic Jacobian information until learning signals arrived during testing (Fig. 4H). We also analyzed the differences in synaptic connectivity learned by pp-prop compared to BPTT and D-RTRL. pp-prop maintained similar weight distributions with D-RTRL and BPTT, yet it displayed distinct column structures, particularly in recurrent weights (Supplementary Fig. 21). Interestingly, similar column structures have been observed in the synaptic connectome of layer 4 in the somatosensory cortex (Fig. 3E in 28), suggesting that pp-prop may capture biologically relevant connectivity patterns.

Finally, we assessed the memory consumption and computational speed of pp-prop compared with other algorithms. Theoretically, for memory consumption, pp-prop requires only O(BN) memory, since it stores just two traces capturing pre- and postsynaptic activities, whereas other online learning methods demand at least O(BN2). Here, B is the batch size. For computation, pp-prop needs O(TBN2) operations for gradient updates but only O(TBN) operations for eligibility traces. In contrast, OSTL, e-prop, and D-RTRL incur quadratic costs for eligibility trace computation due to their reliance on matrix-based updates. Empirically, we ran experiments on the Gesture dataset to compare memory and speed of three training methods, BPTT, D-RTRL, and pp-prop, across different sequence lengths using LIF networks (Supplementary Eq. 70). We found that BPTT’s GPU memory usage grew linearly with time steps and crashed at 400 steps (Fig. 4J). By contrast, both D-RTRL and pp-prop required constant memory regardless of sequence length (Fig. 4J), with pp-prop being 35 times more memory-efficient than D-RTRL (0.5 GB vs. 17.5 GB). For computation time, all algorithms scaled linearly with time steps, but pp-prop ran slightly faster than BPTT and achieved a ten-fold speedup over D-RTRL (Fig. 4K). These memory and speed advantages were consistent on CPU (Supplementary Fig. 19) and TPU platforms (Supplementary Fig. 20). Additionally, we tracked model accuracy during training and found that all three algorithms achieved higher test accuracy as sequence length increased (Fig. 4L).

Performance evaluations on neuromorphic classification tasks
Spiking networks are widely applied in neuromorphic settings1,6. To assess the effectiveness of pp-prop for such applications, we conducted extensive experiments across three standard benchmarks: N-MNIST26, SHD25, and IBM DVS Gesture24. We evaluated multiple SNN architectures: (i) a two-layer RLIF network29 comprising LIF neurons with delta synapses and recurrent connections (Supplementary Eq. 70); (ii) a two-layer RadLIF network29 with adaptive LIF neurons, delta synapses, and recurrence (Supplementary Eq. 71); and (iii) a one-layer event-based gated recurrent unit (EGRU)30 (Supplementary Eq. 82).

Since performance depends on many confounders (such as data preprocessing, network architecture, and training schedule), we varied only the learning algorithm while keeping all other factors fixed. For each algorithm, we performed targeted hyperparameter searches to obtain strong settings (see Supplementary Note L). We report mean ± s.d. accuracy over at least three independent runs (Table 1). Across datasets and architectures, pp-prop matched BPTT on most settings while using an online, memory-efficient updater rule.

Against state-of-the-art online learners, including OSTL9, e-prop8, ETLP31, OSTTP32, OTPE11, S-TLLR33, FPTT34, and OTTT10, pp-prop consistently achieved the strong performance. Notably, on SHD, pp-prop exceeded the strongest baseline by more than 10 percentage points (Table 1).

Compared with state-of-the-art offline learning approaches, our pp-prop algorithm achieves comparable performance to the latest models that employ advanced architectural and training optimizations. For instance, on the SHD dataset, pp-prop even surpasses the best-performing models incorporating delay-based learning35. On the DVS Gesture benchmark, its performance closely matches that of models leveraging batch normalization36 and specialized neural dynamics design37.

Training excitatory-inhibitory spiking networks for cognitive tasks
BrainTrace is theoretically applicable to a broad spectrum of spiking networks, including biologically realistic circuits with complex neuronal dynamics. To demonstrate this, we developed a biologically-informed excitatory-inhibitory spiking network (Fig. 5B) incorporating: four-variable generalized leaky integrate-and-fire (GIF) neurons38 for complex neuronal dynamics, distinct excitatory and inhibitory (E/I) populations, conductance-based synapses39, and sparse connectivity40. Our network comprised 800 GIF neurons, with a 4:1 excitatory-to-inhibitory ratio and a 10% random connection probability, each featuring conductance-based exponential synapse dynamics (Supplementary Eq. 81). We trained this network on an evidence-accumulation task41, where animals navigate a virtual linear track and encounter seven visual cues (Fig. 5A). The task required animals to integrate left/right visual cues and, upon reaching a T-junction, select the corresponding direction. Correct decisions based on majority cue distribution were rewarded.

Our analysis revealed that when the readout layer decodes recurrent membrane potential dynamics for accuracy, both D-RTRL and e-prop achieved stable convergence on this task, whereas BPTT failed without surrogate gradients (Fig. 5C, Supplementary Fig. 25). Despite the network’s small size, BPTT required over 10 GB of device memory, and D-RTRL consumed more than 2 GB. In contrast, pp-prop maintained high performance with a memory footprint of less than 0.2 GB (Fig. 5D). Post-training analysis of pp-prop network dynamics revealed highly selective, temporally precise spiking patterns. During the evidence accumulation phase (0-1100 ms, Fig. 5E), sequentially presented visual cues induced orderly transitions between transient and clustered neuron population activities (Fig. 5F). Neurons in each population predominantly utilized a time-to-first-spike coding strategy for sensory processing. Notably, within each selective population, neuronal responses to external stimuli displayed clear rank ordering, with neurons firing sequentially in response to stimuli (Fig. 5F). Remarkably, these neural dynamics closely align with neuronal activities observed in the posterior parietal cortex of mice performing similar evidence accumulation tasks (Fig. 2a in ref. 41).

Modeling Drosophila resting-state functional connectivity via whole-brain spiking network training
The linear memory complexity of BrainTrace provides a crucial advantage for training large-scale spiking networks. To illustrate this, we trained a whole-brain spiking model of the Drosophila brain constrained by the FlyWire anatomical connectome15, to reproduce both resting-state whole-brain calcium imaging signals and their inter-regional functional connectivity17,18. Although the FlyWire connectome15 offers a detailed anatomical map of neuronal wiring, direct implementation of this connectivity in spiking network models16 fails to accurately capture either the resting-state neural activity in each neuropil (Fig. 6G) or the cross-neuropil correlations (Fig. 6K) recorded in experiments. We therefore trained the synaptic weights of the connectome-based network16 using our pp-prop algorithm.

We first converted the 17 minutes of calcium signals recorded from each brain region (Fig. 6A) into deconvolved firing rates (Fig. 6B), and divided the entire dataset into a training set of 80% and a testing set of 20%. During training, the network received the experimental deconvolved firing rate at time t (ExpFRt) as input and predicted the firing rate at time t+1 (SimFRt+1). The loss function was defined as the mean square error between simulated and experimental firing rates. The pp-prop algorithm then continuously adjusted synaptic weights online according to the error gradients, enabling the network to autonomously learn the complex spatiotemporal patterns underlying resting state dynamics. This model consumes approximately 8.9 GB of GPU memory with pp-prop, and the loss converges slowly over epochs (Fig. 6D). In contrast, training with D-RTRL and BPTT exceeds the 32 GB capacity of a single GPU (Supplementary Fig. 24).

At the neuropil level, we evaluated the model’s dynamics in the mushroom body to assess training-induced improvements (Fig. 6E–G). In the warmup phase, we initialized whole-brain dynamics by feeding the first 84 seconds of experimental firing rates into the model; thereafter, without external inputs, the model generated ongoing neural activity solely from its previous timestep’s neuropil state. Before training, the model failed to reproduce any experimentally observed spontaneous fluctuations (Fig. 6G); after training, however, its output not only exhibited activation patterns highly consistent with the experimental recordings ("Train” phase in Fig. 6F) but also sustained similar oscillation patterns on previously novel test data ("Test” phase in Fig. 6F). These results demonstrate the model’s ability to generalize to untrained neural dynamics. To quantify improvements in all brain regions across the entire Drosophila brain, we computed Pearson correlations between model-generated and experimental firing rates for each of 68 regions, both pre- and post-training. Training produced a marked uplift in regional fit, and statistical tests confirmed that nearly every region exhibited a significant post-training correlation increase (Fig. 6H).

At the whole-brain level, we further compared whole-brain functional connectivity among experimental recordings, the untrained model’s spontaneous activity, and the trained model’s output (Fig. 6I–K). Prior to training, neurons within each neuropil exhibited spontaneous fluctuations with similar activity patterns (Supplementary Fig. 22), leading to high correlations across neuropils (Fig. 6K). After training, the model dynamics generated highly consistent functional connectivity strengths and patterns (Fig. 6J and Supplementary Fig. 23) compared to experimental data, regardless of whether during the training or testing phases.

Methods
Inspecting inherent properties of spiking networks
Property 1: intrinsic neuronal dynamics dominate the hidden Jacobian. A typical SNN exhibit two types of recurrent connections. The first type comprises Dt recurrent spiking connections that propagate action potentials between neurons (see ⊛ in Eq. (1) and examples in Supplementary Note B). The second type Jt consists of intrinsic neuronal dynamics, such as leaky terms in neuronal and synaptic processes (see f Eq. (1) and examples in Supplementary Note B). As a result, the hidden Jacobian is composed of two parts: ∂ht/∂ht-1 = Dt + Jt, where Dt ∈ RH×H is a diagonal matrix resulting from the element-wise intrinsic dynamics, and Jt is the square matrix derived from the recurrent spiking connections. Dt is actually a block diagonal matrix of dimensions RHd×Hd, consisting of H-by-H blocks. Each diagonal block, Dt,ii, is a d-by-d square matrix, where d is the number of state variables aligned to neuron i. For example, refer to Supplementary Note B. Without loss of generality, all subsequent derivations assumed = 1, making Dt a diagonal matrix of dimensions RH×H.

A defining characteristic of SNNs is their sparse activation patterns, with neurons firing at relatively low frequencies. Empirical studies have shown that biological neurons typically emit between 0.2 and 5 spikes per second42, with only 1–2% of neurons in a cortical network active at any given moment43,50. This fact leads to the hypothesis that Jt exerts limited influence and can be ignored in the hidden Jacobian, since Jt is highly sparse in both space and time. Earlier practices, such as e-prop8 and RFLO51, have also demonstrated that neglecting the gradient flow through recurrent spiking connections in specific neuronal dynamics does not significantly hinder learning performance. To further validate this hypothesis, we conducted empirical evaluations across various SNN architectures with different neuronal and synaptic dynamics (Supplementary Note B). Our findings consistently indicated that the hidden Jacobian of most SNNs is predominantly shaped by the diagonal Jacobian Dt induced by intrinsic dynamics. The diagonal approximation demonstrates a remarkably high cosine similarity with the full hidden Jacobian (see Fig. 2B). Notably, under biologically plausible firing regimes (fraction of spiking neurons <10%), the cosine similarity exceeds 99%.

Proposition 1. In sparsely activated SNNs, the hidden Jacobian ∂ht/∂ht-1 is dominated by the diagonal Jacobian Dt induced by intrinsic neural dynamics. This can be expressed as:

∂ht/∂ht-1 ≈ Dt. (9)

Property 2: physical connections yield decomposable hidden-to-weight Jacobian. Biological SNNs are characterized by physically constrained connections, reflecting the actual synaptic formations between pre- and postsynaptic neurons. Specifically, in AlignPre and AlignPost abstractions, this structural property can be mathematically expressed as a parameterized transformation from inputs to hidden vectors, followed by the element-wise activation applied to the resulting hidden vectors (see Eq. (1)). Such mathematical operation leads to an important implication: the hidden-to-weight Jacobian ∂ht/∂θt, which quantifies the sensitivity of the network’s output to changes in synaptic weights, can be decomposed into a Kronecker product of a vector and a diagonal matrix (see Theorem 2 and its proof in Supplementary Note C; this observation has also been illustrated in ref. 52 for rate-based recurrent neural networks).

Theorem 2. Given Eq. (1), let It = θ ⊛ xt, then the gradient of ht with respect to θ at time t is

∂ht/∂θt = Dt ⊗ xt, (10)

where Dt is a diagonal matrix with ∂ht,i/∂It,i as its diagonal entry Dt,ii, and ⊗ is the Kronecker product.

Dt is actually a block diagonal matrix of dimensions RHd×H, consisting of H-by-H blocks. Each diagonal block, Dt,ii, is a d-by-1 matrix. For example, refer to Supplementary Note B. With loss of generality, all subsequent derivations assume d = 1, making Dt a diagonal matrix of dimensions RH×H.

Property 3: sign-preserving neuronal outputs induce low-rank structure. Another unique characteristic of SNNs is that the outputs of each neuron, corresponding to inputs xt, maintain a consistent sign across all time steps. For AlignPost projections, x denotes the binary spike vector (see Supplementary Eq. 63, 64); while for AlignPre projections, xt denotes the synaptic conductance (see Supplementary Eq. 61, 62). In both cases, the elements in xt are sign-preserving and consistently positive across all time steps. Additionally, for the monotonous activation function f, its derivative Dt also maintains a consistent sign throughout all times. This distinctive property leads to an important mathematical implication: the summation of multiple rank-one matrices can be approximated by the rank-one outer product of the summed vectors. Theorem 3 formalizes this concept, with its theoretical proof presented in Supplementary Note D and empirical validation on various spiking dynamics using real-world neuromorphic datasets provided in Supplementary Note F.

Theorem 3. Let a ∈ Rm and b ∈ Rn be vectors where all elements in ai and bi are of the same sign (i.e., they are all positive or negative). Then the cosine similarity between ∑r i=1 ai ⊗ bi and (∑r i=1 ai) ⊗ (∑r i=1 bi) approaches 1 as r → ∞:

lim r→∞ cos(∑r i=1 ai ⊗ bi, (∑r i=1 ai) ⊗ (∑r i=1 bi)) = 1. (11)

Derivation of online learning algorithms
Theoretical basics of RTRL. Let us consider biological SNNs whose dynamics are governed by Eq. (1). At each time step t ∈ T, we assume the objective is to optimize the loss L = ∑t Lt with respect to the parameters θ, utilizing the gradient

∇θ L = ∑t∈T ∂Lt/∂θt = ∑t∈T ∑k=1 t ∂Lt/∂ht ∂ht/∂hk ∂hk/∂θk, (12)

where θk represents the parameters utilized at time t, ∂Lt/∂ht ∂ht/∂hk measures how θ at step k affects the loss at step t≥k, and ∂ht/∂hk transports the error through time from step t to step k. When k ≪ t, the transport captures a long-term dependency; otherwise, a short-term dependency.

To compute Eq. (12), two methods are commonly used, that is BPTT and RTRL13, which iteratively perform the error propagation on the unfolded computation graph through time backwardly and forwardly, respectively (also refer to Supplementary Fig. 16). For the RTRL, the gradient is updated in the forward direction:

∇θ L = ∑t∈T ∂Lt/∂ht ϵt, (13)

ϵt = ∂ht/∂ht-1 ϵt-1 + ∂ht/∂θt, (14)

where the eligibility trace ϵt ∈ RH×H×I calculates the effect of the parameters ϵt. To better illustrate the difference between BPTT and RTRL, let’s examine the gradients of the loss at time t with respect to the parameter at time k, denoted as ∂Lt/∂θk, where t≥k. Supplementary Fig. 16 demonstrates that BPTT’s backward computation of gradients is significantly more efficient than RTRL’s forward computation. This is because the product of Jacobians from the end to the beginning of the composed function only requires O(H) memory storage. In contrast, RTRL’s memory complexity is O(H×H×I) as it computes Jacobian products forwardly. To reduce the intrinsic complexity of RTRL, we should take full consideration of SNN properties.

Derivation of the D-RTRL algorithm from RTRL. In the following, we concentrate on the scenario where the hidden state dimension d = 1 to develop the D-RTRL algorithm.

As stated above, the raw RTRL has a high memory complexity. By substituting Eqs. (9) and (10) into Eq. (14), we can approximate ϵt as:

ϵt ≈ Dt ϵt-1 + Dt ⊗ xt, (15)

Eq. (15) results in a very sparse matrix with O(H×H×I) memory complexity. For compact storage, the above equation is equivalent to the following form:

ϵt ≈ Dtϵt-1 + diag(Dt) ⊗ xt, (16)

where diag(Dt) ∈ RH is the vector storing the diagonal elements of Dt. This results in ϵt having the shape of RH×I.

In practice, Dt is stored as a vector diag(Dt) ∈ RH×1. Consequently, Eq. (16) can be reformulated using element-wise multiplication. Combining with Eq. (13), we can express the learning rule of the D-RTRL algorithm as:

∇θ Lt ≈ ∂Lt/∂ht diag(Dt) ϵt-1 + diag(Dt) ⊗ xt, (17)

where the D-RTRL algorithm is a memory-efficient approach that avoids the O(H3) complexity of RTRL. This efficiency scaling arises because the eligibility trace ϵt only stores historical information for each active synaptic connection. In sparse neural networks, D-RTRL achieves additional memory savings by eliminating traces for zero-weight connections. Convolutional layers further benefit from this approach, as their parameter count is typically smaller than the number of pre- and postsynaptic neurons when using small kernels.

Derivation of the pp-prop algorithm from D-RTRL. Given the learning rule of the D-RTRL algorithm in Eq. (17), let’s consider the influence of θ at step k to the loss Lt at step t (t≥k). It can be approximated by:

∂Lt/∂θk ≈ ∂Lt/∂ht (Dt ⊗ xk), (18)

Therefore, at the time t, the parameter gradient accumulated up to time t has the form of:

∇θ Lt = ∑t k=1 ∂Lt/∂ht (Di ⊗ Dk ⊗ xk), (19)

Here, we assume that a learning target is provided at each time step k ∈ [1,⋯ ,t]. Moreover, ∂Lt/∂ht ∈ RH is defined as a row vector in Eqs. (18) and (19).

Since Di and Dk are diagonal matrices, we can only keep their diagonal elements as vectors, and their matrix multiplication operations are equivalent to the element-wise Hadamard products of diagonal vectors:

∇θ Lt ≈ ∂Lt/∂ht (diag(Dt) ∘ diag(Dt-1) ∘ diag(Dk+1) ∘ diag(Dk) ⊗ xk), (20)

Equation (20) results in significant memory overhead because the summation requires storing t pairs of vectors representing pre- and postsynaptic activity. This memory requirement is equivalent to BPTT, with complexity scaling linearly as time progresses.

We found a solution to this issue by fully leveraging the unique properties of SNNs. The key insight lies in the nature of neuronal output of SNNs: (1) Binary events: neurons in SNNs output binary signals (spikes); (2) Non-negativity: these spike events are always non-negative (see Theorem 3). This theorem allows us to approximate Eq. (20) by

∇θ Lt ≈ ∂Lt/∂ht (∑t k=1 bt ⊗ xk), (21)

where bt = diag(Dt) ∘ ⋯ ∘ diag(Dk+1) ∘ diag(Dk).

To avoid the exploding problem, the eligibility trace should prioritize recent activities rather than retain a complete historical record over a very long sequence. We implemented this by modifying Eq. (21) using two key techniques: exponential smoothing for bt and a low-pass filter for x. The former computes a running average of bt online, while the latter calculates a running summation of x. These adaptations form the foundation of our pp-prop algorithm, whose learning rule can be expressed as follows:

∇θ L ≈ ∑t=T ∂Lt/∂ht ∘ ϵt f ⊗ ϵt x, (22)

ϵt x = αϵt-1 x + xt,

ϵt f = αdiag(Dt) ∘ ϵt-1 f + (1 - α)diag(Dt),

where α (0 < α < 1) is the smoothing factor, corresponding to the time constant τ = -Δt/ln(α).

Generalization to hidden states with multiple variables (d > 1). The derivations above assumed d = 1, but generalize straightforwardly to d > 1. When d > 1, the matrices transform into block structures as follows (see examples in Supplementary Note B):

(cid:129) ht ∈ RHd×1: A block vector with H blocks of shape Rd×1, where block i contains the state variables of neuron i.

(cid:129) Dt ∈ RHd×Hd: A block diagonal matrix with H diagonal blocks of shape d×d, where Dt,ii is the Jacobian between state variables of neuron i.

(cid:129) Dt f ∈ RHd×H: A block diagonal matrix with H diagonal blocks of shape d×1, where Dt,ii = ∂ht,i/∂It,i.

(cid:129) ∂Lt/∂ht ∈ RHd×1: A block vector with H blocks of shape Rd×1.

The online learning framework in Eqs. (17) and (22) remain unchanged; only the ∘ operation differs. For d = 1, ∘ denotes element-wise scalar multiplication. For d > 1, ∘ it becomes block-wise multiplication, where each block undergoes d-dimensional matrix or inner product operations.

To illustrate, consider the learning rule for a single weight θji connecting neurons i and j. For D-RTRL, Eq. (17) becomes:

∇θ ji L ≈ ∑t∈T ∂Lt/∂ht,j ϵt,ji, (23)

ϵt,ji = Dt,jj ϵt-1,ji + Dt,jj f ∘ xi,t,

where xt ∈ R is a scalar, ϵt,ji ∈ RHd×H is a block matrix with H × H blocks where each ϵt,ji ∈ Rd×1, and ⟨(cid:6),(cid:6)⟩: Rd×1 × Rd×1 → R denotes the inner product.

Similarly, for pp-prop, Eq. (24) becomes:

∇θ ji L ≈ ∑t=T ∂Lt/∂ht,j ∘ ϵt f,j ⟨ϵt x,i⟩, (24)

ϵt f,j = αDt,jj f ϵt-1 f,j + (1 - α)Dt,jj f,

where ϵt f,j ∈ RHd×1 is a block vector with each ϵt f,j ∈ Rd×1, and ϵt x,i ∈ R tracks the presynaptic eligibility trace of neuron i.

In summary, extending the derivations from d = 1 to d > 1 simply replaces scalars with blocks and element-wise multiplications with block-wise multiplications. The online learning framework remains structurally identical, while gaining the ability to represent richer neuron dynamics and more complex temporal dependencies.

Learning signals in deep recurrent layers. At time t, we compute the learning signal ∂Lt/∂ht,l at layer l using backpropagation through the stacked hidden layers. To maintain linear memory complexity, we consider two approaches based on when the top-down learning signal reaches layer l.

The first approach uses the learning signal arriving at layer l at time t, capturing instantaneous dependencies across hierarchical layers:

∇θ l L = ∑t∈T ∂Lt/∂ht,l ∘ ϵt,l. (25)

The second approach uses the learning signal arriving at layer l at time t-1, capturing dependencies across time:

∇θ l L = ∑t∈T ∂Lt/∂ht,l-1 ∘ ϵt,l-1. (26)

Experimental details for cognitive tasks
Surrogate gradient function. The spike generation function (Supplementary Eq. 36) uses a Heaviside step function with a transition at the firing threshold vth. Since this function is non-differentiable, we employ a surrogate gradient approach53 to enable gradient-based training.

During the forward pass, spikes are generated using the Heaviside function:

zt = spike(x) = H(vt - vth),

where x = vt - vth. During the backward pass, we approximate the derivative using a triangular surrogate3:

spike'(x) = ReLU(α ⋅ (width - |x|)),

where α = 0.3 controls the gradient amplitude, and width = 1.0 controls the gradient width. This produces a triangular function centered at x = 0, with peak value α ⋅ width, and zero gradient for |x| > width.

Weight initialization. For the networks defined in Supplementary Eqs. 70–78, we initialize input, recurrent, and readout weights using a scaled Gaussian distribution:

Wji ∼ N(0, s/sqrt(Nin)),

where Wji is the weight from neuron i to neuron j, Nin is the number of presynaptic neurons, and s is a scaling factor controlling the overall weight magnitude. This initialization helps keep activation and gradient variances well-behaved, improving training stability.

For networks with explicit excitatory/inhibitory (E/I) separation (Supplementary Eq. 81), we use the same scaled Gaussian distribution but enforce non-negative weights:

Wji ∼ |N(0, s/sqrt(Nin))|.

Taking the absolute value ensures that excitatory neurons produce only positive outgoing weights and inhibitory neurons only negative ones, in accordance with Dale’s principle. This constraint preserves the desired biological structure while maintaining an appropriate weight scale.

Delayed match-to-sample task. We tested the long-term dependency learning capability of our algorithms by training a recurrent spiking network (Supplementary Eq. 73) on a delayed match-to-sample (DMTS) task27,54. This task comprises four phases: fixation, sample, delay, and test (Supplementary Fig. 17). During the sample and test phases, one of eight motion directions is presented for 500 ms. The delay phase (>1000 ms) separates the sample and test, during which no stimulus is presented. The network must maintain sample direction information across this delay, then determine whether the test stimulus matches. A leaky readout layer processes network activity during the test phase to produce the match/non-match decision.

The network receives input from 100 motion-direction-tuned neurons covering 360 degrees, simulating direction-selective cells in the visual cortex. The firing rate ut,i of input neuron i follows a von Mises tuning curve:

ut,i = Aexp(κ(cos(θ - θi,pref) - 1)), (27)

where θ is the stimulus direction, θi,pref is the preferred direction of neuron i, κ = 2.0 is the concentration parameter, and A = 40 Hz is the maximum firing rate. A constant background rate rbg = 1 Hz is maintained throughout all task phases. Spike trains are generated via a Poisson process with rate ut,i + rbg during stimulus periods and rbg otherwise (Supplementary Fig. 17).

Evidence accumulation task. Evidence accumulation is a fundamental paradigm for studying decision-making mechanisms. We adopted a task design from Morcos et al.41, where a mouse navigates a virtual T-maze and encounters visual cues on both sides (Fig. 5A). At the T-junction, the mouse must choose the direction corresponding to the majority of cues, requiring it to accumulate and retain evidence from each side until the decision point.

We simulated this task using 100 input neurons divided into four groups of 25 neurons each: left-cue encoding, right-cue encoding, recall-cue encoding, and background noise. Cue and recall stimuli are modeled as Poisson spike trains at 40 Hz, with each stimulus lasting 150 ms separated by 50 ms intervals. Background neurons fire at a constant 10 Hz throughout the trial. Following the cue presentation, a 1000 ms delay precedes the recall stimulus (Fig. 5A).

The leaky rate readout. The networks in Fig. 4E–H and Fig. 5 use a leaky readout mechanism with output neurons corresponding to each task label. These neurons receive linear projections from the recurrent layer and perform decoding based on membrane potentials rather than spikes.

The readout dynamics follow:

τout dy/dt = -y + Wout rrec + bout, (28)

where y is the output neuron activity, τout is the time constant, Wout and bout are the weights and bias, and rrec concatenates membrane potentials and spike activities from recurrent neurons.

Resting-state whole-brain in vivo imaging data
Our Drosophila whole-brain in vivo imaging data derive from experimental studies by Mann et al. (2017)17 and Turner et al. (2021)18. They recorded resting-state neural activity in Drosophila using whole-brain calcium imaging of adult flies expressing the calcium-sensitive fluorescent protein GCaMP6s and membrane-tagged tdTomato. The central brain was exposed and imaged in vivo using resonant scanning two-photon microscopy at 3 μm isotropic resolution with 1.2 Hz temporal sampling for approximately 17 minutes per recording. They aligned functional data to a standard brain atlas using CMTK computational tools and defined regions of interest (ROIs) based on anatomical boundaries from established atlases (Ito or Branson). From each ROI, they extracted fluorescence time series, applied high-pass filtering to remove slow drift, and converted signals to ΔF/F measurements. The exemplary traces of recording ΔF/F can be obtained in Fig. 6A. We used deconvolution methods to transform these calcium imaging data into the neuropil firing rate (Supplementary Note H).

Whole-brain Drosophila connectome-constrained model
Network architecture. The whole-brain spiking network contains over 125,000 neurons and 50 million synaptic connections, organized into interconnected components (Supplementary Fig. 18).

Training operates on two distinct timescales. The neuropil firing rate is sampled at 1.2 Hz, with each datapoint representing approximately 833 ms of biological time. We denote the neuropil firing rate at timestep T as NeuropilFRT. In contrast, the spiking network simulation uses a 0.2 ms timestep, yielding 4165 simulation steps per neuropil sample. We use t for simulation time and T for neuropil time.

At each simulation step t, the network processes both recurrent input It (see Recurrent encoder) and external input It derived from the previous neuropil firing rate (see Input encoder). The recurrent connectivity matrix W is derived from the FlyWire connectome15. After completing all simulation steps within a neuropil time window, we compute each neuron’s firing rate and transform it to neuropil-level firing rates for comparison with experimental data (see Neuropil firing rate readout).

Network dynamics. The network dynamics follow the LIF model developed by Shiu et al.16:

τmbr dv i/dt = -g i(v i - Vresting) + Iext, (29)

dg i/dt = -g i/τrec + Irec, (30)

where v i is the membrane potential and g i is the synaptic conductance of neuron i. When v i ≥ Vthreshold, the neuron fires and resets to Vreset, followed by a refractory period of Trefractory = 2.2 ms.

Biophysical parameters are set to physiologically realistic values: Vresting = Vreset = -52 mV, Vthreshold = -45 mV, Rmbr = 10 KΩ ⋅ cm2, Cmbr = 2 μF ⋅ cm-2, Tmbr = 20 ms, τ = 5 ms, Tdelay = 1.8 ms, and Wsyn = 0.275 mV.

Input encoder. The input encoder transforms neuropil firing rates from the previous time step into background input for the current time step. Specifically, the encoder takes the neuropil firing rate at timestep T-1 (denoted as NeuropilFRT-1) and generates the external input IT ext for each neuron at timestep T. All firing rates are measured in Hz. During the T-th neuropil timestep, neuron i receives Poisson spike noise based on the transformed firing rate FRT ext,i. To capture temporal dependencies between different time steps, we employ a gated recurrent unit (GRU)55 followed by a linear transformation as the encoder for the Drosophila whole brain model. The transformation is given by:

FRT ext = Linear(GRU(NeuropilFRT-1)) (31)

Then, the neuropil-time-scale firing rate drives Poisson spike noise added to the membrane potential:

Iext(t) = ∑j δ(t - tj) wi, (32)

where δ(t-tj) represents Poisson spike times and wi is the learned synaptic weight for neuron i. During training, NeuropilFRT corresponds to experimental recordings; during inference, it refers to the simulated firing rate at the last step.

Recurrent encoder. The recurrent input to neuron i at time t is:

Irec,i(t) = ∑j δ(t - tspk j) wj,i (33)

The synaptic weight wj,i is computed by scaling the FlyWire connectome strength15 with a sign factor (+1 for excitatory, -1 for inhibitory) and the baseline synaptic weight Wsyn.

Neuropil firing rate readout. Each neuropil’s firing rate is computed as a weighted average of its constituent neurons’ firing rates. Because neuron density and connectivity vary across neuropils, simple averaging would distort the measured activity. Instead, each neuron is weighted by its normalized synaptic projection count into that neuropil:

SimFRi = ∑j rj wj,i / ∑j wj,i (34)

where rj is the firing rate of neuron j and wj,i is its normalized synaptic projection weight to neuropil i, obtained from the FlyWire connectome15. This weighting emphasizes neurons with stronger anatomical influence, yielding a more biologically faithful estimate of neuropil activity.

Robust automatic online learning through the Jaxpr compilation
The backpropagation algorithm owes much of its success to automatic differentiation frameworks like PyTorch21 and TensorFlow22, which abstract away gradient computation details. Inspired by this, we developed a general programming interface for online learning in SNNs that automatically derives eligibility traces and learning signals. Users need only define network dynamics and architecture; the interface handles all online learning computations.

Developing this framework requires addressing three challenges: (1) deriving eligibility traces, (2) generating learning signals, and (3) computing parameter gradients. We leverage Jaxpr, JAX’s intermediate representation56, to inspect computation graphs using abstract symbols. This enables us to analyze user-defined models and regenerate code that encompasses both state evolution and online learning, a process we call “online learning compilation”. During online learning compilation, we can also perform comprehensive optimizations to maximize computational performance and minimize the memory footprint that we need (see Supplementary Note I).

Deriving eligibility traces. Automatic derivation of eligibility traces involves three steps.

First, we analyze hidden state groups h l = [v1, ..., vd] within each recurrent layer to determine how state variables interact. This yields the block diagonal matrices Dt l for each layer l ∈ {1,…,L}.

Second, we analyze weight-to-hidden relationships ∂ht/∂θt by enumerating the computation graph to determine which hidden state each parameter connects to. This yields the block diagonal matrix Dt f.

Third, we generate the Jaxpr representation of eligibility traces (Eqs. (4)–(8)) based on the above analysis.

Generating learning signals. We use backpropagation through deep layers to generate learning signals (Eq. (2)).

Computing Eq. (25) need to save the gradients of intermediate variables. This is not allowed directly in JAX since hidden states computed within a function cannot be passed as function parameters. To overcome this, we employ a perturbation trick: we add zero perturbation δ to each hidden state, yielding ht + δ. By treating this perturbation as a function parameter, we exploit the equivalence ∂Lt/∂ht = ∂Lt/∂(ht + δ). We then use JAX’s jax.vjp interface to compute the ∂Lt/∂(ht + δ) derivatives and extract the result as a Jaxpr representation.

For computing Eq. (26), we directly take derivatives of hidden states at the previous time step ∂Lt/∂ht-1.

Computing parameter gradients. We integrate eligibility traces and learning signals using JAX’s jax.custom_vjp to define a custom vector-Jacobian product rule. This enables weight gradients to be accessed via standard interfaces like jax.grad while internally applying our online learning rules.

Our implementation defines two functions: (1) a forward pass function that computes network propagation and eligibility traces, and (2) a backward pass function that computes learning signals and parameter gradients by combining eligibility traces with learning signals. See Listing S4 for an example.

The resulting interface requires only two lines of code:

# define online computation using the pp-prop algorithm
model = braintrace.pp_prop(net, decay_or_rank=0.98)
# define online computation using the D-RTRL algorithm
model = braintrace.D_RTRL(net)
# compile the online learning graph using the input data
model.compile_graph(inputs)

Peer review rounds
[{"round": 1, "reviewer_comments": {"reviewer_1": "In their manuscript “BrainScale: Enabling Scalable Online Learning in Spiking Neural Networks” Wang et al. present a\nsubstantive amount of work with three main contributions:\n• Derivation and implementation of another (two variants of) approximation of RTRL in spiking neural networks\n• A software platform for SNN simulation and training built on JAX, including methods for using their new RTRL\napproximations\n• Validation of both the learning rules and the software on computational neuroscience models and neuromorphic computing\nbenchmarks\nOverall, the paper is very well structured and written and the results look sensible. However, I have some concerns\nMajor:\n1. In the validation on neuromorphic benchmarks, the authors are not relating their results to the recent literature. To name\njust a few works,\na. On SHD:\ni. Hammouamri, I., Khalfaoui-Hassani, I. and Masquelier, T., 2023. Learning delays in spiking neural networks using dilated\nconvolutions with learnable spacings. arXiv preprint arXiv:2306.17670.\nii. Bittar, A. and Garner, P.N., 2022. A surrogate gradient spiking baseline for speech command recognition. Frontiers in\nNeuroscience, 16, p.865897\niii. Nowotny T, Turner JP, Knight JC. Loss shaping enhances exact gradient learning with eventprop in spiking neural\nnetworks. arXiv preprint arXiv:2212.01232. 2022 Dec 2.\niv. Sun, P., Chua, Y., Devos, P. and Botteldooren, D., 2023. Learnable axonal delay in spiking neural networks improves\nspoken word recognition. Frontiers in Neuroscience, 17, p.1275944.\nv. Pengfei Sun, Ehsan Eqlimi, Yansong Chua, Paul Devos, Dick Botteldooren, Adaptive Axonal Delays in feedforward\nspiking neural networks for accurate spoken word recognition, arXiv:2302.08607\nvi. Yu, C., Gu, Z., Li, D., Wang, G., Wang, A. and Li, E., 2022. Stsc-snn: Spatio-temporal synaptic connection with temporal\nconvolution and attention for spiking neural networks. Frontiers in Neuroscience, 16, p.1079357.\nvii. Yao, M., Gao, H., Zhao, G., Wang, D., Lin, Y., Yang, Z. and Li, G., 2021. Temporal-wise attention spiking neural networks\nfor event streams classification. In Proceedings of the IEEE/CVF International Conference on Computer Vision (pp. 10221-\n10230).\nviii. Patiño-Saucedo, A., Yousefzadeh, A., Tang, G., Corradi, F., Linares-Barranco, B. and Sifalakis, M., 2023, May. Empirical\nstudy on the efficiency of spiking neural networks with axonal delays, and algorithm-hardware benchmarking. In 2023 IEEE\nInternational Symposium on Circuits and Systems (ISCAS) (pp. 1-5). IEEE.\nix. Yin, B., Corradi, F. and Bohté, S.M., 2021. Accurate and efficient time-domain classification with adaptive spiking\nrecurrent neural networks. Nature Machine Intelligence, 3(10), pp.905-913.\nAll these achieve >90% accuracy much beyond the numbers in Table 1.\nb. On DVS Gesture:\ni. She, X., Dash, S. and Mukhopadhyay, S., 2021. Sequence approximation using feedforward spiking neural network for\nspatiotemporal learning: Theory and optimization methods. In International Conference on Learning Representations.\nii. Subramoney, A., Nazeer, K.K., Schöne, M., Mayr, C. and Kappel, D., 2022. Efficient recurrent architectures through activity\nsparsity and sparse back-propagation through time. arXiv preprint arXiv:2206.06178.\niii. Apolinario, M.P.E. and Roy, K., 2023. S-tllr: Stdp-inspired temporal local learning rule for spiking neural networks. arXiv\npreprint arXiv:2306.15220.\niv. Sun, H., Cai, W., Yang, B., Cui, Y., Xia, Y., Yao, D. and Guo, D., 2023. A synapse-threshold synergistic learning approach\nfor spiking neural networks. IEEE Transactions on Cognitive and Developmental Systems, 16(2), pp.544-558.\nv. Xiao, M., Meng, Q., Zhang, Z., He, D. and Lin, Z., 2022. Online training through time for spiking neural networks. Advances\nin neural information processing systems, 35, pp.20717-20730.\nvi. Samadzadeh, A., Far, F.S.T., Javadi, A., Nickabadi, A. and Chehreghani, M.H., 2023. Convolutional spiking neural\nnetworks for spatio-temporal feature extraction. Neural Processing Letters, 55(6), pp.6979-6995.\nvii. She, X., Dash, S. and Mukhopadhyay, S., 2021. Sequence approximation using feedforward spiking neural network for\nspatiotemporal learning: Theory and optimization methods. In International Conference on Learning Representations.\nAll of these outperform what is presented in Table 1.\nThe authors must engage with the recent literature if they want to make credible claims about benchmarking accuracy.\nBased on what is presented in Table 1 and the papers referenced above, I believe the text: “Table 1 presents a summary of\nour model’s training results, juxtaposed with existing state-of-the-art\nonline and offline SNNs. Our findings revealed that the models trained with our online learning algorithms consistently and\nsignificantly outperform SNNs trained using offline BPTT approaches and other online methods,” is not backed up by\nevidence. Similarly in the abstract, “as well as achieves comparable or even superior performance in neuromorphic tasks.”\ndoes not reflect reality.\n2. Some parts of the paper are not necessary, including Figure 1 and Figure 2, which has the additional problem that the\nunderlying concepts (here called “alignPre” and “alignPost”) are well-known in the simulator literature and are used across\nmodern simulator software (e.g. Brian, GeNN, NEST, …). This needs to be said and relevant works need to be cited.\nSimilarly, I would recommend that the Methods section is reviewed and all parts that do not play a direct role in this work are\nremoved. In particular, the paragraph about dendritic computation (and the long list of references with it) seems completely\nout of place. Also, many parts of Figure 4 are more of an advert than a scientific illustration/figure.\nOther concerns:\n1. The name “BrainScale” is very unlucky as it is almost identical with the well-established BrainScaleS neuromorphic\nhardware which is well-known worldwide and has existed for almost 20 years. This clash should be avoided.\n2. Line 83: “Eq. 1 is capable of modeling the exact dynamics of realistic brain networks under …” is unclear. Equation 1 is a\ndiscrete time approximation of a class of neuron models that is usually understood to be defined through ODEs. In what\nsense is it modelling the exact dynamics of brain networks? Besides, in what sense are the dynamics of brain networks\nexact?\n3. Line 90ff : The concepts of “AlignPre” and “AlignPost” are not new and should be referenced to the literature\n4. Line 200ff: “To facilitate seamless integration of online gradient computations with standard AD routines, we have made\nour online algorithms compatible with reverse-mode AD in JAX by customizing the vector-Jacobian product rule. This\nadaptation enables the online computation process to closely resemble that of traditional deep learning models, requiring no\nmore than two lines of additional code (the red box in Figure 4F).” – unclear, how can a forward algorithm be made\n“compatible with reverse-mode AD”?\n5. Line 251: “However, without surrogate gradients, BPTT failed to train successfully on this task, with its training accuracy\nfluctuating around the 50% random baseline (Figure 5F). ” – unclear: how do you train using BPTT without a surrogate\ngradient? What are you using for the non-existing derivative of the Heaviside function? Is that a sensible thing to try?\n6. Line 267: “The capability of our online learning algorithms to learn long-term dependencies raises important questions.\nThe first is how our algorithms achieve long-term dependency learning.” – awkward.\n7. Caption for Figure 6 E: “Sorted raster plot” – sorted according to what criterion?\n8. Line 381: “Spiking networks have recently attracted substantial interest in the field of BIC due to their superior capability to\nreplicate the structure, function, and dynamics of the human brain over traditional artificial neural networks …” – a) define the\nacronym “BIC”; b) I don’t think it’s accurate to say that neuromorphic computing was interested in SNN simply because they\nare more similar to brains. (let alone the human brain!)\n9. Line 470: “learn and adapt in real-time, ultimately bringing digital whole-brain simulations closer to reality.” – there is\nnothing in the paper to suggest we are close to whole-brain simulations. This sentence should be deleted.\n10. Line 490/ equation (8): dv/dt incorrectly used for a discrete-time recurrence equation.\n(Remarks on code availability)\nThe code is not actually available I believe. The manuscript states only that it will be available at certain URLs but I have\ntried them and they are not yet accessible.", "reviewer_2": "The article presents a scalable algorithm for brain dynamics training. The authors defined three postulates regarding\nscalability in relation to wide applicability, effective memory management and prototyping using target hardware\narchitectures. The authors performed tests on a relatively narrow set of LIF and ALIF neurons and various synapse\nconfigurations described in Method 7.2. I positively evaluate narrowing the tests to such a set of models due to their use in\nmodern neuroprocessors. Additionally, the use of a sigma-delta network is worth considering in future analyses. Both the hot\ntopic of training spiking neural networks and the solutions proposed by the authors make the work valuable. I positively\nevaluate the division of work proposed by the authors. The addition of Chapter 7 “Methods” and addition of “Supplementary\nmaterials” make it much easier to navigate through the extensive material. The work is written carefully and is a reliable\nstudy of the proposed approach. However, I think that the work lacks a reliable comparison of the results achieved by the\nauthors with the results of existing frameworks. The functionality of frameworks would also require comparison with existing\ntools. Detailed comments below:\n1. I suggest organizing the numbering of the bibliography, e.g. according to the order of use.\n2. Please verify the content of the Supplementary materials section. Not all of its elements have links in the text, e.g. there\nare no links (and therefore comments) to tables S1-S11. I suggest additionally increasing the font size on some charts in this\nsection, e.g. S12-S14. Please verify the notation in tables S3, S4, S8, S9, S10, S11 with the values V_v, V_rest, V_reset,\nV_th (should be as in table S7).\n3. I suggest adding information about the access date to the used datasets in the Data availability section.\n4. I suggest adding a legend to Figure 2 explaining the meaning of the colors and symbols used in the diagrams. In the case\nof Fig. 2A and 2B, I suggest introducing a nomenclature for the columns (before and after the merging process)\n5. I suggest verifying that all shortcuts are defined on first use. For example, DMTS appears in the caption of Figure 5 on\npage 9 (linked to on page 8), but its explanation appears on page 10. The abbreviation BIC is not explained at all.\n6. I would suggest describing how the online learning program compiler illustrated in Fig. 4 corresponds to the general\nconcept illustrated in Fig. 1.\n7. One of the most recognized SNN learning algorithms is the ANN-SNN conversion (e.g. bootstrap algorithm). It turns out to\nbe competitive with many new and innovative approaches. On the other hand, conversion is often worse according to the\nmemory usage criterion. I believe that the summary of the BPTT, D-RTRL, ES-D-RTRL algorithms shown in Fig. 5JKL (data\nlength vs. memory/time per bach/accuracy) should be compared with the results for the conversion method. It may be\nparticularly valuable to compare the time per batch parameter, which is usually incomparably lower for conversion than for\nbackpropagation. Similarly, I would like to ask to supplement the data in Table 1 with the parameters obtained for the ANN-\nSNN conversion. I would suggest bolding the solution with the best result in the list. It is also worth comparing with the\nconversion algorithm the results of the accuracy for cognitive task performance analysis demonstrated in Fig. 6BC.\n8. The authors describe the compilation functionality for CPU, GPU and TPU. I believe that the results of the analyzes in Fig.\n5JK should be presented separately for all three architectures. It is worth considering presenting such an analysis for the\naccuracy parameter (Fig. 5L) - there are slight discrepancies in this parameter in frameworks depending on CPU vs. GPU.\n9. Please provide information about exactly which devices (CPU, GPU, TPU) the tests were carried out for. The \"time per\nbatch\" parameter is difficult to assess without detailed hardware data.\n10. The work advertises the BrainScale tool as a framework with extensive functionality. The work lacks a comparison in\nterms of functionality of this tool with other existing ones, e.g. Sinabs, Spyketorch, Brain2, Lava-dl, Norse. To begin with, I\nsuggest comparing these frameworks in terms of functionality, e.g. available models and algorithms. It would also be\nvaluable for potential users to compare the performance of the selected dataset (e.g. accuracy and time per batch).\n11. As for the scalability postulates proposed by the authors, there are no hard arguments to prove their fulfillment by\ncomparing BrainScale's performance with existing solutions. And so:\na) in relation to training scope - the statement \"BrainScale significantly reduces model complexity\" would be worth\nsupporting with a comparison of model complexity for models implemented using a conversion algorithm, which usually\nproposes models of very low complexity\nb) in relation to memory efficiency - the statement \"properties that can reduce online learning computations\" would be worth\nsupporting with a comparison of e.g. the avarage or peak memory usage parameter with another algorithm e.g. STDP, which\nfor a small batch size may require relatively small memory resources\nc) with regard to implementation - the statement \"facilitating wide adoptions of SNNs on various applications\" should be\ncarried out through analyzes for various architectures, i.e. CPU, GPU and TPU", "reviewer_3": "Wang et al. introduce BrainScale, a framework for automatic differentiation for spiking networks, that internally uses an\nalgorithm called ES-D-RTRL to compute gradients for spiking networks. ES-D-RTRL is an application of approximate\ndiagonalized real-time recurrent learning (RTRL) to spiking neural networks. Wang et al. present some first results of this\nnew framework to demonstrate its functionality. The manuscript contains a number of nice concepts and is presented clearly.\nHowever, BrainScale lacks novelty over existing methods and should discuss prior work more exhaustively.\nThe need for software tools to automate the process of constructing large and complex spiking networks and suitable\nlearning mechanisms has been extensively explored in prior work, e.g. early on by Chris Eliasmith's Nengo framework [7],\nEshraghian et al., [2] or jaxsnn, Müller et al. 2024 [3], to name just a few (see e.g. Pietrzak et al. 2023 [1] for a review).\nClearly, BrainScales solves many problems in a different way, but in the light of this extensive prior work, the introduction\nand discussion overstate the novelty of the proposed framework. Relevant related work should be cited alongside a\nclarification of what is novel in BrainScale.\nThe idea of using RTRL for spiking neurons is not new either. Note that even e-Prop (ref. 6) is a simple variant of spiking\nRTRL. Mingqing Xiao et al. had a paper in NeurIPS two years ago that is conceptually very similar to the work presented\nhere. It is also cited in the abstract of the manuscript (ref. 82). They used a leaky-integrate-and-fire neuron model and derived\nan approximate version of RTRL. The approximation used a very similar diagonalization to achieve linear scaling (as done\nin the current manuscript). They show theoretical and experimental results for recurrent and feedforward networks. LIF\nneurons are very similar to the exponential synaptic response kernel that is used in the derivation of ES-D-RTRL. In the light\nof this similarity, the manuscript seems a bit unfair towards ref. 82. It is only cited in the abstract in the statement \"... or use\noversimplified models [82], failing to capture the intricate brain dynamics.\" What exactly are the intricate brain dynamics that\nBrainScale can capture and prior methods failed to address? A detailed side-by-side comparison would be in need here.\nAlso the simulation results reported here are misleading in the sense that they miss a fair comparison with state-of-the art.\n94% test accuracy on DVS gesture dataset is not state of the art. Modern machine learning method achieve 100% on this\ndata set, see: [4]. Why this gap? Should be discussed. Furthermore, Xiao et al. [ref 82] also reported superior results on DVS\ngesture dataset (>96%). Please include and discuss their result and compare the difference in detail, given the relevance of\nthis work (see also comment above). They also presented compelling results for the more complex cifar 10 and 100 data\nsets, among others. Would be interesting to include those as well here. Also on Spiking speech commands much better\nresults have been reported, e.g. [5, 6]. Why this gap? Please discuss. Furthermore, the comparison with BPTT without\nsurrogate is not very meaningful, since BPTT for spiking neurons in that case simply doesn't work. Also, why stop learning in\nFig.5E,F at 95%? What is significant about this number, and how would convergence behave after reaching that point?\nReferences:\n[1] Pietrzak et al. 2023: https://www.mdpi.com/1424-8220/23/6/3037\n[2] Eshraghian https://ieeexplore.ieee.org/iel7/5/10251702/10242251.pdf\n[3] Müller et al. 2024. https://arxiv.org/pdf/2401.16841\n[4] Schöne et al. 2024A. https://arxiv.org/pdf/2411.12603\n[5] Schöne et al. 2024B. https://arxiv.org/pdf/2404.18508\n[6] Soydan et al. 2024. https://arxiv.org/abs/2410.03464\n[7] https://www.frontiersin.org/journals/neuroinformatics/articles/10.3389/fninf.2013.00048/full\n(Remarks on code availability)\nCode will be published after acceptance. The provided links to source code are not active yet."}, "author_response": "Dear reviewers,\nThank you for considering our manuscript titled \"BrainScale, Enabling Scalable Online\nLearning in Spiking Neural Networks \" for publication in Nature Communications . We\nare grateful for the criticism, suggestions, and feedbacks provided by the reviewers.\nWe found that these comments have been invaluable in enhancing the quality and scope\nof our work. Over the past eight months, we have diligently made substantial\nrevisions to both the manuscript , figures , software , applications , and\nevaluations , thoroughly addressing all comments and concerns raised by the\nreviewers. Below, we present a detailed, point-by-point response to each of the\nreviewers’ comments.\nResponse to Reviewer #1\nWe appreciate Reviewer #1's positive feedback on the motivation, design, and\ncontribution of BrainScale . We have addressed the following concerns raised to\nstrengthen the manuscript.\nComment 1:\nIn the validation on neuromorphic benchmarks, the authors are not relating their\nresults to the recent literature.\nOn SHD: All these achieve >90% accuracy much beyond the numbers in Table 1.\n1. i. Hammouamri, I., Khalfaoui-Hassani, I. and Masquelier, T., 2023. Learning\ndelays in spiking neural networks using dilated convolutions with learnable\nspacings. arXiv preprint arXiv:2306.17670.\n2. ii. Bittar, A. and Garner, P.N., 2022. A surrogate gradient spiking baseline\nfor speech command recognition. Frontiers in Neuroscience, 16, p.865897\n3. iii. Nowotny T, Turner JP, Knight JC. Loss shaping enhances exact gradient\nlearning with eventprop in spiking neural networks. arXiv preprint\narXiv:2212.01232. 2022 Dec 2.\n4. iv. Sun, P., Chua, Y., Devos, P. and Botteldooren, D., 2023. Learnable axonal\ndelay in spiking neural networks improves spoken word recognition. Frontiers\nin Neuroscience, 17, p.1275944.\n5. v. Pengfei Sun, Ehsan Eqlimi, Yansong Chua, Paul Devos, Dick Botteldooren,\nAdaptive Axonal Delays in feedforward spiking neural networks for accurate\nspoken word recognition, arXiv:2302.08607\n6. vi. Yu, C., Gu, Z., Li, D., Wang, G., Wang, A. and Li, E., 2022. Stsc-snn:\nSpatio-temporal synaptic connection with temporal convolution and attention\nfor spiking neural networks. Frontiers in Neuroscience, 16, p.1079357.\n7. vii. Yao, M., Gao, H., Zhao, G., Wang, D., Lin, Y., Yang, Z. and Li, G., 2021.\nTemporal-wise attention spiking neural networks for event streams\nclassification. In Proceedings of the IEEE/CVF International Conference on\nComputer Vision (pp. 10221-10230).\nNo. 1 / 22\n8. viii. Patiño-Saucedo, A., Yousefzadeh, A., Tang, G., Corradi, F., Linares-\nBarranco, B. and Sifalakis, M., 2023, May. Empirical study on the efficiency\nof spiking neural networks with axonal delays, and algorithm-hardware\nbenchmarking. In 2023 IEEE International Symposium on Circuits and Systems\n(ISCAS) (pp. 1-5). IEEE.\n9. ix. Yin, B., Corradi, F. and Bohté, S.M., 2021. Accurate and efficient time-\ndomain classification with adaptive spiking recurrent neural networks. Nature\nMachine Intelligence, 3(10), pp.905-913.\nOn DVS Gesture: All of these outperform what is presented in Table 1.\n1. i. She, X., Dash, S. and Mukhopadhyay, S., 2021. Sequence approximation using\nfeedforward spiking neural network for spatiotemporal learning: Theory and\noptimization methods. In International Conference on Learning Representations.\n2. ii. Subramoney, A., Nazeer, K.K., Schöne, M., Mayr, C. and Kappel, D., 2022.\nEfficient recurrent architectures through activity sparsity and sparse back-\npropagation through time. arXiv preprint arXiv:2206.06178.\n3. iii. Apolinario, M.P.E. and Roy, K., 2023. S-tllr: Stdp-inspired temporal\nlocal learning rule for spiking neural networks. arXiv preprint\narXiv:2306.15220.\n4. iv. Sun, H., Cai, W., Yang, B., Cui, Y., Xia, Y., Yao, D. and Guo, D., 2023. A\nsynapse-threshold synergistic learning approach for spiking neural networks.\nIEEE Transactions on Cognitive and Developmental Systems, 16(2), pp.544-558.\n5. v. Xiao, M., Meng, Q., Zhang, Z., He, D. and Lin, Z., 2022. Online training\nthrough time for spiking neural networks. Advances in neural information\nprocessing systems, 35, pp.20717-20730.\n6. vi. Samadzadeh, A., Far, F.S.T., Javadi, A., Nickabadi, A. and Chehreghani,\nM.H., 2023. Convolutional spiking neural networks for spatio-temporal feature\nextraction. Neural Processing Letters, 55(6), pp.6979-6995.\n7. vii. She, X., Dash, S. and Mukhopadhyay, S., 2021. Sequence approximation\nusing feedforward spiking neural network for spatiotemporal learning: Theory\nand optimization methods. In International Conference on Learning\nRepresentations.\nThe authors must engage with the recent literature if they want to make credible\nclaims about benchmarking accuracy. Based on what is presented in Table 1 and the\npapers referenced above, I believe the text: “Table 1 presents a summary of our\nmodel’s training results, juxtaposed with existing state-of-the-art online and\noffline SNNs. Our findings revealed that the models trained with our online\nlearning algorithms consistently and significantly outperform SNNs trained using\noffline BPTT approaches and other online methods,” is not backed up by\nevidence. Similarly in the abstract, “as well as achieves comparable or even superior\nperformance in neuromorphic tasks.” does not reflect reality.\nThanks for the valuable and detailed comment.\nIn the revised manuscript, we have reported strong performance on neuromorphic\ndatasets such as SHD and DVS Gesture, and have made the following improvements to\naddress your concerns:\n1. Controlled comparisons under unified conditions\nNo. 2 / 22\nWe performed controlled experiments using a standardized setting across models,\ndatasets, and training hyperparameters, to isolate the effect of the learning\nalgorithms. Specifically, we re-implemented four high-performing spiking neuron\nmodels (LIF, RLIF, adLIF, RadLIF) for SHD dataset, including those used in [Bittar &\nGarner, 2022], under a unified training framework. However, we intentionally removed\ncertain components such as Batch Normalization that are incompatible with the\nconstraints of online learning. We trained them under strictly controlled conditions\n- same hyperparameters, same data preprocessing, same hardware. Under these strictly\ncontrolled and fair conditions, our ES-D-RTRL algorithm achieves training performance\ncomparable to BPTT, and significantly outperforms existing online learning algorithms\n(see Table 1).\nMoreover, following the suggestion of Reviewer #3, we re-implemented the EGRU model\nfrom Subramoney et al. (ICLR, 2023) and evaluated it on the DVS gesture dataset. Our\nES-D-RTRL algorithm achieves over 97% accuracy on this benchmark, matching the\ncurrent state-of-the-art performance (see Table 1).\nThese new results, we believe, more clearly reflect the intrinsic advantages of our\nlearning rule in a realistic and practical online learning context.\n2. Clarified wording in the manuscript\nWe have also carefully revised the wording in both the abstract and main text to\navoid overly strong claims. For instance, the sentence:\n“Table 1 presents a summary of our model’s training results, juxtaposed with\nexisting state-of-the-art online and offline SNNs. Our findings revealed that the\nmodels trained with our online learning algorithms consistently and significantly\noutperform SNNs trained using offline BPTT approaches and other online methods,”\nhas been revised to (Lines 266-268):\n“Table 1 presents the mean and standard deviation of classification accuracy over\nat least three independent runs for each dataset. The results indicate that ES-D-\nRTRL matches BPTT’s performance across most tested datasets and network\narchitectures.”\nSimilarly, the abstract now reads (Lines 15-18):\n“Extensive evaluations across diverse SNN dynamics and computational tasks\ndemonstrate that BrainScale consistently delivers strong training performance\nwhile maintaining extremely low memory usage and high computational efficiency.”\nWe acknowledge that some recent works, particularly those leveraging extensive\narchitectural optimizations, delayed learning mechanisms, and hybrid training\nstrategies, achieve higher absolute performance (often exceeding 90% on SHD or\nsignificantly outperforming on DVS Gesture). However, many of these advances, such as\ndilated convolutions, temporal attention, and adaptive delays, are are not directly\napplicable to strict online learning settings.\nTherefore, to ensure a fair and meaningful evaluation, we conducted an extensive\ncomparison in Table 1 against the latest state-of-the-art online learning methods\n(including OSTL, E-prop, ETLP, OSTTP, OTPE, S-TLLR, FPTT, and OTTT), deliberately\nexcluding offline approaches such as BPTT with many architectural optimizations. This\nallows us to better isolate and highlight the strengths of our proposed algorithm\nwithin the context of online learning, rather than conflating performance gains\narising from architectural enhancements or offline training techniques.\nNo. 3 / 22\nComment 2:\nSome parts of the paper are not necessary, including Figure 1 and Figure 2, which\nhas the additional problem that the underlying concepts (here called “alignPre”\nand “alignPost”) are well-known in the simulator literature and are used across\nmodern simulator software (e.g. Brian, GeNN, NEST, …). This needs to be said and\nrelevant works need to be cited. Similarly, I would recommend that the Methods\nsection is reviewed and all parts that do not play a direct role in this work are\nremoved. In particular, the paragraph about dendritic computation (and the long\nlist of references with it) seems completely out of place. Also, many parts of\nFigure 4 are more of an advert than a scientific illustration/figure.\nThanks for the suggestions.\nWe have carefully revised the manuscript in response to your comments, and made the\nfollowing specific changes:\n1. Removal of Figure 1\nWe agree that the original Figure 1 was not essential to the main scientific content\nof the paper. Its primary purpose was to help organize the narrative and clarify the\nmotivation and structure of the work. To improve focus and eliminate redundancy, we\nhave now removed this figure entirely.\n2. Revision of Figure 2 (now Figure 1) and clarification of AlignPre/AlignPost\nconcepts\nIndeed, the method of exponential synapses, referred to here as AlignPost , has been\nextensively discussed and implemented in various simulators such as BrainPy, Brian,\nGeNN, and NEST. However, the AlignPre abstraction has not been as clearly formalized\nor systematically explored in the existing literature.\nTo address this gap, the revised manuscript now explicitly acknowledges this prior\nwork and cites representative references accordingly. Nonetheless, our main\ncontribution lies in tackling a critical yet often overlooked issue in the field of\nonline learning algorithms: the lack of a systematic summary and formal definition of\nthe model classes these algorithms support. For instance, methods such as OTTT are\nlimited to simple neuron models like LIF and rely on assumptions that do not\ngeneralize to more complex neural dynamics. While more general algorithms such as e-\nprop and OSTL exhibit some degree of model flexibility, their original presentations\ndo not clearly define the scope of supported network architectures and neuronal\ndynamics. As a result, these algorithms offer limited guarantees when applied to\nuntested or more intricate models, raising concerns about their generalizability and\nstability in practice.\nTo overcome this limitation, we propose a set of explicit abstractions  AlignPre and\nAlignPost   that systematically delineate the class of neural network models\nsupported by online learning rules. This abstraction not only encompasses the\ndominant spiking neuron models used today but also provides a formal foundation for\napplying online learning algorithms to a broader range of network structures and\ndynamic systems. The clear presentation of these model boundaries is a cornerstone of\nthe extensibility and generalizability of our framework, and it plays a foundational\nrole in advancing the application of online learning methods to complex neural\nmodeling and neuromorphic systems. Moreover, this formalism will offer valuable\nguidance for the future design and refinement of online learning algorithms.\n3. Removal of dendritic computation paragraph and associated references\nNo. 4 / 22\nWe agree that the section on dendritic computation was not directly relevant to the\nscope of this manuscript. This paragraph and the accompanying references have been\nremoved for better focus and conciseness.\n4. Reworking of Figure 4 (now Figure 3)\nThanks for the feedback that the previous version of Figure 4 was overly promotional\nin tone. In the revised manuscript, we have substantially reworked this figure into a\nmore technical and explanatory illustration (now Figure 3). The updated figure\nfocuses on the hierarchical compilation mechanism used in BrainScale for online\nlearning, and clearly delineates how our approach differs from traditional BPTT-based\npipelines. The new version emphasizes how the compilation structure and system-level\ninsight enables the automation of online learning implementation, rather than just\nadvertising the software framework.\nComment 3:\nThe name “BrainScale” is very unlucky as it is almost identical with the well-\nestablished BrainScaleS neuromorphic hardware which is well-known worldwide and\nhas existed for almost 20 years. This clash should be avoided.\nWe are indeed aware that the name BrainScale bears some resemblance to the existing\nhardware platform BrainScaleS (brain-inspired multi-scale computation in\nneuromorphic hybrid systems). This similarity has prompted several rounds of internal\ndiscussion and exploration of alternative naming options. However, after careful\nconsideration, we believe that BrainScale remains the most appropriate name to\ncapture the core vision and technical essence of our system.\nFirst, BrainScale is tightly coupled with the broader Brain Modeling Ecosystem (htt\nps:  brainmodeling.readthedocs.io/) we are building. This ecosystem includes tools\nsuch as BrainUnit, which provides physical unit support; BrainEvent, which implements\nevent-driven computing capabilities; BrainCell, which enables biologically detailed\nneuron modeling; and BrainState, which supports state-based compilation. BrainScale\nplays a critical role in this ecosystem by delivering scalable online learning\ncapabilities, thereby completing the stack needed for advanced brain dynamics\nmodeling.\nSecond, the name BrainScale directly reflects the central objective of our system:\nenabling scalability in brain dynamics simulation and modeling during training. The\nterm concisely encapsulates our project's overarching mission, supporting the scaling\nup of brain-inspired computational models. Concepts such as \"scale up\" and\n\"scalability\" are not auxiliary features but are foundational to our work. As such,\nthe name BrainScale clearly and efficiently communicates this core innovation.\nMoreover, while the names are similar, BrainScale and BrainScaleS differ\nsignificantly in scope and technological focus. BrainScaleS is renowned for its\ncontributions to neuromorphic hardware. In contrast, BrainScale operates at the\nsoftware and algorithmic levels, representing a fundamentally different layer of the\ntechnology stack. This distinction highlights the divergence in application domains\nand intended usage between the two systems.\nThat said, we deeply value feedback from the community and will continue to assess\nthe appropriateness and potential impact of our naming choices to ensure they do not\ninadvertently cause confusion or undermine existing contributions.\nNo. 5 / 22\nComment 4:\nLine 83: “Eq. 1 is capable of modeling the exact dynamics of realistic brain\nnetworks under …” is unclear. Equation 1 is a discrete time approximation of a\nclass of neuron models that is usually understood to be defined through ODEs. In\nwhat sense is it modelling the exact dynamics of brain networks? Besides, in what\nsense are the dynamics of brain networks exact?\nThanks for the careful comment.\nWe agree that the original phrasing  “Eq. 1 is capable of modeling the exact dynamics\nof realistic brain networks under …”  was imprecise and potentially misleading.\nOur intended meaning was to highlight the consistency between simplified models (such\nas those based on the AlignPre and AlignPost abstractions) and the original, non-\nsimplified dynamical models. Specifically, we meant that the simplified formulation\nused in Eq. 1 preserves key dynamical features when compared to the original,\ncontinuous-time models, thereby enabling faithful simulation under certain\nassumptions.\nHowever, even with such consistency, Eq. 1 remains a discrete-time approximation and\ntherefore cannot be said to model the exact dynamics in the strict mathematical or\nbiological sense. A more accurate description is that, under a given time step and\nspecific numerical integration scheme, Eq. 1 provides a numerically faithful\napproximation to the solutions of differential equation-based neuron models. It\neffectively captures the essential dynamical behavior of these models, but only\nwithin the limits of its discretization.\nAlso, the term \"exact dynamics of brain networks\" is problematic in itself, as the\ndynamics of real brain networks are highly complex and inherently uncertain. No\nmathematical model can claim to capture them with absolute precision; all modeling\nefforts are necessarily performed at a particular level of abstraction.\nTo improve the rigor and clarity of the manuscript, we have removed the old phrase\nfrom the revised version. Instead, we now explicitly emphasize the rationale behind\nthe use of the simplified discrete form (see Lines 86-94) and discuss its advantages\nmore clearly (see Lines 95-99).\nComment 5:\nLine 90ff : The concepts of “AlignPre” and “AlignPost” are not new and should be\nreferenced to the literature\nIn the revised manuscript, we have added appropriate references to acknowledge this\nprior literature. Although AlignPost is widely recognized in the community through\nits use in exponential synapse models, to the best of our knowledge, our work is the\nfirst to provide a systematic formulation and unified abstraction of these two\nclasses of methods, AlignPre and AlignPost . We believe this structured presentation\nhelps clarify their conceptual distinctions and practical implications, and serves as\na foundation for more boundary-defined and extensible online learning systems.\nComment 6:\nNo. 6 / 22\nLine 200ff: “To facilitate seamless integration of online gradient computations\nwith standard AD routines, we have made our online algorithms compatible with\nreverse-mode AD in JAX by customizing the vector-Jacobian product rule. This\nadaptation enables the online computation process to closely resemble that of\ntraditional deep learning models, requiring no more than two lines of additional\ncode (the red box in Figure 4F).” – unclear, how can a forward algorithm be made\n“compatible with reverse-mode AD”?\nBy \"compatible with reverse-mode AD,\" we do not mean that our online gradient\nalgorithm (a forward-mode computation) internally uses reverse-mode automatic\ndifferentiation. Rather, we refer to the fact that we have customized the interface\nof our online learning algorithm so that it can be seamlessly invoked via standard\nreverse-mode AD tools, such as jax.grad . That is to say, by using familiar APIs\nlike:\ngrads = jax.grad(loss_fn)(params)\nusers can easily compute the forward gradients at time : . This compatibility\nsignificantly lowers the engineering overhead for integrating our methods into\nexisting deep learning pipelines.\nThis is precisely what we meant by saying that \"no more than two lines of additional\ncode\" are required (highlighted in red in Figure 4F in previous manuscript).\nWe have removed this wording from the main text and revised the corresponding\nparagraph to clarify the point and prevent potential misunderstandings (see\nSupplementary Note L.3).\nComment 7:\nLine 251: “However, without surrogate gradients, BPTT failed to train\nsuccessfully on this task, with its training accuracy fluctuating around the 50%\nrandom baseline (Figure 5F). ” – unclear: how do you train using BPTT without a\nsurrogate gradient? What are you using for the non-existing derivative of the\nHeaviside function? Is that a sensible thing to try?\nThanks for this question.\nIn our experiments on brain simulation models (Figure 4E–H, Figure 5), the outputs of\nthe recurrent layers are based on membrane potentials rather than spike events. This\ncontrasts with the brain-inspired computational models (Figure I–L, Table 1), where\nthe outputs are discrete spike trains. The main reason for this design choice is that\nspikes are extremely sparse in realistic neuronal dynamics, whereas membrane\npotentials provide richer and more continuous information, making them more suitable\nas decoding targets in these models. Moreover, since it remains unclear how the brain\ndecodes information, we opted to use membrane potential as the decoding target\n(Supplementary Note I.5). Nevertheless, decoding based on spikes can also\nsuccessfully accomplish this task, but converges slowly (data not shown).\nAlthough the decoding is based on membrane potentials, neuronal spikes play a\ncritical role in the internal forward dynamics, since spikes are still used to\nmediate synaptic interactions. Therefore, it is not a trivial thing that BPTT without\nsurrogate gradients cannot successfully train this tasks (as shown in Fig. 4F and\nFig. 5C).\nNo. 7 / 22\nWe have clarified this explanation in the revised manuscript to better articulate the\nrationale that allows this comparison to be made (see Supplementary Note I.5).\nComment 8:\nLine 267: “The capability of our online learning algorithms to learn long-term\ndependencies raises important questions. The first is how our algorithms achieve\nlong-term dependency learning.” – awkward.\nWe have revised the sentence to read:\n“We investigated how our algorithm achieves long-term dependency learning.”\nto improve clarity and readability.\nComment 9:\nCaption for Figure 6 E: “Sorted raster plot” – sorted according to what\ncriterion?\nThe spike raster plot in Figure 6E in the previous manuscript (now Figure 5E in the\nrevised manuscript) is sorted according to the timing of the first spike emitted by\neach neuron. Specifically, neurons that spike earlier are assigned lower indices in\nthe plot, making the temporal structure of the spiking activity more visually\ninterpretable.\nComment 10:\nLine 381: “Spiking networks have recently attracted substantial interest in the\nfield of BIC due to their superior capability to replicate the structure,\nfunction, and dynamics of the human brain over traditional artificial neural\nnetworks …” – a) define the acronym “BIC”; b) I don’t think it’s accurate to say\nthat neuromorphic computing was interested in SNN simply because they are more\nsimilar to brains. (let alone the human brain!)\nThanks for your constructive feedback.\na) Regarding the term \"BIC\", we originally used it as an acronym for Brain-Inspired\nComputing , but we agree it was a mistake not to define it upon first use. To improve\nclarity and consistency, we have removed all instances of \"BIC\" and replaced them\nwith the more standard term \"neuromorphic computing\" in the revised manuscript.\nb) We acknowledge that the previous wording was misleading. While the biological\nplausibility of spiking neural networks (SNNs) is one motivating factor, it is by no\nmeans the sole or primary reason for their adoption in neuromorphic computing. SNNs\noffer several practical advantages, including event-driven computation, temporal\ncoding, sparse activity, and high energy efficiency—especially beneficial in real-\ntime and resource-constrained scenarios.\nTo address this, we have removed the problematic sentence and replaced it with a more\nneutral and accurate statement (Line 256):\n“Spiking networks have been widely used in neuromorphic tasks [3, 12].”\nNo. 8 / 22\nThis revision avoids overgeneralization and better reflects the multi-faceted\nmotivations behind the use of SNNs. Thank you again for helping us improve the\naccuracy and rigor of our presentation.\nComment 11:\nLine 470: “learn and adapt in real-time, ultimately bringing digital whole-brain\nsimulations closer to reality.” – there is nothing in the paper to suggest we are\nclose to whole-brain simulations. This sentence should be deleted.\nWe have deleted this line, and totally revised the Discussion section.\nComment 12:\nLine 490/ equation (8): dv/dt incorrectly used for a discrete-time recurrence\nequation.\nWe have fixed this typo.\nComment 13:\nThe code is not actually available I believe. The manuscript states only that it\nwill be available at certain URLs but I have tried them and they are not yet\naccessible.\nWe have now made the full codebase publicly available at the URLs provided in the\nmanuscript. To ensure transparency and reproducibility, the repository includes all\nsource code, documentation, and scripts necessary to reproduce the main experiments\nand figures. We have also double-checked that the links are active and accessible\nwithout restrictions.\nNo. 9 / 22\nResponse to Reviewer #2\nWe thank the reviewer for the constructive feedback and revision suggestions. We are\nparticularly grateful for the positive assessment of the proposed framework and its\nfocus on scalability, the appreciation of our architectural modularity and model\nselection, and the recognition of the overall writing. Below we provide a detailed\nresponse to each point raised, and we have revised the manuscript accordingly.\nComment 1:\nI suggest organizing the numbering of the bibliography, e.g. according to the\norder of use.\nWe have organized the bibliography according to the order of citation.\nComment 2:\nPlease verify the content of the Supplementary materials section. Not all of its\nelements have links in the text, e.g. there are no links (and therefore comments)\nto tables S1-S11. I suggest additionally increasing the font size on some charts\nin this section, e.g. S12-S14. Please verify the notation in tables S3, S4, S8,\nS9, S10, S11 with the values V_v, V_rest, V_reset, V_th (should be as in table\nS7).\nThanks for your careful review regarding the Supplementary Notes section.\n1. Missing links to supplementary tables (S1–S11):\nThe primary purpose of these tables is to provide a comprehensive summary of model\nparameters used in the main text, serving as a reference for reproducing the\nexperimental results (although the source code can be a better way for\nreproducibility). While we acknowledge that the tables are not explicitly cited in\nthe main text, their role is to complement the source code in supporting\nreproducibility. Since one figure or table in the main text may correspond to\nmultiple supplementary tables, we chose not to add direct links, but rather to\npresent these tables collectively as parameter summaries.\n2. Chart readability (e.g., Figures S12–S14):\nWe agree that the font sizes in some supplementary figures are too small and may\naffect readability. We have revised Supplementary Figures by increasing the font size\nof axis labels, legends, and annotations to ensure better clarity, especially for\nprint and screen viewing. Please see Supplementary Fig. S9-S11.\n3. Consistency of notation in tables (S3, S4, S8–S11):\nWe have revised all Supplementary Tables to ensure that the notations for parameters\nsuch as , , , and are consistent with the standardized notation used\nin Table S7. Any inconsistencies in symbol formatting, subscripts, or definitions\nhave been corrected to maintain uniformity throughout the supplementary section.\nComment 3:\nNo. 10 / 22\nI suggest adding information about the access date to the used datasets in the\nData availability section.\nWe have added the access date.\nComment 4:\nI suggest adding a legend to Figure 2 explaining the meaning of the colors and\nsymbols used in the diagrams. In the case of Fig. 2A and 2B, I suggest\nintroducing a nomenclature for the columns (before and after the merging process)\nWe have revised Figure 2 (now presented as Figure 1 in the current manuscript) in\nresponse to your suggestions regarding the use of colors and symbols.\nThe original Figure 2A and 2B, which illustrated the theoretical motivation behind\nthe simplification of the AlignPre and AlignPost merging process, have been moved\nto the supplementary material as Supplementary Figure S6. To facilitate\ninterpretation, we have introduced a clear nomenclature for each column and added\nexplanatory labels in the figure legend to clarify the meaning of the different\nstages.\nThe original Figure 2C and 2D have been updated and are now presented as Figure 1A\nand 1B in the revised manuscript. These diagrams have been redrawn and revised to\nexplicitly align each visual component with the corresponding mathematical symbols\nin Eq. 1, including , , , and the parameter . Following your suggestion,\nwe have also added a legend to the figure and used color to distinguish the types\nof variables: presynaptic variables are shown in orange, postsynaptic variables in\nblue, and synaptic interactions in green. These changes aim to improve the clarity\nand interpretability of the figure.\nComment 5:\nI suggest verifying that all shortcuts are defined on first use. For example,\nDMTS appears in the caption of Figure 5 on page 9 (linked to on page 8), but its\nexplanation appears on page 10. The abbreviation BIC is not explained at all.\nWe have carefully reviewed the manuscript to ensure that all abbreviations and\nacronyms are defined at their first appearance. Specifically:\n\"DMTS (Delayed Match-to-Sample)\" is now defined at its first occurrence in the\nmain text, prior to its appearance in figure captions or references.\nThe term \"BIC (Brain-Inspired Computing)\" has been removed from the manuscript and\nreplaced with \"neuromorphic computing\", which is now used consistently throughout\nthe text.\nComment 6:\nI would suggest describing how the online learning program compiler illustrated\nin Fig. 4 corresponds to the general concept illustrated in Fig. 1.\nNo. 11 / 22\nIn the revised manuscript, the original Figure 1 has been removed. To address the\nconcern regarding the connection between the general concept and the compiler\nimplementation, we now provide a clearer and more detailed explanation of how our\nonline learning compiler operates. This is presented in the new section titled\n\"BrainScale facilitates automated online learning compilation\" , where we explicitly\ndescribe the compilation workflow and its underlying principles.\nParticularly, we describe how the compiler abstracts the user-define SNN models,\nextracts relationships between model parameters, connections, and states, and\ngenerates the code of online learning algorithms. We believe that this revised\npresentation can better convey the relationship between the conceptual framework and\nthe concrete implementation.\nComment 7:\nOne of the most recognized SNN learning algorithms is the ANN-SNN conversion\n(e.g. bootstrap algorithm). It turns out to be competitive with many new and\ninnovative approaches. On the other hand, conversion is often worse according to\nthe memory usage criterion. I believe that the summary of the BPTT, D-RTRL, ES-D-\nRTRL algorithms shown in Fig. 5JKL (data length vs. memory/time per\nbach/accuracy) should be compared with the results for the conversion method. It\nmay be particularly valuable to compare the time per batch parameter, which is\nusually incomparably lower for conversion than for backpropagation. Similarly, I\nwould like to ask to supplement the data in Table 1 with the parameters obtained\nfor the ANN-SNN conversion. I would suggest bolding the solution with the best\nresult in the list. It is also worth comparing with the conversion algorithm the\nresults of the accuracy for cognitive task performance analysis demonstrated in\nFig. 6BC.\nWe appreciate your insight regarding the potential comparison with ANN-SNN conversion\nmethods. However, we respectfully maintain our current experimental scope for the\nfollowing reasons:\nFirst, our team's expertise and focus lie in the area of online learning for SNNs,\nrather than offline training or ANN-SNN conversion methods. Incorporating and\nevaluating conversion-based baselines would require significant additional effort,\nincluding non-trivial reimplementation and tuning. We lack in-house expertise with\nANN-SNN conversion protocols, and reproducing these experiments would require a\nsubstantial investment of time and resources  far exceeding the eight months already\ndevoted to the present revision.\nSecond, our primary contribution lies in co-optimizing algorithms and frameworks for\nscalable online SNN learning. Therefore, we compare against state-of-the-art online\nlearning methods (see OSTL, E-prop, ETLP, OSTTP, OTPE, S-TLLR, FPTT, and OTTT in\nTable 1) that share the same learning constraints and objectives. Including ANN-SNN\nconversion would introduce a different class of methods that, while valuable,\naddresses distinct research questions related to knowledge transfer rather than\nonline adaptation.\nTo acknowledge this valuable suggestion, we have:\nAdded a paragraph in the “Limitations and Future Work” section to discuss how ANN-\nSNN conversion could complement our study. See Lines 413-419.\nBolded the best-performing solution in Table 1, as requested.\nNo. 12 / 22\nComment 8:\nThe authors describe the compilation functionality for CPU, GPU and TPU. I\nbelieve that the results of the analyzes in Fig. 5JK should be presented\nseparately for all three architectures. It is worth considering presenting such\nan analysis for the accuracy parameter (Fig. 5L) - there are slight discrepancies\nin this parameter in frameworks depending on CPU vs. GPU.\nIn the revised manuscript, we have conducted a comprehensive evaluation of our\nframework on CPU, GPU, and TPU architectures. Among these, the results on GPU are\npresented in Fig. 4J–L (corresponding to Fig. 5J–L in the previous manuscript). The\nresults on CPU and TPU are included in the supplementary material.\nSpecifically, Supplementary Fig. S13 presents the evaluation results obtained on a\nKaggle TPU-v3 instance, while Supplementary Fig. S14 shows the results from an Intel\nCPU server. The analyses across these two architectures are consistent with those\nshown in Fig. 4J–K: the ES-D-RTRL algorithm demonstrates the lowest memory consumption\nand fastest computation time, owing to its linear memory complexity.\nHowever, we regret that it is difficult for us to provide full training performance\nevaluations (e.g., accuracy metrics as in Fig. 4L) on CPU and TPU. Due to access\nlimitations, we rely on free Kaggle TPU-v3 instances, which provide only 20 hours of\nusage per month—insufficient for running the full-scale training benchmarks presented\nin Fig. 4L. Similarly, the extremely slow training speed on CPU makes it impractical\nto complete the extensive evaluations required for Fig. 4L. For these reasons, we did\nnot include accuracy comparisons across all three architectures.\nComment 9:\nPlease provide information about exactly which devices (CPU, GPU, TPU) the tests\nwere carried out for. The \"time per batch\" parameter is difficult to assess\nwithout detailed hardware data.\nWe have provided detailed hardware specifications used in our evaluations in Appendix\nO of the supplementary material.\nSpecifically, all experiments were conducted on a Linux server, with Intel Xeon W-\nKaggle-free TPU v3-8 cloud instance.\nComment 10:\nThe work advertises the BrainScale tool as a framework with extensive\nfunctionality. The work lacks a comparison in terms of functionality of this tool\nwith other existing ones, e.g. Sinabs, Spyketorch, Brain2, Lava-dl, Norse. To\nbegin with, I suggest comparing these frameworks in terms of functionality, e.g.\navailable models and algorithms. It would also be valuable for potential users to\ncompare the performance of the selected dataset (e.g. accuracy and time per batch).\nThanks for the feedback. In response, we have substantially revised the manuscript\nand added new comparisons and figures in Supplementary Note N to address this point.\nNo. 13 / 22\nWe would like to clarify that comparing BrainScale directly with other existing\nframeworks (e.g., Sinabs, Spyketorch, Brian2, Lava-dl, Norse) is not straightforward,\nas these tools differ significantly in their design goals, scopes, and underlying\nabstractions.\nFirst, BrainScale is not a standalone framework for spiking network modeling or\ntraining. Rather, it is a key component within a broader ecosystem we have developed.\nThis ecosystem is centered around BrainPy, which provides high-performance tools for\nSNN simulation, BPTT-based training, and automatic analysis of dynamical systems. In\nterms of simulation and training functionality, BrainPy is comparable to established\nframeworks such as Brian2, NEST, snnTorch, and SpikingJelly. BrainScale complements\nBrainPy by introducing a general and efficient online learning layer. Specifically,\nit enables the integration of biologically plausible and task-driven online learning\nalgorithms directly into arbitrary SNN models defined in BrainPy. This distinguishes\nour framework ecosystem from others, where online learning support is often either\nlimited in scope or tied to specific architectures.\nSecond, while some existing frameworks (e.g., Brian2, NEST) support online synaptic\nplasticity rules, these are typically limited to local rules (e.g., STDP) and are not\noptimized for large-scale cognitive tasks. In contrast, BrainScale supports general\nonline learning algorithms, including online variants of gradient-based methods. This\nenables scalable, efficient, and biologically inspired learning in networks with\ncomplex neuronal and synaptic dynamics.\nWe have added detailed functional comparisons (including available models, training\nalgorithms, and learning paradigms) between existing SNN frameworks (including\nNEURON, NEST, Nengo, jaxsnn, etc.) and our \"BrainPy + BrainScale\" system in\nSupplementary Note N . We hope this expanded comparison will provide valuable\nguidance for potential users and better position BrainScale within the current SNN\nmodeling landscape.\nComment 11:\nAs for the scalability postulates proposed by the authors, there are no hard\narguments to prove their fulfillment by comparing BrainScale's performance with\nexisting solutions. And so:\na) in relation to training scope - the statement \"BrainScale significantly\nreduces model complexity\" would be worth supporting with a comparison of model\ncomplexity for models implemented using a conversion algorithm, which usually\nproposes models of very low complexity\nb) in relation to memory efficiency - the statement \"properties that can\nreduce online learning computations\" would be worth supporting with a\ncomparison of e.g. the avarage or peak memory usage parameter with another\nalgorithm e.g. STDP, which for a small batch size may require relatively small\nmemory resources\nc) with regard to implementation - the statement \"facilitating wide adoptions\nof SNNs on various applications\" should be carried out through analyzes for\nvarious architectures, i.e. CPU, GPU and TPU\nWe have carefully addressed each point and clarified the relevant claims in the\nrevised manuscript. Below, we respond in detail to each of the concerns raised:\na) Training Scope\nNo. 14 / 22\nWe agree with the reviewer that the statement \"BrainScale significantly reduces model\ncomplexity\" was misleading. Our goal is not to reduce the complexity of models, but\nrather to support greater diversity in model structures, including complex neuronal\nand synaptic dynamics. We have revised the manuscript to remove this inaccurate\nphrasing.\nIn terms of training scope, BrainScale is designed to support a wide range of neuron\nand synapse types, including adaptive and non-linear dynamics (e.g., ALIF, GIF, NMDA-\nbased synapses), which are generally not supported by ANN-to-SNN conversion\napproaches. Conversion-based methods typically yield models composed of simple LIF\nneurons with low temporal resolution, and they cannot easily be extended to support\nricher dynamics or biologically plausible learning rules.\nIn contrast, BrainScale provides a general mechanism for online learning that is\ncompatible with a broad class of neuron and synapse models, abstracted using our\nproposed AlignPre and AlignPost interfaces. This enables the flexible training\nof models with complex dynamics  a key requirement for neuroscience modeling.\nWe have emphasized this distinction more clearly in the revised manuscript. Lines\n364-366.\nb) Memory Efficiency\nWe appreciate the suggestion to compare memory usage with STDP. As the reviewer\nnoted, STDP requires memory for pre- and post-synaptic traces, which is\ncomparable to the memory footprint of our ES-D-RTRL algorithm. However, the key\ndistinction lies in computational cost: STDP requires weight updates at every time\nstep, while ES-D-RTRL performs updates only when a learning signal is delivered. This\nleads to significantly lower computational overhead during online learning,\nespecially for long sequences or sparse feedback settings.\nFurthermore, unlike STDP, which is a purely local plasticity rule, ES-D-RTRL\napproximates gradient-based credit assignment. It achieves strong task performance\nwhile maintaining linear memory complexity  matching STDP in efficiency but with\nsignificantly greater generalization power. See supplementary Note M.\nc) Implementation and Portability\nThe core technical contribution of BrainScale lies in the automated compilation of\nonline learning rules into efficient, hardware-compatible implementations. To\nevaluate implementation scalability, we conducted experiments on diverse SNN\narchitectures running on CPUs, GPUs, and TPUs (Fig. 4, Fig. S13, and Fig. S14). Our\nframework is implemented in JAX, which enables seamless execution on all major\nhardware platforms via XLA compilation.\nWe hope these clarifications and additions address the reviewer's concerns and more\nclearly support our claims about BrainScale's scalability in terms of training scope,\nmemory efficiency, and implementation portability.\nNo. 15 / 22\nResponse to Reviewer #3\nWe thank Reviewer #3 for the feedback, as well as for recognizing the clarity of our\npresentation and the usefulness of our concepts. Below, we address the reviewer's\nconcerns in detail, particularly regarding the novelty of our work, comparison to\nprior art, and experimental validation.\nComment 1:\nWang et al. introduce BrainScale, a framework for automatic differentiation for\nspiking networks, that internally uses an algorithm called ES-D-RTRL to compute\ngradients for spiking networks. ES-D-RTRL is an application of approximate\ndiagonalized real-time recurrent learning (RTRL) to spiking neural networks. Wang\net al. present some first results of this new framework to demonstrate its\nfunctionality. The manuscript contains a number of nice concepts and is presented\nclearly. However, BrainScale lacks novelty over existing methods and should\ndiscuss prior work more exhaustively.\nThanks for the feedback and summary. To respond to the criticism of novelty, we have\nsignificantly expanded our discussion of prior work in the revised manuscript to\naddress your concerns. For a detailed comparison of frameworks , please refer to\nSupplementary Note N; for a comparison of algorithms , please see Supplementary Note\nM.\nWe disagree with the statement that \"BrainScale lacks novelty over existing\nmethods.\" As elaborated in the Introduction and Discussion sections, BrainScale\nintroduces innovation across four distinct levels:\n1. Model Universality: BrainScale is the first to systematically integrate the\nAlignPre and AlignPost mechanisms into a general-purpose online learning\nframework for spiking neural networks.\n2. Algorithm: We introduce ES-D-RTRL, the first algorithm to reduce the complexity of\nReal-Time Recurrent Learning (RTRL) to in the context of spiking networks\nby fully leveraging their event-driven nature. Please see Theorem C.3 for the\ntechnical details.\n3. Compiler Implementation: BrainScale is the first framework to provide an online\nlearning compiler that automates the transformation of high-level SNN models into\nefficient, low-level online learning implementations. This significantly lowers\nthe barrier to prototyping and deploying online learning algorithms.\n4. Application to Whole-Brain Modeling: To our knowledge, this is the first work to\napply an online-learning-based spiking framework to fit the whole-brain activity\nof a real animal  Drosophila   demonstrating the framework's scalability and\npractical utility.\nWe hope this revised manuscript and supplementary material more clearly convey the\ncontributions of BrainScale in relation to existing work.\nComment 2:\nNo. 16 / 22\nThe need for software tools to automate the process of constructing large and\ncomplex spiking networks and suitable learning mechanisms has been extensively\nexplored in prior work, e.g. early on by Chris Eliasmith's Nengo framework [7],\nEshraghian et al., [2] or jaxsnn, Müller et al. 2024 [3], to name just a few\n(see e.g. Pietrzak et al. 2023 [1] for a review). Clearly, BrainScales solves many\nproblems in a different way, but in the light of this extensive prior work, the\nintroduction and discussion overstate the novelty of the proposed framework.\nRelevant related work should be cited alongside a clarification of what is novel\nin BrainScale.\nThanks for the comment. We think there may be some misunderstanding regarding the\nfocus and scope of BrainScale .\nWe clarify the distinctions below and have also revised the Introduction and\nDiscussion sections accordingly to better articulate BrainScale's contributions\nrelative to prior work.\nFirst , BrainScale is not a simulation tool.\nBrainScale is not a framework for constructing or simulating spiking neural networks\nper se. That functionality is already provided by our previously developed BrainPy\nframework [Wang et al., eLife 2023], which offers general-purpose and efficient\nsimulation for large-scale brain dynamics. Instead, BrainScale builds on BrainPy to\nprovide  for the first time  a compiler-based system for automated online learning in\ngeneral SNN models. This is a fundamentally different goal from existing simulators\nor frontend frameworks.\nSecond , lack of support for automatic online learning in prior work.\nWe agree that tools like Nengo, jaxsnn, NEST, NEURON, and Brian2 have made\nsignificant contributions to the construction and simulation of SNNs. However, none\nof these tools currently supports automated online learning, which is the central\ntechnical contribution of BrainScale . For instance:\nNengo focuses on functional modeling using the Neural Engineering Framework and\nsupports only limited learning capabilities (e.g., PES, Voja) that are manually\ndesigned.\njaxsnn explicitly states that support for delays, or the incorporation of online\nplasticity remains future work.\nWe have added a detailed comparative discussion in Supplementary Note N, outlining\nthe limitations of these tools with respect to online learning.\nThrid , key novelty: automated online learning via compilation.\nWhile many online learning algorithms have been proposed for SNNs, no existing\nframework provides a compiler infrastructure to automate their implementation  \ncomparable to how deep learning frameworks automate training via automatic\ndifferentiation. BrainScale introduces this capability, filling a major gap in the\nSNN tooling ecosystem. It allows researchers to specify learning rules at a high\nlevel, which are then automatically compiled into efficient, event-driven,\ndifferentiable code for SNNs. This represents a conceptual and technical advance\nbeyond what current frameworks provide.\nWe have updated the manuscript to cite all relevant prior work more comprehensively,\nwhile precisely positioning the novelty of BrainScale in this broader context.\nNo. 17 / 22\n[1] Wang C, Zhang T, Chen X, et al. BrainPy, a flexible, integrative, efficient, and\nextensible framework for general-purpose brain dynamics programming[J]. elife, 2023,\n12: e86365.\nComment 3:\nThe idea of using RTRL for spiking neurons is not new either. Note that even e-\nProp (ref. 6) is a simple variant of spiking RTRL. Mingqing Xiao et al. had a\npaper in NeurIPS two years ago that is conceptually very similar to the work\npresented here. It is also cited in the abstract of the manuscript (ref. 82).\nThey used a leaky-integrate-and-fire neuron model and derived an approximate\nversion of RTRL. The approximation used a very similar diagonalization to achieve\nlinear scaling (as done in the current manuscript). They show theoretical and\nexperimental results for recurrent and feedforward networks. LIF neurons are very\nsimilar to the exponential synaptic response kernel that is used in the\nderivation of ES-D-RTRL. In the light of this similarity, the manuscript seems a\nbit unfair towards ref. 82. It is only cited in the abstract in the statement\n\"    or use oversimplified models [82], failing to capture the intricate brain\ndynamics.\" What exactly are the intricate brain dynamics that BrainScale can\ncapture and prior methods failed to address? A detailed side-by-side comparison\nwould be in need here.\nThanks for the question, which help us to clarify the novelty of our algorithm in\nlight of existing work, particularly that of Xiao et al.\nWe fully agree that using RTRL-like methods for training spiking neural networks is\nnot new; however, to solve the key scalability question, we should not focus on\nwhether methodology has not been explored before, but on the practical utility of our\nalgorithm in enabling scalable online learning for general SNNs.\n1. The contribution and novelty\nWhile the general idea of applying RTRL to spiking neurons has been explored before,\nour key contribution lies in deriving a theoretically principled, event-driven RTRL\nalgorithm (ES-D-RTRL) for general SNN models that achieves linear memory complexity .\nThis extends beyond previous approximations in both scope and rigor.\nIn contrast to OTTT (Xiao et al.), our algorithm:\n(A) Generalization to arbitrary neuron & synapse models\n– Xiao et al. (OTTT) focus on LIF neurons with single exponential synapses. In ES-D-\nRTRL, we derive an online gradient for any combination of neuron and synapse dynamics\n—LIF, ALIF, conductance-based channels, multi-exponential kernels, adaptation, reset\nand refractoriness—without altering the core algorithm.\n– This breadth is crucial for brain simulation and neuromorphic applications, where\ndiverse biophysical mechanisms interplay.\n(B) Exact O(n) memory & computation with theoretical guarantees\n– OTTT reduces complexity by dropping post-synaptic reset/adaptation terms from the\nsurrogate gradient and then empirically showing negligible impact on small networks.\n– In contrast, ES-D-RTRL performs a rank-1 decomposition of the full Jacobian,\nincluding all post-synaptic dynamics, and reformulates the double sum as a product of\ntwo vector sums. This yields true memory and per-step updates with\nrigorous bounds on gradient fidelity and convergence.\nNo. 18 / 22\n(C) Intricate brain dynamics captured by ES-D-RTRL\nPrior methods typically ignore or oversimplify:\nreset and refractoriness\nmulti-timescale adaptation\nconductance-based synaptic kinetics (e.g., NMDA vs. AMPA)\ndendritic compartment interactions\nneuromodulatory gating\nES-D-RTRL, by contrast, seamlessly integrates these elements into the gradient\nupdate, enabling large-scale, biologically realistic network training.\nTo ensure we accurately represented the distinction, we also discussed this directly\nwith Dr. Mingqing Xiao (ref. 82's first author), who confirmed the following reply\nabout the differences:\nTip\n第⼀，考虑的神经元模型不⼀样，这篇⽂章的算法针对更⼴泛的神经元和突触模型，⽬的包括了脑模拟等等，\nOTTT主要还是考虑ai的任务和简单的模型；\nFirst, the neuron models considered in the two works differ. The algorithm\nproposed in this paper is designed for a broader class of neuron and synapse\nmodels, with objectives that include brain simulation and related applications. In\ncontrast, OTTT primarily focuses on AI tasks and employs relatively simple neuron\nmodels.\n第⼆，在从O(n^2)降到O(n)复杂度那⾥的⽅法也不⼀样，OTTT是略去post-synaptic中reset或者其他可\n能的dynamics的影响（可看作对这些项不使⽤surrogate gradient），然后在相对简单的神经元设定下分析\n这个对梯度下降⽅向的影响，这篇⽂章的算法对post-synaptic的dynamics部分是有考虑的，通过秩1矩阵的\n和转成向量和的乘积进⾏简化并提供了相应的理论分析。\nSecond, the approaches used to reduce computational complexity from to\nare also distinct. OTTT achieves this by omitting the influence of post-\nsynaptic dynamics such as reset mechanisms—effectively ignoring these components\nin the surrogate gradient computation—and then analyzing the impact of this\nsimplification on the gradient descent direction under a relatively simple\nneuronal configuration. In comparison, the algorithm presented in this paper\nexplicitly accounts for post-synaptic dynamics. It simplifies the computation by\nexpressing the relevant terms as products of vector sums derived from rank-1\nmatrix decompositions, and provides corresponding theoretical analysis to support\nthis approach.\n2. Side-by-Side Comparison\nHere we provide a detailed side-by-side algorithmic comparison.\nNo. 19 / 22\nAspect OTTT (Xiao et al.) ES-D-RTRL (This Work)\nLIF, ALIF, adaptive\nNeuron models LIF only exponential, multi-\ncompartment, etc.\nExponential, Two-exponential,\nSynapse models None conductance-based, short-term\nplasticity\nTreatment of\nIgnored in surrogate Fully accounted via Jacobian\npost-synaptic\ngradient decomposition\ndynamics\nComplexity True O(n) with all dynamics\nO(n²)→O(n) by omission\nreduction retained\nFormal proofs of gradient\nTheoretical Empirical validation on\nerror bounds and convergence\nanalysis simple models/tasks\ntheorems\nLarge recurrent nets and\nLIF-based\nScalability brain-scale simulations with\nfeedforward/recurrent nets\ncomplex dynamics\n3. Exponential kernel vs. LIF dynamics\nRegarding the similarity between LIF neurons and the exponential kernel in Eq. (7):\nWhile both may involve exponential terms, their roles are fundamentally different. In\nOTTT, the exponential kernel arises directly from the LIF neuron dynamics and is\ntightly coupled to the structure of the neuron model. In contrast, our use of\nexponential kernels in Eq. (7) is part of a generic temporal filter applied to the\noutput spikes, and is independent of the internal dynamics of the neuron model (e.g.,\nLIF, ALIF, Izhikevich, Hodgkin–Huxley, etc).\nMoreover, the influence of neuron dynamics in our method is captured by Eq. (8),\nwhich models post-synaptic responses. This separation of spike filtering (Eq. 7) and\ndynamical response (Eq. 8) allows our framework to generalize well across neuron and\nsynapse types.\nComment 4:\nAlso the simulation results reported here are misleading in the sense that they\nmiss a fair comparison with state-of-the-art. 94% test accuracy on DVS gesture\ndataset is not state of the art. Modern machine learning method achieve 100% on\nthis data set, see: [4]. Why this gap? Should be discussed. Furthermore, Xiao et\nal. [ref 82] also reported superior results on DVS gesture dataset (>96%). Please\ninclude and discuss their result and compare the difference in detail, given the\nrelevance of this work (see also comment above). They also presented compelling\nresults for the more complex cifar 10 and 100 data sets, among others. Would be\ninteresting to include those as well here. Also on Spiking speech commands much\nbetter results have been reported, e.g. [5, 6]. Why this gap? Please discuss.\nFurthermore, the comparison with BPTT without surrogate is not very meaningful,\nsince BPTT for spiking neurons in that case simply doesn't work. Also, why stop\nlearning in Fig.5E,F at 95%? What is significant about this number, and how would\nconvergence behave after reaching that point?\n[4] Schöne et al. 2024A. https:  arxiv.org/pdf/2411.12603\nNo. 20 / 22\n[5] Schöne et al. 2024B. https:  arxiv.org/pdf/2404.18508\n[6] Soydan et al. 2024. https:  arxiv.org/abs/2410.03464\nThanks for the comment. We think there may be some misunderstandings on our\ncomparisons.\n1. On the comparison with state-of-the-art accuracy\nWe acknowledge that the absolute test accuracy numbers in the initial submission\n(e.g., 94% on DVS Gesture) fall short of the current state-of-the-art achieved by\nsome recent machine learning methods. However, our intention was not to claim a new\nSOTA accuracy, but to provide a fair and controlled comparison among different\nlearning algorithms    particularly BPTT, D-RTRL, and our proposed ES-D-RTRL    under\nidentical conditions. This includes using the same network architecture, same\noptimizer, and same hyperparameters, so that differences in performance are\nattributable solely to the learning algorithm itself. These results are reported in\nTable 1. We now emphasize this point more clearly in the revised text to Table 1.\nLines 265-268.\n2. Model architecture and fairness of results\nIn the initial submission, the observed performance gap was primarily due to the use\nof simple network architectures and neuron models, which were intentionally selected\nto ensure clarity and controlled analysis rather than to maximize performance. In the\nrevised version, we have incorporated more advanced neuron models (e.g., adLIF,\nRadLIF, and EGRU) and demonstrated that our learning algorithm scales effectively,\nachieving over 97% accuracy on DVS Gesture and over 93% on SHD. These improvements\nbring our results closer to those reported in prior online and offline learning\nbenchmarks, notably surpassing the results in Xiao et al. This demonstrates that our\nonline learning method supports high-performance training of advanced spiking network\narchitectures.\nThe reviewer noted that modern machine learning methods achieve 100% on this\ndataset, see: [4]. The referenced work employs deep state-space models (SSMs) for\nprocessing neuromorphic datasets. While SSMs have shown strong performance on tasks\ninvolving long-term sequential dependencies, they are rate-based rather than spiking\nmodels, and thus fall outside the scope of this manuscript. Moreover, such models are\nincompatible with the ES-D-RTRL framework we explore. For these reasons, we did not\nevaluate SSMs within our online learning system.\nFinally, our current evaluations  covering a variety of SNN models and benchmarks  \nalready demonstrate the strong performance of our learning algorithms. Therefore, we\nchose not to include additional datasets such as CIFAR-10 or CIFAR-100, which are\nbeyond the primary focus of this study.\n3. Stopping at 95% in Fig. 5E,F\nThe learning was stopped at 95% in Fig. 5E,F (the previous manuscript) or at 90% in\nFig. 5C (the revised manuscript) since the real animals that participated in this\ntask do not exhibit 100% accuracy [1]. Models show 90% accuracy is enough to study\nthe neural mechanism that how the mice completes this task.\nOn the other hand, we also evaluated the model's learning performance over very long\ntime scales. Supplementary Fig. S21 shows that the training accuracy steadily\nimproves as the number of learning epochs increases.\nNo. 21 / 22\n[1] Morcos A S, Harvey C D. History-dependent variability in population dynamics\nduring evidence accumulation in cortex[J]. Nature neuroscience, 2016, 19(12): 1672-\n1681.\nComment 5:\nCode will be published after acceptance. The provided links to source code are\nnot active yet.\nIn this version, we open source all our code. BrainScale now can be freely available\non the GitHub.\nNo. 22 / 22"}, {"round": 2, "reviewer_comments": {"reviewer_1": "The authors have made large strides in addressing the previous round of reviews. One of the main issues was to relate their\nwork to the recent literature on training SNNs on datasets such as SHD. The results of the proposed methods are now\nreported to be much more competitive, presumably by using more advanced neuron models.\nHowever, the authors have chosen to rerun methods of competing works while removing some aspects of their methods -\nmost notably batch norm. They argue that this is a fairer comparison of learning methods as batch norm is incompatible with\nonline learning, which is their own aim.\nThis is a problematic argument to make. Fair comparisons are notoriously difficult, but in this case, the offline, back-\npropagation-through-time methods yield trained SNNs with the stated performance and any normalisation with batch norm\nduring the training phase can be absorbed into weights during inference. So, in essence, there is no fundamental argument\nagainst the use of batch norm. Removing it to make the comparison with a method that can't employ batch norm more\nfavourable to that other method doesn't seem to improve fairness.\nAdditionally, the authors are still missing at least one other training method which also achieves competitive performance\n(without batch norm), such as the results based on EventProp (e.g. 10.1088/2634-4386/ada852). Eventprop is also missing\nin table S1, though it is very competitive, and in their otherwise nice simulator comparison (table S2), (ml)GeNN is strangely\nabsent.", "reviewer_2": "I would like to thank the authors for their hard work in preparing the new version of the article and for taking my comments\ninto account. I appreciate the corrections that have been made. Thank you for enriching the Supplementary Notes with the\nsections \"N Comparison with other frameworks\" and \"M Comparison with other learning algorithms\", which allow the\nauthors' achievements to be assessed against other existing solutions.\nThe authors were unable to carry out some of the comparative analyses suggested in the previous review. This is mainly\ndue to limited access to specialised computing equipment. I do not consider this to be a shortcoming, as devices dedicated\nto SNN are niche solutions and access to them is severely limited in terms of time. I would like to thank the authors for their\nefforts in preparing the comparative analyses included in the article. They are entirely satisfactory.\nThe authors did not follow the suggestion to compare the results obtained with their solutions with the results obtained using\nthe conversion method (comment no. 7). The authors gave two reasons for not doing so. I do not entirely agree with the first\nargument, which suggests that such analyses are very costly. A counterargument could be that conversion algorithms are\nimplemented in commonly used frameworks, where examples using the N-MNIST dataset, which the authors also use, are\noften implemented. Therefore, conducting appropriate comparative tests should not be a very time-consuming task.\nHowever, I agree with the second argument concerning the limitation of the research to online learning algorithms. I accept\nthe authors' decision and have no comments in this regard. I recommend, however, that future research consider the steps\nsuggested by the authors in the Limitations and Future Work section in lines 413-419.\nFinally, I have two editorial comments:\n1. Re: comment 2 point 3:\nI have doubts about the values V_reset, V_rest, V_th presented in tables S3, S4, S8, S9, S10, S11, S12. Values such as “1.”\nor “0.” are given. This may be a source of confusion. I think it would be better to write “1.0 V” and “0.0 V”.\n2. The meaning of TPU is not explained in the text. For people outside the field of SNN, it may not be as obvious as, for\nexample, CPU or GPU. I think it would be a good idea to consistently explain all three abbreviations in the text of the article.\nIn my opinion, after taking into account the above two minor comments, the work can be published.", "reviewer_3": "I want to thank the authors for the clarifications. The presentation has improved, but my main concerns remain.\n- For the claimed main innovation, the online learning abilities, the approximations use in this work are already well known.\nEsp. the approach taken earlier in Marschall et al. 2020 is quite similar and seems very relevant to contrast with the\nproposed method, which is now presented as if it were original work.\n- I still disagree with the argumentation of the authors regarding the choice of benchmarks and baseline algorithms. If the\nfocus of this work is to demonstrate scaling, which the name of the model, the abstract, introduction and discussion\nsuggest, then scaling to state-of-the-art learning performance is an integral part of that. Being able to train large(r) models does not\nprove that the algorithm scales, if the performance is not matched on hard tasks. Hence, if scalability of existing algorithms\nwas the problem, then it should be solved by the proposed method, which has to be demonstrated to support the many\nclaims regarding the superiority of the method.\n- To match common scientific standards, the comparison between algorithms should involve individual hyperparameter\ntuning (not as the authors suggest, use the same hyperparameters), since choice of parameters may influence algorithmic\nperformance. I.e. if hyperparameter tuning was done for one model, the same hyperparameter tuning has to be done for all\nothers. Tuning parameters for one algorithm and then applying the same hyperparameters for comparisons is considered\nbad practice.\n- As a side remark: I agree with the other reviewers that the name clash with BrainScaleS is very unfortunate. I also don't\nthink BrainScale is a particularly good and descriptive name for the proposed model.\nMarschall et al. 2020: https://jmlr.org/papers/volume21/19-562/19-562.pdf"}, "author_response": "Dear reviewers,\nThank you for considering our manuscript, titled as A general, efficient, and\nautomated online learning system for spiking network modeling, for publication in\nNature Communications . We are grateful for the thoughtful and constructive feedback\nprovided from all reviewers; your comments substantially improved the quality and\nclarity of our work.\nOver the past three months, we carefully addressed every point through substantial\nrevisions and comprehensive experiments. The attached, point-by-point document\ndetails our responses to each comment. For ease of review, all textual changes are\nhighlighted with tracked changes in the revised manuscript.\nBefore addressing each comment in detail, we would like to summarize the major\nchanges in the current revision:\n1. Learning performance and hyperparameter tuning: In response to concerns about\nlearning performance and unfair hyperparameter setting, we have conducted\nextensive parameter searches on both BPTT and pp-prop methods using the SHD\ndataset. We achieved a testing performance of 95.33 0.11% on SHD, demonstrating\nthe robustness of our approach.\n2. Novelty compared to Marschall et al. (2020) : We have provided a detailed\ncomparison demonstrating how our methods differ from Marschall et al. (2020) ,\nclarifying the distinct contributions of our work.\n3. Clarification of \"scaling\" and \"scalable\" claims: In response to concerns about\npotentially ambitious interpretations of these terms, we have carefully revised\nthe manuscript (including title, introduction, and discussion) to clarify the\nthree core innovations of our framework: a general , efficient , and automated\nonline learning system for spiking network modeling, rather than summarize them\nas a misleading term \"scalable\".\n4. Package naming: To avoid potential confusion between the names \"BrainScaleS \"\nand \"BrainScale \", we have renamed our package BrainTrace, highlighting its\nfocus on eligibility trace–based online learning algorithms for brain modeling.\n5. Algorithm naming: Furthermore, we designate the algorithm as \"pp-prop \" (pre–\npost eligibility trace propagation) algorithm to explicitly convey its key\nmechanism based on pre- and post-synaptic eligibility trace tracking and\npropagation.\nNo. 1 / 24\nResponse to Reviewer #1\nComment 1\nThe authors have chosen to rerun methods of competing works while removing some\naspects of their methods - most notably batch norm. They argue that this is a\nfairer comparison of learning methods as batch norm is incompatible with online\nlearning, which is their own aim.\nThis is a problematic argument to make. Fair comparisons are notoriously\ndifficult, but in this case, the offline, back-propagation-through-time methods\nyield trained SNNs with the stated performance and any normalisation with batch\nnorm during the training phase can be absorbed into weights during inference. So,\nin essence, there is no fundamental argument against the use of batch norm.\nRemoving it to make the comparison with a method that can't employ batch norm\nmore favourable to that other method doesn't seem to improve fairness.\nWe thank you for raising this important point. We fully understand the concern\nregarding the use of batch normalization (BN) and have carefully reconsidered our\nexperimental design choices in the comparative study.\nFirst, we would like to further clarify why BN is fundamentally incompatible with a\ntrue online learning paradigm during training. BN is inherently a mini-batch–based\nnormalization operation that depends on the global statistics (mean and variance)\ncomputed over a batch of samples and across all time steps. This violates the core\nassumption of online learning, which processes data sequentially and cannot access\nfuture samples or global statistics. Therefore, BN cannot be applied in a strictly\nonline setting, and we intentionally excluded this component from our online learning\nframework.\nTo comprehensively address this concern, in this revision we performed systematic\nablation experiments comparing the effects of different normalization methods on both\nBPTT and pp-prop (see SI. P for details).\nFinding 1: The effect of BN depends on the complexity of neuronal\ndynamics\nAs shown in Figure S6, in BPTT training, batch normalization improved performance\nonly for the simpler RLIF network. In contrast, for more complex neuron models such\nas RadLIF, BN actually degraded performance. This observation provides an important\ninsight: when neuron dynamics are sufficiently rich, the network can achieve strong\nperformance without relying on global normalization operations such as BN. This also\nexplains why our online method performs well without BN when using complex neuron\nmodels.\nNo. 2 / 24\nFinding 2: pp-prop achieves competitive performance without BN\nIn RLIF networks, we observed that with a properly chosen surrogate-gradient\nfunction, the pp-prop online learner without BN can reach performance comparable to\nBPTT with BN (90.1% vs 90.7%; see Supplementary Figure S11 vs Figure S6).\nMoreover, as shown in Figure S6, even within the BPTT framework, the combination of\nBN + RLIF yields lower accuracy than the RadLIF network without BN. This suggests\nthat increasing neuronal complexity may be a more effective strategy than introducing\nBN. Our experiments further demonstrate that pp-prop benefits more from this strategy\nand achieves better training performance under such conditions.\nNo. 3 / 24\nFinding 3: Complex neuronal dynamics confer unique advantages for\nonline learning\nMore importantly, pp-prop shows a natural affinity for complex neuronal dynamics (see\nFigure 4 for approximation accuracy evaluation in the main text). On the SHD dataset,\npp-prop combined with complex neurons readily achieves > 94% accuracy without BN or\ndata augmentation (Figures S9 and S11). In contrast, BPTT typically requires\nadditional data augmentation to consistently surpass the same threshold (Figure S8).\nConclusion\nOverall, our intention was not to “beautify” online learning by removing BN, but\nrather to present fair and transparent comparisons between two realistic usage\nscenarios. The new experiments on the SHD dataset clearly indicate that BN is not a\nprerequisite for high performance. With appropriately chosen parameters and complex\nneuronal dynamics, pp-prop can achieve performance comparable to BPTT + BN while\nfully adhering to the constraints of true online learning.\nComment 2\nAdditionally, the authors are still missing at least one other training method\nwhich also achieves competitive performance (without batch norm), such as the\nresults based on EventProp (e.g. 10.1088/2634-4386/ada852). Eventprop is also\nmissing in table S1, though it is very competitive, and in their otherwise nice\nsimulator comparison (table S2), (ml)GeNN is strangely absent.\nWe thank you for the insightful comments. We highly acknowledge the\nrepresentativeness of EventProp and (ml)GeNN in the field of spiking neural network\n(SNN) training and simulation. In this revision, we have added a detailed discussion\nof EventProp in SI. M and of (ml)GeNN in SI. N. Furthermore, we have included a\nNo. 4 / 24\ncomparison of EventProp’s best training performance [4] on the SHD dataset using the\nRLIF model (see Table 1 in the main text).\n1. EventProp follows a different training paradigm from this work\nEventProp represents an event-driven, offline backpropagation algorithm. Although\nits memory usage scales with the number of spikes rather than timesteps, it still\nrequires full backward propagation at each time step. In contrast, our work focuses\non online learning algorithms , which perform real-time, forward-only updates with\nO(N) memory complexity, without backward time unfolding or full activity storage.\nIn practice, EventProp achieves roughly a fourfold improvement in memory efficiency\ncompared with BPTT on the SHD dataset [4]. However, our pp-prop algorithm typically\nprovides tens to hundreds of times higher memory efficiency (see Figure 4 memory and\nspeed evaluation in the main text).\nTherefore, to avoid conceptual confusion, Table 1 primarily compares algorithms that\nsatisfy the online constraint, such as e-prop, OSTL, OTTT, and our pp-prop.\nNevertheless, following the reviewer’s suggestion, we have also added several\nonline/offline spiking models that achieve state-of-the-art (SOTA) performance for\nbroader context.\nFurthermore, the competitive results of EventProp on the SHD dataset are largely\nattributed to extensive data augmentation and machine learning tricks, rather than\nintrinsic superiority of the training algorithm itself. As reported in the original\npaper, before these enhancements, the accuracies were modest -- ( for\nfeedforward and for recurrent models, Fig. 4 in [4]) -- slightly below\nprior BPTT-based SNN results. After employing data augmentation and architectural\noptimization (Section 2.3 of the original paper [4]), performance on the SHD test set\nimproved to .\nTo ensure fair comparison, we also tested our online learning algorithms under\nsimilar data augmentation settings. We found that RLIF networks trained with pp-prop\nachieved around 94% accuracy, demonstrating that our online learning approach\nachieves even stronger performance under comparable experimental conditions.\nNo. 5 / 24\n2. (ml)GeNN does not implement a general-purpose online learning\ncompiler\nWe agree that our earlier comparison of simulation frameworks omitted (ml)GeNN.\n(ml)GeNN is a high-performance, GPU-optimized SNN machine learning library built upon\nGeNN. It provides a high-level Keras-like interface that simplifies model prototyping\nand leverages GeNN’s sparse data structures for efficient simulation. The framework\nsupports several representative training algorithms (e.g., e-prop and EventProp) and\noffers excellent engineering value for rapid prototyping.\nHowever, (ml)GeNN’s design does not address the “general and automated online\nlearning compilation ” that our work targets. From both its source code and\ndocumentation, it is evident that mlGeNN’s training compilers currently support only\npredefined neuron dynamics and loss functions. For example, the EventProp compiler\nsupports a fixed LIF hidden-layer model and a limited set of task types\n(classification via SparseCategoricalCrossentropy, regression via MeanSquareError)\n[2]; similarly, the e-prop compiler is restricted to LIF and ALIF architectures [3].\nThus, mlGeNN does not yet provide a universal compiler that can generate online\nlearning updates for arbitrary user-defined neuron, synapse, or learning-rule\ndefinitions.\nIn contrast, our system is designed around a unified model abstraction and compiler\narchitecture, which can automatically generate online learning code for any spiking\nnetwork model under the same abstraction. This is the fundamental conceptual\ndistinction between our framework and mlGeNN. We have clarified mlGeNN’s functional\nscope and design philosophy in the revised table notes (Table S2) and expanded the\nrelated discussion and citations in the supplementary materials (Lines 1762-1764,\nLines 1811-1823) to ensure a more complete and fair comparison.\nReferences\n[1] Wunderlich, T.C., Pehle, C. Event-based backpropagation can compute exact\ngradients for spiking neural networks. Sci Rep 11, 12829 (2021). https:-/doi.org/10.\n1038/s41598-021-91786-z\n[2] https:-/github.com/genn-team/ml_genn/blob/master/ml_genn/ml_genn/compilers/event_\nprop_compiler.py\n[3] https:-/github.com/genn-team/ml_genn/blob/68dfb63ef6fb824e06554f4050b89acab422f21\nc/ml_genn/ml_genn/compilers/eprop_compiler.py#L127-L221\n[4] Nowotny T, Turner J P, Knight J C. Loss shaping enhances exact gradient learning\nwith eventprop in spiking neural networks [J]. Neuromorphic Computing and\nEngineering, 2025, 5(1): 014001.\nNo. 6 / 24\nResponse to Reviewer #2\nThank you for your feedback. We are grateful for your understanding regarding the\nlimitations we faced with specialized computing equipment and comparative analyses.\nComment 3\nRegarding other two editorial comments, we have addressed them immediately:\nFinally, I have two editorial comments:\n1. Re: comment 2 point 3:\nI have doubts about the values V_reset, V_rest, V_th presented in tables\nS3, S4, S8, S9, S10, S11, S12. Values such as \"1.\" or \"0.\" are given. This\nmay be a source of confusion. I think it would be better to write \"1.0 V\"\nand \"0.0 V\".\n2. The meaning of TPU is not explained in the text. For people outside the\nfield of SNN, it may not be as obvious as, for example, CPU or GPU. I think\nit would be a good idea to consistently explain all three abbreviations in\nthe text of the article.\nIn my opinion, after taking into account the above two minor comments, the work\ncan be published.\n1. Value notation in tables: We have revised all voltage values in Tables S3, S4, S8,\nS9, S10, S11, and S12 to use the clearer format \"1.0 mV\" and \"0.0 mV\" as you\nsuggested. Actually, these values are normalized to the range [0, 1] from their\noriginal biophysical magnitudes (e.g., a resting potential of –70 mV and a spiking\nthreshold of –50 mV) to ensure numerical stability during training. Moreover, we have\nalso added corresponding physical units to other symbols in all tables in SI, such as\n, , , .\n2. Explanation of abbreviations: We agree that TPU (Tensor Processing Unit) should be\nexplained for readers outside the immediate field. We have added clear definitions\nfor all three abbreviations (CPU - Central Processing Unit, GPU - Graphics Processing\nUnit, and TPU - Tensor Processing Unit) at their first occurrence in the main text to\nensure consistency and accessibility for all readers (see Lines 53-54).\nNo. 7 / 24\nResponse to Reviewer #3\nComment 4\nFor the claimed main innovation, the online learning abilities, the\napproximations use in this work are already well known. Esp. the approach taken\nearlier in Marschall et al. 2020 is quite similar and seems very relevant to\ncontrast with the proposed method, which is now presented as if it were original\nwork.\nMarschall et al. 2020: https:-/jmlr.org/papers/volume21/19-562/19-562.pdf\nWe thank you for pointing out the connection between this study and Marschall et\nal. (2020) [1]. We believe that Marschall et al. (2020) actually demonstrates the\ndistinct originality of our algorithm. While both studies belong to the broad\ncategory of \"achieving efficient online learning through RTRL approximation,\" they\ndiffer fundamentally in three dimensions: research targets, approximation principle,\nand final results. These differences clearly demonstrate that this study is not a\nrepetition of existing approximation methods or a simple application to SNNs, but\nrather represents a substantial breakthrough specifically addressing the requirements\nof spiking neural networks and large-scale brain simulation.\n1. Different research targets: from continuous-valued RNNs to event-\ndriven SNNs\nMarschall et al. (2020) focuses on conventional, non-spiking recurrent neural\nnetworks. The unified framework they proposed and the algorithms analyzed therein\n(e.g., UORO [2], KF-RTRL [3] —— all reviewed and compared in our SI. M, Lines 1692-\n1697) are specifically designed for RNNs with continuous-valued hidden states (see\nSections 3–4 of the original paper [1]).\nThe essence of these methods lies in reducing the memory complexity of full\nRTRL to through tensor or low-rank decompositions — for example, the\nfactorization in UORO [2], or Kronecker decompositions in KF-RTRL [3]\nor OK [4]. However, their derivation remains confined to non-spiking networks that\ndo not incorporate biophysical neuronal dynamics (see Table 1 in Marschall et al. ,\n2020).\nIn contrast, our study explicitly targets spiking neural networks , which are\ncharacterized by event-driven computation and intrinsic temporal dynamics. These core\nfeatures fall entirely outside the assumptions and scope of the RNN framework\nanalyzed by Marschall et al. Leveraging these unique biophysical properties of SNNs,\nour proposed pp-prop algorithm achieves innovations that go beyond the Marschall\nframework (as detailed below).\nNo. 8 / 24\n2. Different approximation principles: from low-rank assumptions to\nSNN-specific mechanisms\nThe various methods summarized by Marschall et al. (2020) (such as UORO [2], KF-RTRL\n[3], OK [4] and RFLO) all rely on the premise that \"eligibility traces possess a low-\nrank structure.\" Generally speaking, these methods employ the following general\napproximation approach :\nNote\nA. Basic form of RTRL\nIn RTRL, we define\nThe core recursive formulation of RTRL is:\nThis update requires maintaining and updating a matrix of size ( ),\nresulting in a computational complexity of and memory complexity of\n, making it extremely difficult to scale.\nB. General approximation framework for RTRL\nTo reduce RTRL’s computational and memory costs, Marschall et al. (2020) proposed\na unified approximation framework adopted by most subsequent methods. The general\nsteps are as follows:\nChoose a low-rank representation (e.g., Kronecker product) to approximate\n;\nExpress both and in the same representation format;\nCombine the two terms while maintaining unbiasedness;\nMinimize variance as much as possible.\nThe core properties of this framework can be summarized as:\nUnbiasedness:\nLow-rank structure:\nMinimum variance: Among all unbiased approximations, choose the one with the\nsmallest variance.\nC. Comparison of UORO, KF-RTRL, and OK approximation methods\nThe online learning algorithms summarized by Marschall et al. (2020) generally\nconform to this structure.\n(a) UORO — Unbiased Online Recurrent Optimization\nUORO [2] adopts a rank-1 Kronecker product approximation:\nNo. 9 / 24\nIt computes the first term as:\nand approximates the second term as:\nThe two components are then merged using the sign trick :\nThe advantage of UORO is simple computation and low memory usage, but its\ndisadvantage is high variance and slow convergence.\n(b) KF-RTRL — Kronecker-Factored RTRL\nKF-RTRL [3] also uses a rank-1 Kronecker product form:\nLeveraging structured networks such as LSTM or RHN, it represents exactly as:\nThen it computes:\nand again combines terms using the sign trick :\nKF-RTRL has lower variance than UORO, but it remains a 1-Kronecker approximation\nwith room for variance improvement.\n(c) OK — Optimal Kronecker-Sum Approximation\nOK [4] has the following approximation form:\nIts core idea is to represent as an r-Kronecker sum rather than a single\nKronecker product. Therefore, during mixing, the problem is transformed into:\napproximating an Kronecker sum with Kronecker products. OK uses Minimum\nVariance Unbiased Estimator (MVUE) theory to construct the optimal low-rank\napproximation.\nFirst, given:\nSecond, construct an orthogonal basis:\nNo. 10 / 24\nNext, construct matrix .\nThen, perform SVD on to obtain its optimal rank-r approximation .\nFinally:\nis the minimum variance unbiased approximation.\nThe advantage of OK is minimum variance, theoretically optimal, and can be\nextended to arbitrary , balancing accuracy and efficiency. However, its\ncomputational complexity is O( ), higher than KF-RTRL.\nNow, let us revisit the core principles of the pp-prop algorithm presented in the\nmain text section BrainTrace achieves linear-memory online learning algorithm (Lines\n85-120) and Appendix SI. G Derivation of online learning algorithms. pp-prop does not\nimpose prior assumptions on the eligibility trace; rather, it performs stepwise\nreduction based on the network's inherent properties. Specifically, the key\ninnovation of pp-prop lies in:\nNote\nA. Leveraging spike sparsity and input sign consistency\npp-prop exploits a unique property of SNNs: the input signals (e.g., spikes)\nhave consistent signs (nonnegative). This allows the sum of multiple rank-1 outer\nproducts\nto be approximated as\nThis yields a rank-1 approximation that is deterministic and does not rely on\nrandom projections or the sign trick used by methods such as UORO or KF-RTRL.\nB. Decomposition of the eligibility trace into pre- and post-synaptic components\nBuilding on this property, pp-prop introduces a novel factorization of the\neligibility trace:\nNo. 11 / 24\nwhere tracks presynaptic activity (input) and tracks\npostsynaptic Jacobian dynamics.\nBoth components are updated through exponential smoothing, requiring only\nmemory—effectively linear in network size.\nC. Compatibility with complex neuron and synapse models (e.g., ALIF, GIF, STP,\nSTD)\nUORO, KF-RTRL, and OK were designed primarily for standard RNNs or LSTMs .\nHowever, pp-prop is designed to be applicable to: multi-variable neuron models\n(e.g., membrane potential, adaptation current, synaptic conductance), synaptic\nplasticity mechanisms (short-term facilitation/depression), excitatory–inhibitory\n(E/I) balanced networks, and sparse or structured connectivity.\nThe systematic comparison presented above reveals that pp-prop is fundamentally\ndistinct from the general approximation methods discussed previously in three key\naspects: approximation strategy, structural assumptions, and applicability scope.\nNote\nA. Distinction 1: fundamental shift in approximation target\nMethods such as UORO/KF-RTRL/OK directly approximate the gradient matrix\nitself, consistently focusing on how to preserve the matrix addition\nin an unbiased manner within a low-rank subspace.\nIn contrast, pp-prop does not approximate . Instead, it redefines the\ncomputation of the eligibility trace, approximating the temporal integration of\nthe trace itself as separable into two independent presynaptic and postsynaptic\ncomponents.\nThe fundamental difference lies in: the former employs \"matrix factorization\",\nwhile the latter adopts \"temporal compression\".\nB. Distinction 2: specialized vs. general structural assumptions\nUORO/KF-RTRL/OK design approximation frameworks for general RNNs, utilizing only\nlinear algebraic structures (Kronecker products) without dependence on specific\nneuron types.\nIn contrast, pp-prop is specifically designed for SNNs and deeply leverages three\nbiological properties inherent to spiking networks:\nSparse activation → diagonal-dominant Jacobian\nPhysical connectivity → Kronecker decomposability\nSign preservation (non-negative spikes) → temporal low-rank theorem (Theorem\nC.3, our original contribution)\nTherefore, the fundamental distinction is that pp-prop's approximation is non-\ntransferable to non-spiking RNNs.\nC. Distinction 3: breakthrough in memory complexity\nNo. 12 / 24\nThe memory complexity of UORO/KF-RTRL/OK remains super-linear ( ), as\nthey require storing the Kronecker factors of .\npp-prop achieves, for the first time, linear memory complexity by storing\nonly two vectors (pre/post traces), independent of the number of connections .\nD. Distinction 4: trade-off strategies for unbiasedness\nUORO/KF-RTRL/OK pursue statistical unbiasedness and are theoretically strictly\nunbiased, with as a core design constraint (achieved through sign\ntricks or MVUE).\npp-prop pursues a balance between biological plausibility, computational\nfeasibility, and accuracy. It permits biased approximation but guarantees\nasymptotic exactness through theoretical proof (Theorem C.3) and empirical\nvalidation (cosine similarity > 0.99).\nConclusion\nIn summary, the fundamental innovation of pp-prop lies in its departure from the\n\"matrix approximation\" paradigm of the Marschall et al. framework, instead\nexploiting the physical constraints of SNNs (sparsity, non-negativity, local\nconnectivity) to reformulate the mathematical structure of the problem. Rather\nthan pursuing gradient precision and unbiasedness approximations to RTRL as in\nUORO/KF-RTRL/OK, pp-prop is a custom-tailored RTRL approximation for SNNs,\nemphasizing: linear memory, biological plausibility, long-sequence online\nlearning, and system automation.\nIn short: pp-prop is not a “more accurate” RTRL approximation — it is a “more\nsuitable” gradient approximation for SNNs.\n3. Different computational complexity and engineering practicality\nThe RNN online learning approximation methods summarized by Marschall all remain at\nmemory, whereas our pp-prop algorithm exploits SNN properties to reduce\ncomplexity to , requiring only two vector-level traces (presynaptic and\npostsynaptic traces). This fully demonstrates the unique advantages of SNNs and\nrepresents a more biologically reasonable implementation.\nFurthermore, we applied the pp-prop algorithm to a Drosophila whole-brain SNN\nsimulation task (125,000 neurons, 50 million synapses) and found that it requires\nonly 8.9GB of GPU memory, whereas approximation methods and BPTT require over\n100GB of memory. This makes whole-brain SNN training possible for the first time.\nFinally, Marschall's work focuses on theoretical frameworks and algorithmic\ncomparisons, without providing scalable implementation tools. Algorithms such as\nUORO, KF-RTRL, and OK often suffer from low computational efficiency due to complex\ntechniques such as low-rank approximation and SVD decomposition, making large-scale\nengineering implementation and application difficult. In contrast, the pp-prop\nalgorithm proposed in this study is simpler and facilitates the development of a\nNo. 13 / 24\ncomplete online learning compiler, enabling automation from model definition to\nhardware deployment. This fills the gap of lacking system-level support in the work\nsummarized by Marschall.\nIn short: pp-prop is not a proof-of-concept experimental algorithm, but an\nadvance that can be truly implemented and widely applied.\n4. Citation of prior work\nIn the revised manuscript, we have incorporated explicit citations to Marschall et\nal. (2020) [1] at key locations in the text.\nSpecifically, at Lines 47-48, we cite Marschall et al. (2020) in the context of RTRL\nreformulation to acknowledge the generality of existing approximation approaches:\nFor efficiency, BrainTrace, leverages the sparse, event-driven nature of spiking\nneurons to reformulate real-time recurrent learning (RTRL) [13, 14], ...\nFurthermore, at Lines 321-322, we reference Marschall et al. (2020) when\nhighlighting our algorithm's advantages, noting that conventional RNN approximations\nincur at least quadratic memory complexity:\nApproximating online gradient computations in recurrent neural networks has been\nextensively studied, yet existing effective approximations exhibit at least\nquadratic memory complexity [14].\nReferences\n[1] Owen Marschall, Kyunghyun Cho, and Cristina Savin. A unified framework of online\nlearning algorithms for training recurrent neural networks . JMLR, 21(135):1–34, 2020.\n[2] Tallec C, Ollivier Y. Unbiased online recurrent optimization[J]. arXiv preprint\narXiv:1702.05043, 2017.\n[3] Mujika, Asier, Florian Meier, and Angelika Steger. \"Approximating real-time\nrecurrent learning with random kronecker factors.\" Advances in neural information\nprocessing systems 31 (2018).\n[4] Benzing F, Gauy M M, Mujika A, et al. Optimal kronecker-sum approximation of real\ntime recurrent learning[C]-/International Conference on Machine Learning. PMLR, 2019:\n604-613.\nComment 5\nI still disagree with the argumentation of the authors regarding the choice of\nbenchmarks and baseline algorithms. If the focus of this work is to demonstrate\nscaling, which the name of the model, the abstract, introduction and discussion\nsuggest, then scaling to state-of-the-art learning performance is an integral\npart of that. Being able to train large(r) models does not prove that the\nalgorithm scales, if the performance is not matched on hard tasks. Hence, if\nscalability of existing algorithms was the problem, then it should be solved by\nNo. 14 / 24\nthe proposed method, which has to be demonstrated to support the many\nclaims regarding the superiority of the method.\nWe thank you for raising this important point.\nWe agree that one ambitious interpretation of \"scaling\" is scaling to state-of-the-\nart performance on the hardest benchmarks . That is a valuable goal for the field. Our\nwork, however, targets a different/complementary dimension of scaling: making online\ntraining general, memory- and compute-efficient, and automatically applicable to\narbitrary user-defined spiking models. To eliminate any room for over-interpretation,\nwe have revised the paper to make this scope explicit and to temper language that\ncould be read as claims of benchmark superiority.\n1、What we mean by \"scaling\" (clarified in the revised manuscript).\nIn our manuscript we define scaling along three axes:\nGenerality. The framework applies to arbitrary SNN architectures and learning\nfunctions (user-defined neurons/synapses, heterogeneous dynamics, custom\nlocal/global objectives), rather than a fixed catalog.\nEfficiency. The online algorithm maintains linear memory in the number of states\nand constant memory in sequence length, avoiding the memory growth of\nbackpropagation through time or of SOTA online learning methods. This\nenables training at long horizons and with larger models where prior\nonline/offline methods become infeasible.\nAutomation. An online-learning compiler eliminates manual\nderivations/implementations of eligibility traces for new models, lowering the\nbarrier to adoption and reducing error-prone re-engineering.\nWe have adjusted the manuscript to consistently use the phrasing \"a general,\nefficient, and automated online learning system for spiking network modeling\"\nrather than the umbrella term \"scalability\" used before.\n2、High-performing and linear-memory system enables whole-brain\ntraining of Drosophila models.\nThe most compelling demonstration of the scalability of our system is its ability to\nperform end-to-end training of a Drosophila whole-brain spiking neural network\nmodel.\nWith the rapid progress in large-scale synaptic connectome and neural activity\nimaging, integrating multi-scale and multi-modal neuroscience data to construct\nfunctional whole-brain models has become a central challenge. Taking Drosophila as\nan example, its complete connectome has been reconstructed at single-synapse\nresolution [1], while whole-brain calcium imaging experiments have captured neural\nactivity across nearly the entire brain [2, 3]. Based on these advances, a natural\nand crucial objective is to train a connectome-constrained whole-brain SNN capable of\nreproducing experimentally recorded brain-wide activity. However, this task is\ncomputationally formidable:\nNo. 15 / 24\n(1) Whole-brain calcium imaging in Drosophila typically has a sampling rate of ~1\nHz, implying that fitting a single data point requires simulating approximately one\nsecond of network activity, i.e., about 1,000 time steps (with a 1 ms integration\nstep). Consequently, training such a model with backpropagation through time (BPTT)\nwould demand roughly 1.1 TB of GPU memory (Fig. S25).\n(2) The Drosophila connectome (FlyWire v630) contains over 12,500 neurons (more than\n14,000 in the latest v783 release). Using an online learning algorithm with\nmemory complexity —— such as D-RTRL —— would still require over 120 GB of GPU memory\n(Fig. S25).\nIn contrast, our proposed pp-prop algorithm achieves only memory complexity in\ntheory, enabling efficient training of the entire Drosophila whole-brain SNN model\non a single GPU. This achievement makes connectome-constrained whole-brain training\nfeasible for the first time under limited computational resources, fully\ndemonstrating the core value of our algorithm: making whole-brain simulation and\ntraining truly practical.\nFurthermore, our BrainTrace framework automatically compiles the pp-prop algorithm,\nallowing users to focus solely on defining the whole-brain SNN model without\nconcerning themselves with implementation details. This automation mechanism\nsignificantly enhances the convenience and generalizability of model construction,\nestablishing a scalable foundation for whole-brain training of more complex network\narchitectures and neural dynamics models.\n3、SOTA performance on neuromorphic benchmarks\nIn our latest experiments, by leveraging well-established neuronal dynamics models\n[4] and machine learning techniques [5] from the field, we have achieved state-of-\nthe-art results on common neuromorphic tasks. Notably, on the SHD dataset, our online\nlearning approach achieved an average accuracy of 95.33%. This result substantially\nexceeds existing online learning methods and is even comparable to offline learning\nmethods that employ extensive network architectures and training techniques (Table 1,\nor see the screenshot below).\nNo. 16 / 24\nReferences\n[1] Dorkenwald, Sven, et al. \"Neuronal wiring diagram of an adult brain.\" Nature\n634.8032 (2024): 124-138.\n[2] Mann, Kevin, Courtney L. Gallen, and Thomas R. Clandinin. \"Whole-brain calcium\nimaging reveals an intrinsic functional network in Drosophila.\" Current Biology\n27.15 (2017): 2389-2396.\n[3] Turner, Maxwell H., Kevin Mann, and Thomas R. Clandinin. \"The connectome predicts\nresting-state functional connectivity across the Drosophila brain.\" Current biology\n31.11 (2021): 2386-2394.\n[4] Bittar, Alexandre, and Philip N. Garner. \"A surrogate gradient spiking baseline\nfor speech command recognition.\" Frontiers in Neuroscience 16 (2022): 865897.\n[5] Nowotny, Thomas, James P. Turner, and James C. Knight. \"Loss shaping enhances\nexact gradient learning with eventprop in spiking neural networks.\" Neuromorphic\nComputing and Engineering 5.1 (2025): 014001.\nComment 6\nTo match common scientific standards, the comparison between algorithms should\ninvolve individual hyperparameter tuning (not as the authors suggest, use the\nsame hyperparameters), since choice of parameters may influence algorithmic\nperformance. I.e. if hyperparameter tuning was done for one model, the same\nhyperparameter tuning has to be done for all others. Tuning parameters for one\nalgorithm and then applying the same hyperparameters for comparisons is\nconsidered bad practice.\nWe sincerely thank you for raising the important point regarding hyperparameter\noptimization. We fully agree that, for a fair and rigorous comparison of algorithmic\nperformance, each method should undergo independent hyperparameter tuning to avoid\npotential biases introduced by parameter settings that may favor a particular\napproach.\nIn the previous version of the manuscript, we did not deliberately fine-tune our\nproposed method and then transfer its optimized hyperparameters to competing\nalgorithms such as BPTT. The main reason is that BPTT, D-RTRL, and pp-prop are\nfundamentally different implementations of the same underlying gradient computation\nprinciple, all applicable to identical network architectures and optimizers.\nTherefore, we aimed to maintain consistent training configurations across methods to\nensure comparability.\nIn this paper revision, we have conducted systematic and comprehensive hyperparameter\noptimization for both the BPTT and pp-prop algorithms on the SHD dataset.\nThe objective of this tuning is to identify the optimal parameter combination for\neach algorithm to achieve the best learning performance under specific network\narchitectures. The optimization encompasses\nlearning rate,\nsurrogate gradient function types,\nnormalization methods,\nand data augmentation strategies,\nNo. 17 / 24\namong other aspects.\nAll experimental details and results are systematically reported in Supplementary\nMaterials SI. P, and please run the attached code for reproduction.\nStep 1: Normalization and Learning Rate Optimization for BPTT\nWe first systematically investigated the impact of learning rate and batch\nnormalization (BN) on the training performance of RLIF and RadLIF networks. Since\nBPTT is an offline training algorithm capable of processing complete time series,\nfollowing commonly adopted practices, we applied BN operations across both the batch\nand temporal dimensions. Experimental results demonstrate that BN significantly\nimproves the training stability and convergence speed of RLIF networks, thereby\nachieving higher test accuracy (Figure S6).\nHowever, in RadLIF networks, the introduction of BN slightly degraded performance,\nwhich we hypothesize may be due to nonlinear interactions between BN and the adaptive\nneuronal dynamics.\nRegardless of whether BN was employed, both networks achieved optimal performance\nwhen the learning rate was set to 0.005.\nTherefore, in subsequent experiments, we retained BN for RLIF networks while removing\nit for RadLIF networks, and uniformly adopted a learning rate of 0.005 as the default\nsetting.\nStep 2: Surrogate Gradient Function Optimization for BPTT\nNext, we investigated the impact of different surrogate gradient functions on model\ntraining performance. We selected five representative functions: Boxcar, Gaussian,\nMulti-Gaussian, ReLU, and Sigmoid, and conducted comparative training on both RLIF\nand RadLIF networks.\nThe results demonstrate (Figure S7) that different surrogate functions led to\nmarkedly different learning curves and final accuracies.\nNo. 18 / 24\nOverall, the Sigmoid-type surrogate function exhibited the highest accuracy and most\nstable convergence across both network architectures, and was therefore selected as\nthe default surrogate gradient function for BPTT in subsequent experiments.\nStep 3: Impact of Data Augmentation Strategies on BPTT\nSubsequently, we examined the effect of data augmentation on improving training\nperformance with BPTT.\nRecent research by Thomas et al. (2025) [1] demonstrated that multiple data\naugmentation strategies, including\nglobal temporal shift of input channels,\nrandom redistribution of spikes,\ncompression and stretching of stimulus temporal scales,\nintra-class sample mixing,\nand introduction of multi-delay input replicas,\ncan significantly enhance the generalization performance of SNNs.\nInspired by this work, we selected two representative augmentation methods for our\nexperiments: global temporal shift and intra-class blending.\nAs shown in Figure S8, both strategies consistently improved model performance.\nFor RLIF networks, temporal shift augmentation increased test accuracy by\napproximately 3%, with intra-class blending providing an additional 0.1%–0.5%\nimprovement on top of this; for RadLIF networks, the magnitude of improvement was\nslightly smaller (approximately 1%), but the trend remained consistent, indicating\nthat data augmentation can universally enhance model performance across different\nneuronal dynamics.\nNo. 19 / 24\nStep 4: Impact of Smoothing Factor α on pp-prop\nWhen performing hyperparameter optimization for the pp-prop algorithm, we first\ninvestigated the effect of the smoothing factor on model performance.\nThe results demonstrate (Figure S9) that pp-prop exhibits good robustness across a\nwide range of α values in both RLIF and RadLIF networks.\nSpecifically, RLIF networks achieved optimal performance in the range of = 0.90–\n0.93, while RadLIF networks maintained high accuracy across the range of = 0.90–\n0.98.\nUltimately, we selected = 0.93 (RLIF) and = 0.97 (RadLIF) as the default\nparameters for subsequent experiments, respectively.\nNo. 20 / 24\nStep 5: Impact of Learning Rate on pp-prop\nWe further analyzed the impact of learning rate on pp-prop training dynamics and\nperformance.\nAs shown in Figure S10, learning rate plays a decisive role in model convergence and\nfinal accuracy. Learning rates that are either too low or too high lead to training\ninstability or slow convergence.\nExperiments demonstrate that across both network architectures, a learning rate of\n0.01 consistently achieved the highest average accuracy, and we therefore adopted it\nas the standard setting for all subsequent experiments.\nStep 6: Impact of Surrogate Gradient Functions on pp-prop\nConsistent with the BPTT experiments, we also tested five surrogate gradient\nfunctions (Boxcar, Gaussian, Multi-Gaussian, ReLU, Sigmoid) in pp-prop.\nThe results (Figure S11) demonstrate that the type of surrogate function has a\nsignificant impact on model performance, with the Multi-Gaussian function achieving\nthe highest accuracy in both RLIF and RadLIF networks.\nThis finding indicates that pp-prop exhibits strong sensitivity to the choice of\nsurrogate function, and appropriate selection can significantly enhance learning\neffectiveness.\nNo. 21 / 24\nStep 7: Impact of Data Augmentation Strategies on pp-prop\nFinally, we explored the impact of data augmentation on pp-prop learning performance.\nWe similarly evaluated two augmentation methods: global temporal shift and intra-\nclass blending.\nThe results demonstrate (Figure S12) that temporal shift augmentation significantly\nimproved test accuracy for both RLIF and RadLIF networks, indicating that temporal\nperturbations help the model learn time-invariant features and enhance generalization\ncapability.\nIn contrast, adding intra-class blending on top of this did not yield additional\nimprovements, suggesting that under pp-prop's online training mechanism, sample\nmixing may not further regularize the model.\nNo. 22 / 24\nSummary\nSynthesizing the above results, we have identified optimal hyperparameter\ncombinations for BPTT and pp-prop respectively on the SHD dataset.\nUnder identical network architectures and data augmentation conditions, pp-prop's\nfinal performance is comparable to that of BPTT (Table 1).\nSpecifically, in RLIF networks, BPTT achieved a test accuracy of 94.01%, while pp-\nprop achieved 93.93%; in RadLIF networks, BPTT achieved 94.72%, while pp-prop reached\nan average accuracy of 95.33%, even surpassing some state-of-the-art offline\nalgorithms that incorporate specialized structural designs.\nThese results provide strong evidence that, after systematic parameter optimization,\npp-prop can achieve learning performance comparable to BPTT while maintaining linear\nmemory complexity, thereby providing a solid algorithmic foundation for large-scale\nonline learning.\nReferences\n[1] Nowotny, Thomas, James P. Turner, and James C. Knight. \"Loss shaping enhances\nexact gradient learning with eventprop in spiking neural networks.\" Neuromorphic\nComputing and Engineering 5.1 (2025): 014001.\nComment 7\nAs a side remark: I agree with the other reviewers that the name clash with\nBrainScaleS is very unfortunate. I also don't think BrainScale is a particularly\ngood and descriptive name for the proposed model.\nWe thank Reviewer #1 (prior round) and Reviewer #3 (this round) for the clear and\nconstructive feedback on the system name. We fully agree that BrainScale was too\nclose to BrainScaleS , the well-known neuromorphic hardware platform from Heidelberg,\nand that it was not sufficiently descriptive of our contribution.\nNo. 23 / 24\nTo remove any potential confusion and to better reflect the algorithmic core, we have\nrenamed the system to BrainTrace: eligibility trace-based online learning algorithms\nfor brain dynamics modeling . The new name explicitly signals our method’s foundation,\nand motivations for scalable, efficient, biologically plausible online learning in\nspiking neural networks, while remaining consistent with our broader BrainPy\necosystem.\nWe have replaced \"BrainScale\" with \"BrainTrace\" throughout the title, abstract, main\text, figures, and captions. The implementation is now available on PyPI as\nbraintrace:\npip install braintrace\nand exposed under the braintrace namespace:\nimport braintrace\nbraintrace.ETraceParam\nbraintrace.ETraceAlgorithm\nWe believe \"BrainTrace\" removes the name clash with BrainScaleS and communicates the\ntechnical essence of our work more clearly. Thank you again for prompting this\nimprovement.\nNo. 24 / 24"}, {"round": 3, "reviewer_comments": {"reviewer_1": "The authors have thoroughly revised their manuscript and underlying work in response to the concerns raised by the\nreviewers. They have now included detailed benchmarking comparisons to other relevant online learning algorithms and\nback-propagation through time algorithms that illustrate that their new forward learning method is competitive with previous\nwork.\nThe authors have also taken on board the other concerns and have renamed their software to avoid the unfortunate name\nclash with the BrainScaleS neuromorphic system.\nIn its current form, the work presents a significant contribution in the area of online learning in spiking neural networks, and\nthe manuscript presents a very substantial body of work.\nThe potential flaws in the validation approach (tuning parameters or using batch-norm or not) have been addressed\nthoroughly and to my full satisfaction.\nThe code is publicly available on Github, including scripts to reproduce benchmarking results. I note that the code repos still\nhave the previous BrainScale naming. Ideally, this should be updated to avoid confusion.\nOtherwise, I have no further concerns.\n(Remarks on code availability)\nI have opened the repository for a quick look, but have not had the time to review the code in any detail or check its function.", "reviewer_2": "Thank you for addressing my comments.\nUnfortunately, tables S7 and S8 still contain incorrect notations for the values V_reset, V_rest, and V_th. These have been\ncorrected in the other tables. I would suggest making the omitted corrections. However, I do not believe this requires a new\nreview stage, as it is purely an editorial matter.\n(Remarks on code availability)"}, "author_response": "[Review submitted as PDF attachment - text not found in document]"}]

Decision letter