Supervised training of spiking neural networks for robust deployment on mixed-signal neuromorphic processors

Mixed-signal analog/digital circuits emulate spiking neurons and synapses with extremely high energy efficiency, an approach known as “neuromorphic engineering”. However, analog circuits are sensitive to process-induced variation among transistors in a chip (“device mismatch”). For neuromorphic implementation of Spiking Neural Networks (SNNs), mismatch causes parameter variation between identically-configured neurons and synapses. Each chip exhibits a different distribution of neural parameters, causing deployed networks to respond differently between chips. Current solutions to mitigate mismatch based on per-chip calibration or on-chip learning entail increased design complexity, area and cost, making deployment of neuromorphic devices expensive and difficult. Here we present a supervised learning approach that produces SNNs with high robustness to mismatch and other common sources of noise. Our method trains SNNs to perform temporal classification tasks by mimicking a pre-trained dynamical system, using a local learning rule from non-linear control theory. We demonstrate our method on two tasks requiring temporal memory, and measure the robustness of our approach to several forms of noise and mismatch. We show that our approach is more robust than common alternatives for training SNNs. Our method provides robust deployment of pre-trained networks on mixed-signal neuromorphic hardware, without requiring per-device training or calibration.


Scientific Reports
| (2021) 11:23376 | https://doi.org/10.1038/s41598-021-02779-x www.nature.com/scientificreports/ In contrast to current mainstream Deep Neural Networks (DNNs), spiking networks suffer from a severe configurability problem. The backpropagation algorithm permits configuration of extremely deep NNs for arbitrary tasks 22 , and is effective also for network models with temporal state 23 , but is difficult to apply to the discontinuous dynamics of SNNs [24][25][26] . Methods to approximate the gradient calculations by using surrogate functions 27 , eligibility traces 28 or adjoint networks 29 have provided a way to adapt backpropagation for spiking networks. Non-local information is required for strict implementation of the backpropagation algorithm, but random feedback 30 and local losses 31 have been employed with some success to train multi-layer spiking networks. Alternative approaches using initial random dynamics coupled with error feedback and spike-based learning rules can permit recurrent SNNs to mimic a teacher dynamical system 32,33 . Strictly-local spike-timing-based learning rules, inspired by results in experimental neuroscience 34 , have been implemented in digital and mixed-signal neuromorphic devices, as they provide a better match to the distribution of information across neuromorphic chips 35 . Unfortunately, local spike-dependent rules such as Spike-Timing Dependent Plasticity (STDP) are themselves not able to perform supervised training of arbitrary tasks, since they do not permit error feedback or error-based modification of parameters. In both cases, implementing strictly local or backpropagation-based learning infrastructure on-chip adds considerable complexity, size and therefore cost to neuromorphic hardware designs. This cost makes it impractical to use on-chip learning and adaptation to solve the mismatch problem on mixed-signal architectures.
Robustness to noise and variability can be approached from the architectural side. For example, a network architecture search approach can identify networks that are essentially agnostic to precise weight values 36 . However, these networks rely on complex combinations of transfer functions which do not map to neuromorphic SNN designs.
Alternatively, a class of analytically-derived network architectures have been proposed for spiking networks, known as Efficient Balanced Networks [37][38][39][40][41][42][43] , relying on a balance between excitation and inhibition to provide robustness to sources of noise including spike-time stochasticity and neuron deletion. These networks can be derived to mimic an arbitrary linear dynamical system through an auto-encoding architecture 38 or can learn to represent and mimic dynamical systems 37,[40][41][42] . We propose to adapt the learning machinery of this spiking architecture to produce deployable SNN-based solutions for arbitrary supervised tasks that are robust to noise and device mismatch.
In this work we present a method for training robust networks of Leaky Integrate and Fire (LIF) spiking neurons that can solve supervised temporal signal regression and classification tasks. We adopt a knowledge distillation approach, by first training a non-spiking Recurrent Neural Network (RNN) to solve the desired supervised task using Back-Propagation Through Time (BPTT) 23 . By then interpreting the activations of the RNN as a teacher dynamical system, we train an SNN using an adaptation of the learning rule from Ref. 41 to mimic the RNN. We show that the resulting trained SNN is robust to multiple forms of noise, including simulated device mismatch, making our approach feasible for deployment on to mixed-signal devices without post-deployment calibration or learning. We compare our method with several other standard approaches for configuring SNNs, and show that ours is more robust to device mismatch.

Results
We assume a family of tasks defined by mappings c(t) →ŷ(t) , where c(t) ∈ R d1 and ŷ(t) ∈ R d2 are temporal signals with arbitrary dimensionality ( Fig. 1a; see "Methods"). For simplicity of notation we do not write the temporal dependency "(t)" for the remainder of the paper. This definition encompasses any form of deterministic temporal signal processing or classification task without loss of generality. We refer to our network architecture as ADS (Arbitrary Dynamical System) spiking networks.
Our approach begins by training a non-spiking rate network to implement the arbitrary task mapping by learning the dynamical system through modification of the recurrent weights ˆ ∈ RN ×N ; encoding and decoding weights F ∈ R d1×N and D ∈ RN ×d2 ; biases b ∈ RN ; time constants τ ∈ RN ; and non-linear transfer function f (·) = tanh(·) . BPTT or any other suitable approach can be used to obtain the trained rate network.
We subsequently train a network of spiking neurons to emulate x , with leaky membrane dynamics defined by with spike trains o = V > V thresh produced when exceeding threshold voltages V thresh ; leak rate ; and fast and slow recurrent weights f and s ( Fig. 1b; see "Methods"). The decoded dynamics x ≈x are obtained from the filtered spiking activity r with x = Fr . By feeding back an error signal e =x −x under the control of a decaying feedback rate k, the spiking network is forced to remain close to the desired target dynamics. f is initialised to provide fast balanced feedback 39 , and s is learned using the rule under learning rate η (see "Methods" and Ref. 41 ). Note that we do not require complex multi-compartmental neurons or dendritic nonlinearities in our neuron model, but use a simple leaky integrate-and-fire neuron that is compatible with compact mixed-signal neuromorphic implementation 2 . Once the spiking network has learned to represent x ≈x with high accuracy, we replace the rate network entirely with the spiking network (Fig. 1c). see "Methods"). This task requires memory of past inputs to produce a delayed output, as well as a nonlinear mapping between the memory state and the output variable. A network receives a single input channel where pulses of varying width (100-230 ms) and sign are presented in sequence. The network must report the XOR of the two input pulses by delivering an output pulse of appropriate sign after the second of the two input pulses. A non-spiking RNN ( N = 64 ) was trained to perform the temporal XOR task, using BPTT with Mean-Squared Error (MSE) loss against the target output signal (target and output signals shown in Fig. 2a). After 20 epochs of training with 500 samples per epoch, the RNN reached negligible error on 200 test samples ( ≈ 100% accuracy). A spiking ADS network ( N = 320 ) was then trained to perform the task, reaching equivalent accuracy ( Fig. 2a,b).
Wake-phrase detection. The temporal XOR task demonstrates that one-dimensional nonlinear tasks requiring memory can be learned through our method through supervised training. To show that our approach also works on more realistic tasks with complex input dynamics, we implemented an audio wake-phrase detection task ( Fig. 3; see "Methods"). Briefly, real-time audio signals were extracted from a database of spoken wake phrases ("Hey Snips" dataset 44 ), or from a database of noise samples ("DEMAND" dataset 45 ). The target wake phrase data was augmented with noise at an SNR of 10 dB, then passed through a bank of 16 Butterworth filters with central frequencies spaced between 0.4 and 2.8 kHz (Fig. 3b). We trained a non-spiking RNN ( N = 128 ) to perform the task with high accuracy, using BPTT under an MSE loss function against a smooth target classification signal (Fig. 3d,e). We then trained a spiking ADS network ( N = 768 ) to implement the audio classification task. The non-spiking RNN achieved a testing accuracy of ≈ 90% , and our spiking imitator achieved ≈ 87% after training for 10 epochs on 1000 training samples.
Training considerations. We found that slower input, internal and target dynamics in the RNN were easier for the SNN to reconstruct than very rapid dynamics, depending on the neuron and synaptic time constants in the SNN. Longer and slower target responses yielded smoother ANN dynamics, which were easier for the spik- www.nature.com/scientificreports/ ing ADS network to learn. Our approach did not assume any dendritic non-linearities, or multi-compartmental dendrites with complex basis functions. Instead, the non-linearity of the spiking neuron dynamics is sufficient to learn the dynamics of a non-spiking ANN using the tanh nonlinearity. We found that including a learning schedule for the error feedback rate k was important to achieve low reconstruction error. The factor k must drop to close to zero before the end of training, or else the SNN learns to rely on error feedback for accuracy, and generalisation will be poor once error feedback is removed. Conversely, if k drops too rapidly during training, the SNN is not held close to the desired target dynamics, and is unable to correctly learn the slow feedback weights s . For these reasons, a well-chosen schedule for k is important during learning. In this work we chose a progressive stepping function that decrements k by a fixed amount after some number of signal iterations (see "Methods"). Setting k to a fixed value for some number of trials enables the SNN to adapt to the corresponding scale of error feedback by updating s .
Robustness to noise sources. The slow learned recurrent feedback connections s in the spiking network enable the SNN to reproduce a learned task. In contrast, the balanced fast recurrent feedback connections f are designed to enable the SNN to encode the dynamic variables x in a way that is robust to perturbation 38,39 . We examined the robustness of our trained networks to several sources of noise (Fig. 4).
Device mismatch. We first introduced frozen parameter noise as a simulation of device mismatch present in mixed-signal neuromorphic implementation of event-driven neuron and synapses. We measured distributions of neuronal and synaptic parameters induced in silicon spiking neurons by device mismatch (see "Methods"; Fig. S1). Measurements were performed on 1 core of 256 analog neurons and synapses, on fabricated mixedsignal neuromorphic DYNAP-SE processors 2 . We observed a consistent relationship between the mean and variance of parameter distributions: the variance of the measured parameters increased linearly with the magnitude of the set parameter. We used this experimentally-recorded relationship to simulate mismatch in our spiking network implementations, simulating deployment of the networks on mixed-signal neuromorphic hardware. Mismatched parameters ′ were generated with � ′ ∼ N (�, δ�) , where δ determines the level of mismatch, which we found experimentally to be between 10-20%. Under 20% simulated mismatch on weights, thresholds, biases, synaptic and neuronal time constants, our networks compensated well for the frozen parameter noise present in mixed-signal deployment (Fig. 4b).
Quantisation noise. In contrast to 64-bit floating point precision used by the non-spiking RNN, deployment of NN architectures in memory-constrained systems often uses low bit-depth precision for weights and neuron state. Mixed-signal neuromorphic architectures use analog voltages or currents to represent internal neural state, We imposed weight quantisation constraints on our spiking model, and found that our networks compensated well for the resulting frozen quantisation noise ( Fig. 4c; see "Methods").
Thermal noise. Due to the analog representation of neuron and synapse states in mixed-signal neuromorphic chips, these state variables are subject to thermal noise. Thermal noise appears as white-noise stochastic fluctuations of all states. We simulated thermal noise by adding noise ζ ∼ N (0, σ ) to membrane potentials V, with σ = 1%, 5%, 10% scaled to the range between reset and threshold potentials V reset and V thresh . The spiking ADS network performed well in the presence of thermal noise (Fig. 4d).
Sudden neuron failure. The fast recurrent feedback connections f present in spiking balanced networks have been shown to be able to compensate for neuron loss, where a subpopulation of spiking neurons is silenced during a trial 38,41,46 . We examined this property in our spiking ADS networks that include fast balanced feedback, and found that indeed our networks compensated well for neuron loss (Fig. 4e). In the absence of fast recurrent feedback (i.e. f = 0 ), neuron silencing degraded the performance of the spiking ADS networks (Fig. S3).  implementations of arbitrary tasks, defined through supervised training. We compared our approach against several alternative methods for supervised training of SNNs, and evaluated the performance of these methods under simulated deployment on mixed-signal neuromorphic hardware: • Reservoir Computing, in the form of a Liquid State Machine 18 , relies on the random dynamics of an SNN to project an input over a high-dimensional temporal basis. A readout is then trained to map the random temporal basis to a specified target signal, using regularised linear regression. Since perturbation of the weights and neural parameters will directly modify the temporal basis, we expect the Reservoir approach to perform poorly in the presence of mismatch. • The spiking FORCE algorithm 32 trains an SNN to mimic a teacher dynamical system. We applied this algorithm to a trained non-spiking RNNs to produce a trained SNN, similarly as in our spiking ADS approach. • We implemented the BPTT algorithm to train an SNN end-to-end, using a surrogate gradient function similar to Ref. 25 . During training, these networks received input and target functions identical to those presented to the non-spiking RNN.
We first examined simulated deployment of all architectures by simulating parameter mismatch ( Fig. 5; see "Methods"). We trained 10 networks for each architecture, and evaluated each network at three levels of mismatch ( δ = 5%, 10%, 20% ) for 10 random mismatch trials of 500 samples each. We quantified the effect of mismatch on the performance of each network architecture by measuring the MSE between the SNN-generated output ỹ and the training target for that architecture.  www.nature.com/scientificreports/ a direct test of inherent robustness to quantisation noise. The Reservoir architecture broke down for any quantisation level (chance task performance accuracy ≈ 50% ). The FORCE architecture performed well down to 5 bits (median accuracy 85%), beyond which MSE increased and performance decayed to chance level at 3 bits (med. acc. 50%). Both the ADS spiking network and BPTT architectures maintained good performance down to 4 bits of precision (med. acc. ADS 81%; BPTT 87%), decaying to chance level at 2 bits (med. acc. ADS 52%; BPTT 54%). We compared the effect of thermal noise of the four architectures, simulated as membrane potential noise ( Fig. S4; see "Methods"). The FORCE architecture was most robust to thermal noise, performing best at all noise levels (higher accuracy, p < 5 × 10 −2 ; lower MSE, p < 1 × 10 −3 except for highest level of noise; U test). All other architectures degraded progressively with increasing noise levels. Our spiking ADS network architecture showed the smallest degradation in general over increasing noise levels (MSE 0.0083-0.115; acc. 82-67%). The BPTT architecture also fared well, while dropping in accuracy for the largest noise level (med. acc. 56%).

Power comparison for mixed-signal and traditional implementations.
We estimated and compared the power requirements between a direct implementation of the recurrent non-spiking network dynamics on commodity and ASIC hardware, against our mixed-signal spiking implementation of the network dynamics. We performed the power comparison for the real-time audio processing task outlined above, for varying recurrent network dimensions. Computation on the DYNAP-SE1 processor occurs continuously in real-time, with no clock. We selected the slowest clock speeds for the commodity hardware that are sufficient to support real-time operation.
We estimated the power requirements for an ultra-low-power digital microcontroller from ST Microelectronics (STM32L552xx) 47 (see Table S1). When operating at a 16 MHz clock frequency and efficiently implementing only the recurrent dynamics required by a N = 64-neuron non-spiking RNN, the low-power MCU was estimated to require 260 µW when simulating with a time-step dt = 10 ms , increasing to 1130 µW for dt = 1 ms . For the equivalent spiking network with N = 768 spiking neurons the DYNAP-SE1 processor requires 288 µW when fabricated at 180 nm process, and 38 µW when fabricated at 65 nm process. For larger non-spiking RNNs, the DYNAP-SE1 processor has an increasing energy advantage over the low-power MCU.
We also considered the implementation of the non-spiking RNN on an ultra-low-power ASIC, EIE 48 . When implementing the dynamics required by a N = 64-neuron non-spiking RNN, the ASIC required 11 µW when simulating with a time-step dt = 10 ms , increasing to 105 µW for dt = 1 ms . The ASIC displays a power advantage when simulating dynamics for extremely small RNNs with N < 35 , or with large time-steps dt = 10 ms and N < 200 . For larger networks and with more accurate temporal dynamics, the mixed-signal SNN implementation using our approach is more energy-efficient. For further details of the power estimations see "Methods" and Table S1.

Discussion
We propose a method for supervised training of spiking neural networks that can be deployed on mixed-signal neuromorphic hardware without requiring per-device retraining or calibration. Our approach interprets the activity of a non-spiking RNN as a teacher dynamical system. Using results from dynamical systems learning theory, our spiking networks learn to copy the pre-trained RNN and therefore perform arbitrary tasks over temporal signals. Our method is able to produce spiking networks that perform both simple and complex nonlinear temporal detection and classification tasks. We show that our networks are considerably more robust to several forms of parameter and state noise, compared with several other common techniques for training spiking networks. Our networks are by design robust to common sources of network and parameter variation, both intra-and inter-chip, which must be compensated for when deploying to mixed-signal neuromorphic hardware. For levels of mismatch measured directly from neuromorphic devices, we show that common SNN network architectures break down badly. Usual approaches for compensating for mismatch-induced parameter variation on neuromorphic hardware employ either on-device training [49][50][51][52][53] or per-device calibration 14,15,[53][54][55] , entailing considerable additional expense in hardware complexity or testing time. In contrast, our method produces spiking networks that do not require calibration or retraining to maintain performance after deployment. As a result, our approach provides a solution for cost-efficient deployment of event-driven neuromorphic hardware.
The coding scheme used by our spiking networks has been shown to promote sparse firing 38 . For mixed-signal neuromorphic hardware, power consumption is directly related to the network firing rate. Our method therefore produces networks that consume little power compared with alternative architectures that use firing-rate encoding or do not promote sparse activity 18,23,32 .
Our approach to obtain high-performing SNNs is at heart a knowledge-transfer approach, relying on copying the dynamics of a highly-performing non-spiking RNN. This two-step approach is needed because the learning rule for our SNN requires a task to be defined in terms of a dynamical system, and is not able to learn the dynamics of an arbitrary input-output mapping (see Supplementary Methods). Consequently, our spiking networks can only perform as well as the pre-trained non-spiking RNN, and require multiple training steps to build a network for a new task. Nevertheless, training non-spiking RNNs is efficient when using automatic differentiation, justin-time compilation and automatic batching 56 , and can be performed rapidly on GPUs. Our approach trades off between training time on commodity hardware, and immediate deployment on neuromorphic hardware with no per-device training required.
The robustness of our spiking ADS networks comes partially from the fast balanced recurrent feedback connections, which ensure sparse encoding and compensate in real-time for encoding errors 38,39 . These weights also degrade under noise, but can be adapted in a local untrained fashion using local learning rules that are compatible with HW implementation 57 .
Our supervised training approach is designed for temporal tasks, where input and target output signals evolve continuously. This set of tasks encompasses real-time ML-based signal processing and recognition, but is a poorer fit to high-resolution frame-based tasks such as frame-based image processing. These "one-shot" tasks can be mapped into the temporal domain by serialising input frames 58 or by using temporal coding schemes 59 . We found anecdotally that temporal discontinuities in input and target time series made training our ADS networks more difficult, with the implication that a careful matching between task and network time constants is important.
Our approach builds single-population recurrent spiking networks, in contrast to deep non-recurrent network architectures which are common in 2021 60 . Recurrent spiking networks such as Liquid State Machines (LSMs) have been shown to be universal function approximators 61 , but RNNs do not perform the progressive task decomposition that can appear in deep feed-forward networks 62 . Interpretability of the internal state of recurrent networks such as ours is therefore potentially more difficult than for deep feedforward architectures.
Neuromorphic implementation of spiking neural networks has been hailed as the next generation of computing technology, with the potential to bring ultra-low-power non-von-Neumann computation to embedded devices. However, parameter mismatch has been a severe hurdle to large-scale deployment of mixed-signal neuromorphic hardware, as it directly attacks the reliability of the computational elements -a problem that commodity digital hardware generally does not face. Previous solutions to device mismatch have been impractical, as they require expensive per-device calibration or training prior to deployment, or increased hardware complexity (and therefore cost) in the form of on-device learning circuits. We have provided a programming method for mixed-signal neuromorphic hardware that frees application developers from the necessity to worry about computational unreliability, and does not require per-device handling during or after deployment. Our approach therefore removes a significant obstacle to the large-scale and low-cost deployment of neuromorphic devices.

Methods
We trained and simulated ANNs and SNNs using Rockpool 63 , an open-source Python package for machine learning of SNNs. We implemented a liquid state machine SNN 18 ; spiking FORCE network 32 ; and a BPTT-trained SNN 23 using Jax 56 and custom-written forward-Euler solvers. Parameters for all architectures are given in the Supplementary Material. Code to generate all models, analysis and figures in this paper are available from https:// github. com/ synse nse/ Robust-Class ifica tion-EBN.
Temporal XOR task. We created signals of a total duration of 1 second, of which the first two thirds were dedicated to the input and the last third to the target (Fig. 2). During the input time-frame, two activity bumps were created on a single input channel representing the binary inputs to the logical XOR operation. The bumps had varying length (uniformly drawn between 66-157 ms) and magnitude ±1 , and were smoothed with a Gaussian filter to produce smooth activity transitions. In the final third of the signal we defined a target bump of www.nature.com/scientificreports/ magnitude ±1 , indicating the true output of the XOR operation. The target bump was also smoothed with a Gaussian filter. We trained a rate network ( N = 64 ) to high performance on the XOR task, then subsequently trained a spiking model ( N = 320 ) to follow the dynamics of the trained rate network . We used a fixed learning rate η = 1 × 10 −5 and fixed error feedback rate k = 75 during SNN training. Output classification from both networks was determined by the network output passing the thresholds ±0.5.
Speech classification task. We drew samples from the "Hey Snips" dataset 44 , augmented with noise samples from the DEMAND dataset 45 , with a signal-to-noise ratio of 10 dB . Each signal had a fixed length of 5 s and was pre-processed using a 16-channel bank of 2nd-order Butterworth filters with evenly-spaced centre frequencies ranging 0.4-2.8 kHz. The output of each filter was rectified with abs(·) , then smoothed with a 2nd order Butterworth low-pass filter with cut-off frequency 0.3 kHz to provide an estimate of the instantaneous power in each frequency band. The rate network for the speech classification task ( N = 128 ) was trained for 1 epoch on 10 000 samples to achieve roughly the same performance as the spiking network trained with BPTT. We trained spiking networks ( N = 768 ; τ mem = 50 ms ; τ fast = 1 ms ; τ slow = 70 ms ) for 5 epochs on 1000 training samples, validated on 500 validation samples and 1000 test samples. To perform a classification we integrated the output of the network when it passed a threshold of 0.5. We then applied a subsequent threshold on this integral, determined by a validation set, to determine the final prediction. We used a fixed learning rate η = 1 × 10 −4 and a decaying step function for the error feedback factor k (from 200-25 in 8 evenly-spaced steps).

Spiking neuron model and initialisation.
We used an LIF neuron model with a membrane time constant τ mem = 50 ms ; reset potential V reset = 0 ; resting potential V rest = 0.5 ; and spiking threshold V thresh = 1 . The membrane potential dynamics for the neuron model were given by with input current I inp ; fast and slow recurrent post-synaptic potentials (PSPs) I fast and I slow ; error current I e = kD T e ; and noise current η . Output spikes from a neuron are given by o(t) = V > V thresh . Synaptic dynamics were described by with input synaptic weights W; synaptic time constants τ syn = 1 ms and 70 ms for fast and slow synapses, respectively; and simulation time step t . Feed-forward and decoding weights were initialised using a standard normal distribution scaled by the number of input/output dimensions ( N ). Fast balanced recurrent feedback connections were initialised and rescaled according to the threshold and reset potential, as described in Ref. 38 . The spiking network was simulated using a forward Euler solver with a simulation time step of 1 ms.
Non-spiking network. The dynamics of a neuron in the non-spiking RNN were described by with input c(t); encoding weights F ; recurrent weights ˆ ; non-linearity f (·) = tanh(·) ; bias b; and noise term ǫ . Time constants τ were initialised with linearly spaced values ( 10 − 100 ms ). The trainable parameters in this network are the time constants τ ; the encoding and recurrent weights F and ˆ ; and the biases b. No noise was applied during training or inference ( ǫ = 0).

Measurements of parameter mismatch.
Using recordings from fabricated mixed-signal neuromorphic chips we measured levels of parameter mismatch (i.e. fixed substrate noise pattern) present in hardware. In particular, for DYNAP-SE 2 , a neuromorphic processor which emulates LIF neuron, AMPA and NMDA synapse models with analog circuits, we measured neuron and synaptic time constants, and synaptic weights for individual neuron units, by recording and analysing the voltage traces produced by these circuits. We observed levels of mismatch in the order of 10-20% for individual parameters, with widths of the distributions being proportional to the means (see Fig. S1).
Power estimates. Since the input and output weighting differs between the spiking and non-spiking network, and comprises only a small portion of the parameters, we limited our power comparison to the recurrent portion of the network. Updating the recurrent dynamics for the non-spiking rate network requires multiply-accumulate operations for the recurrent input r t =�f (x t ) (neglecting the transfer function f (·) ); multiply-accumulate operations for the Euler solver update x t+1 = x t +ẋ t * dt/τ ; and accumulate operations for ẋ t = −x t + i t + b + r t . With N = 64 neurons, these amount to 8576 OPs, with MACs counted as two OPs. With a time-step of dt = 1 ms , this corresponds to 8.58 GOPS (Giga-OPs per second). We estimated the power to implement our RNN on non-neuromorphic NN accelerators by using previously reported power as GOPS/W. We examined only chips with published data for total power, and where we could identify the fabrication node for the published chip. We re-scaled power estimates to normalise against the fabrication node, providing estimates for 65 nm nodes in all cases. For the ultra-low-power microcontroller (STM32L552xx), we assumed that the MCU switched to a low-power sleep mode once the dynamics for a given time-step were computed. This τ mem ∂V ∂t = V rest − V + I inp + I fast + I slow + I e + η n τ syn ∂I * ∂t = −I * + (Wo(t)τ syn )/�t τ jẋj = −x j +Fc j (t) +�f (x) + b j + ǫ j Scientific Reports | (2021) 11:23376 | https://doi.org/10.1038/s41598-021-02779-x www.nature.com/scientificreports/ permits the MCU to save power when only a portion of computing resources is required to simulate real-time dynamics. Again neglecting synaptic operations required for input and output, we estimate the energy for routing a single recurrent spike on the DYNAP-SE1 mixed-signal neuromorphic processor as 3.3 nJ . We found that the firing rate of the spiking population is upper-bounded by approximately 3 Hz per neuron during simulation. For the spiking recurrent population with N = 768 , this corresponds to energy usage of 7.6 µW dynamic power consumption. Static power consumption for the DYNAP-SE1 processor is estimated at 30 µW . Table S1 compares the energy consumption of running the ANN on an efficient ASIC 48 , and a low-power general purpose MCU 64 to the energy consumption of the DYNAP-SE1 using the spiking network with 12 times more neurons.
Simulated mismatch. To simulate parameter mismatch in mixed-signal neuromorphic hardware we derived a model where the values for each parameter follow a normal distribution with the standard deviation depending linearly on the mean. The mismatched parameters ′ are obtained with � ′ ∼ N (�, δ�) where δ determines the level of mismatch. We considered three levels of mismatch: 5, 10 and 20%.
Simulated thermal noise. Thermal noise is inherent in neuromorphic devices and can be modeled by Gaussian noise on the input currents. We applied three different levels of thermal noise ( σ = 0.01, 0.05, 0.1 ) that was scaled according to the difference between V reset and V thresh to assure equal amounts of noise for neuron model and network architecture.
Neuron silencing. We created four network instances grouped into two pairs: One pair was trained with the fast recurrent feedback connections f as described above, and the other pair with f = 0 . We then clamped 40% of the neurons of one instance of both pairs to V reset while evaluating 1000 test samples.

Benchmark network architectures.
We investigated the robustness to simulated noise for four different learning paradigms, including the FORCE method 32 , BPTT 23 and reservoir computing 18 . All parameters for these methods are given in Supplementary Material.

Statistical tests.
All statistical comparisons were double-sided Mann-Whitney U tests unless stated otherwise.

Data availability
Code to generate all models, analysis and figures in this paper are available from https:// github. com/ synse nse/ Robust-Class ifica tion-EBN.