Online Learning Of Neural Computations From Sparse Temporal Feedback
Year: 2022

Abstract:
Neuronal computations depend on synaptic connectivity and intrinsic electrophysiological properties. Synaptic connectivity determines which inputs from presynaptic neurons are integrated, while cellular properties determine how inputs are filtered over time. Unlike their biological counterparts, most computational approaches to learning in simulated neural networks are limited to changes in synaptic connectivity. However, if intrinsic parameters change, neural computations are altered drastically. Here, we include the parameters that determine the intrinsic properties, e.g., time constants and reset potential, into the learning paradigm. Using sparse feedback signals that indicate target spike times, and gradient-based parameter updates, we show that the intrinsic parameters can be learned along with the synaptic weights to produce specific input-output functions. Specifically, we use a teacher-student paradigm in which a randomly initialised leaky integrate-and-fire or resonate-and-fire neuron must recover the parameters of a teacher neuron. We show that complex temporal functions can be learned online and without backpropagation through time, relying on event-based updates only. Our results are a step towards online learning of neural computations from ungraded and unsigned sparse feedback signals with a biologically inspired learning mechanism.

heading: 1 Introduction
text: In biological and artificial spiking neurons, information processing depends on synaptic connectivity as well as intrinsic variables. How a neuron integrates synaptic inputs over time changes its inputoutput function, leading to a large variety of neural computations [1]. In biology, the input-output function of a neuron is modified not only by synaptic plasticity [2], but also by alterations of the neuron’s electrophysiology, e.g., through changes in the composition of trans-membrane ion channels, so-called intrinsic plasticity [3, 4]. Synaptic interactions influence the amount of neurotransmitter released from the presynaptic neuron or the amount of receptors in the postsynaptic terminal and therefore how strongly neurons are connected to each other. Intrinsic properties of the neuron, on the other hand, affect the neuron’s threshold and resting potential and how strongly postsynaptic potentials are amplified [4].
Correspondingly, the computational properties of an artificial neuron depend on the values of its synaptic weights and intrinsic parameters. For instance, the neural computation performed by a leaky integrate-and-fire neuron [5] or a leaky resonate-and-fire neuron [6] with identical synaptic inputs and synaptic weights but different intrinsic parameters varies strongly within and across neuron models (Figure 1A). Learning in artificial neurons is thus the process of modifying a set of synaptic weights and intrinsic parameters to reach a target input-output function. However, unlike their biological counterparts, most computational approaches to learning in spiking neurons are limited to effects of
1. Department of Experimental Psychology, University of Oxford, Oxford, United Kingdom 2. Bernstein Center for Computational Neuroscience Berlin, Berlin, Germany 4. Institute of Science and Technology Austria, Klosterneuburg, Austria
35th Conference on Neural Information Processing Systems (NeurIPS 2021).
alterations in synaptic connectivity. Here, we derive a gradient-based online learning rule for leaky integrate-and-fire neurons and leaky resonate-and-fire neurons that can learn the synaptic connectivity along their intrinsic parameters from sparse temporal feedback in the teacher-student paradigm.

heading: 1.1 Related Work
text: Applying gradient-based optimisation and the backpropagation algorithm [7, 8] to spiking neurons is difficult since the gradient of spike trains is ill-defined: it is zero everywhere except at spike times, where it is infinite. Methods to overcome the problem of non-differentiable spike trains include using gradients with respect to the probability density of spike-times [9], minimising the temporal distance between output and target spike times [10], temporal smoothing of synaptic-currents [11] or the discrete error signals caused by target-spikes [12, 13, 14] or by temporal gradient accumulation using the adjoint method [15]. In contrast to approaches that are directly optimising output spike times, there are also approaches that optimise the activity of spiking neurons through their network activity [16, 17], firing rate response function [18] or via the auxiliary objective of generating a fixed amount of output spikes for certain input spike patterns [19, 20]. Alternatively, so-called "surrogate gradient" methods replace the non-differentiable threshold function with a smooth surrogate function during the backpropagation of errors [21, 22, 23]. Using surrogate gradients and backpropagation through time [24], multilayer spiking neural networks can be successfully trained on a wide range of supervised learning tasks [25, 26, 27], even when errors are backpropagated at spike times only [28].
Our learning rule is closely related to the work by Gütig and Sompolinsky [19] and Memmesheimer et al. [29]. The tempotron learning rule [19] trains leaky integrate-and-fire neurons to emit an output spike for some input spike patterns, but to remain silent for others. Memmesheimer et al. [29] adapted this learning rule to learn synaptic weights that map input spike trains to a target output spike train. In contrast to the tempotron learning rule, the feedback signal does not indicate the class membership of an input spike pattern but target spike times: on the occurrence of each missed or false output spike the distance to the spiking threshold is minimised or respectively, maximised. If an output spike does not fall within a tolerance window around a target spike, a single parameter update is performed and the training procedure and neuron model are reset to the beginning of the training sequence. Consequently, the target spike train is learned one by one, from the first to the last spike of the target sequence, eventually requiring many replays of the sequence until the algorithm converges.
We extend the learning rule by Memmesheimer et al. [29] by an event-dependent scaling factor, such that synaptic weights and intrinsic parameters can be learned mutually and online, i.e., from a continuous stream of inputs without resetting the training procedure after each error and without relying on a tolerance window. Like Memmesheimer et al. [29], we train without backpropagating errors through time [24]. We hypothesise that a learning signal that indicates when a neuron should spike, paired with a gradient-based online parameter update rule, is sufficient to learn the complex temporal input-output relationships of a neuronal computation.

heading: 2 Methods

heading: 2.1 Teacher-student paradigm
text: Learning in artificial neurons relies on finding the set of parameters that facilitates a certain target computation. One approach to test an algorithm’s learning capability is to test its ability to map randomly sampled input spike patterns to randomly sampled output spike trains (e.g. used in [12, 13, 23] and in [30] to compare algorithms). However, under this testing paradigm, the difficulty of the task is directly dependent on the statistics (e.g. dimensionality and firing rate) of the input spike patterns and target spikes. Therefore, the size and shape of the solution space and thus the difficulty of the task are unknown and strongly fluctuate with the choice of hyperparameters. As a consequence, the set of parameters that solves the problem may be large and thus even suboptimal learning algorithms can still yield good performance. For example, some of the solutions to a learning task may not be anticipated, like networks with biologically implausible high firing rates [27]. Further, using random input-output spike patterns of fixed length is unsuitable for evaluating an online-learning algorithm since, as the length of the random pattern increases, the probability of a solution existing is vanishing to zero.
The teacher-student paradigm [31, 32] circumvents this problem by generating target spike times using a model (i.e. the teacher neuron). The learning task is then to recover the teacher’s parameter that determine its neural computation with a student neuron (e.g. in [29] “Reconstruction of Synaptic Weights”). Like this, the learning problem is well defined, as a single point in parameter space. In our case, learning signals are derived from comparing the output spike train of the teacher neuron with the output spike train of the student (Figure 1B). The teacher-student paradigm has several advantages, e.g., over the test performance in a standardised supervised learning task. First, if the teacher and student neuron have the same architecture, it is guaranteed that a unique solution to the training problem exists. Second, the student’s parameters, and thus performance, can be compared to the solution during and after the training process. Finally, since the teacher is used to generate training samples, an infinite amount of unique and continuous training data can be sampled. In summary, by controlling the solution space of the task, the teacher-student paradigm can be used to precisely test the capability of performing spike-pattern recognition and the validity of assumptions of a learning algorithm.

heading: 2.2 Error signals
text: As the dynamics of the teacher and student neuron evolve over time, three possible events can occur when comparing their output spike trains (Figure 1B). First, the student may generate a spike at the same time as the teacher (a hit). Second, the teacher may generate a spike when the student does not (a miss, denoted by d− = −1). Third, the student may generate a spike when the teacher does not (a false positive, denoted by d+ = 1). Since false positives can be detected by the student neuron without feedback, only missed spikes require error signals between teacher and student neurons that are ungraded, unsigned and therefore "unary". Intuitively, a miss occurs when the student’s membrane potential V is too low and a false positive occurs when the student’s membrane potential is too high. We base learning on this observation and update the student’s parameters θ such that the distance to the threshold is reduced for each miss, but increased for each false positive spike [19, 29].

heading: 2.3 The event-dependent scaling online learning rule
text: Our learning rule consists of a gradient-based parameter update like in Memmesheimer et al. [29] (for a derivation see appendix A) and an additional event-dependent scaling factor λ(∆d±). We thus call it the Event-Dependent Scaling (EDS) rule. More specifically, we use the partial derivative of the membrane potential with respect to a parameter and perform a gradient ascent step on each miss (d− = −1) and a gradient descent step on each false positive (d+ = 1) spike:
∆θ ∝ −λ(∆d±)d± ∂V (t) ∂θ . (1)
Here, the EDS factor λ(∆t) scales parameter updates according to the truncated sigmoidal function
λ(∆d±) = 1000− 1000 exp ( ln(0.5) ( min(∆d±,∆θ)
m
)b) , (2)
where ∆θ is a cutoff parameter, m determines the midpoint and b the slope of the sigmoidal and ∆d± denotes the time since the last parameter update. Thus, the size of the parameter update is event-dependent. In practice, we use ∆θ = 75, m = 500 and b = 4, which is plotted in figure 2A. In contrast to the EDS rule, surrogate gradient methods [23, 26, 27] use a scaling function λ(V (t)) that is based on the current membrane potential. We compare the EDS factor λ(∆d±) to the the SuperSpike surrogate gradient λ(V ) = (β|V − VThr|+ 1)−2 [23]. Using the EDS rule, parameter updates are temporal, stochastic and online. Updates are temporal since the underlying function that we optimise is a function evolving over time i.e. the parameterised function that describes the membrane potential. They are stochastic since updates are performed for a single teacher-student pair and not as averages over multiple initialisations. And lastly, they are online since gradients are calculated for the current point in time only, without propagating errors through time [24] and therefore without taking the recurrent nature of neural dynamics into account (cf. [26]). To speed up learning, we use the Adam optimiser [33] with default parameters (β1 = 0.9, β2 = 0.999, = 1e−8) and per parameter learning rates, which are scaled according to the range from which target values are sampled (ηw = 35e−6, ητs = 7e −4, ητm = 28e −4 , ηVr = 7e
−5 and ηw̃ = 8e
−5, ηb = 15e−6, ηω = 33e−7, ηṼr = 8e −5 ηĨr = 8e −5).

heading: 2.4 Neuron models
text: Throughout the experiments, we use two neuron models: the leaky integrate-and-fire (LIF) neuron [5] and the leaky resonate-and-fire (LRF) neuron [6] with current-based synaptic inputs. Both of the models are based on linear subthreshold dynamics, a fixed spiking threshold and a post spike reset rule: Whenever the linear dynamics reach the threshold value, the neuron elicits an output spike and its dynamics are reset according to the reset rule. For simplicity, the membrane potential V of both models is rescaled such that the distance from the resting potential (V0 = 0) and the spiking threshold Vthr is 1. Input, output and target spikes are modelled as Dirac delta spike trains S(t) = ∑ {t̂} δ(∆t̂), where ∆t̂ = t− t̂ substitutes the time that has elapsed since a spike event t̂, and {t̂} is a shorthand notation for the set of n spike times of a spike train {tk | k = 1, 2, ..., n and tk ≤ t}. The condition tk ≤ t preserves temporal causality. Solving the linear system of differential equations that determines the subthreshold voltage dynamics of the leaky integrate-and-fire (LIF) neuron and adding the spike reset rule (appendix B.2) yields
the spike-response model [34]. It describes how the membrane potential V of neuron j evolves over time:
Vj(t) = N∑ i wij ∑ {t̂i} K(∆t̂i) + (Vr − Vthr) ∑ {t̂j} exp ( −∆t̂j τm ) , (3)
with double exponential voltage kernel
K(∆t) = exp ( −∆t τm ) − exp ( −∆t τs ) , (4)
input spike times t̂i, output spike times t̂j and spiking threshold Vthr. The synaptic weights wij determine how strongly the membrane potential V is perturbed by synaptic inputs and if the perturbation is excitatory or inhibitory. The synaptic time constant τs and membrane time constant τm determine how quickly synaptic currents and the membrane potential return to rest (Fig. 2C). Together with the reset potential Vr, they determine the intrinsic dynamics of the neuron. During training, the intrinsic parameters (τs, τm, Vr) are optimised along with synaptic weights.
The temporal dynamics of the leaky resonate-and-fire (LRF) neuron [6] provide a model for phenomena observed in biological neurons that can not be generated by a LIF neuron [1], most prominently sub-threshold oscillations of the membrane potential [35]. This makes LRFs sensitive to inputs of a certain frequency, whereas LIF neurons are sensitive to the rate of inputs. Solving the system of the two coupled linear differential equations (appendix B.3) that describe the oscillating subthreshold potential Ṽ and the resonant current Ĩ yields
Ṽj(t) = N∑ i wij ∑ {t̂i>t̂j} exp(∆t̂ib) sin(∆t̂iω) + exp(∆t̂jb)(Ṽr cos(∆t̂jω) + Ĩr sin(∆t̂jω)), (5)
with input spike times t̂i and most recent output spike time t̂j (Fig. 2D). Akin to a pendulum, the temporal dynamics of the LRF neuron are perturbed by incoming spikes and scaled by the synaptic weights wij . The frequency of the subthreshold oscillations is determined by ω ≥ 0, and the dynamics return to rest based on the damping factor b ≤ 0 (Fig. 2E). In contrast to the LIF neuron, the synaptic currents are reset after each spike, i.e., the dynamics are set to Ĩr, Ṽr on the phase plane. Consequentially, in equation 5 we sum only over input spikes since the most recent output spike {t̂i > t̂j}. In summary, the learnable parameters of the LRF neuron are the synaptic weights wij , the damping factor b, the frequency of the subthreshold oscillations ω and the reset values Ĩr and Ṽr.

heading: 2.5 Parameter initialisations
text: Input spike trains are homogeneously Poisson distributed. In total, there are 100 synaptic inputs and the ratio between excitatory and inhibitory inputs is 4 : 1 [36]. Excitatory and inhibitory input firing rates differ accordingly, i.e., r+in = 10Hz and r − in = 40Hz.
Synaptic weights are sampled independently from a log-normal distribution [37]. The mean µw and variance σ2w of the distribution are set so that approximately 99% of the cumulative distribution function is ≤ 0.2 and thus below one-fifth of the distance between resting potential and threshold (appendix C). Further, we use the bias term β+ to fit the target output firing rate rout by creating a small imbalance between excitatory and inhibitory weights (Fig. 2F):
wexc ∼ β+κ Lognormal(µw, σ2w) and winh ∼ κ Lognormal(µw, σ2w). (6)
If a weight is > 0.3, it is resampled. Weights are then scaled once with the normalisation factor κ [38, p. 143], in order for the peak amplitude of the postsynaptic potential to be equivalent to the size of the sampled weight, independent of the values of the intrinsic parameters (Fig. 2G). Note, that κ is different for leaky integrate-and-fire and leaky resonate-and-fire neurons (appendix D). To ensure balance between excitation and inhibition, target values are resampled as described in section 2.5 until there is a β+ < 2.5 that reaches the target output firing rate rout. Importantly, the position of excitatory and inhibitory weights is randomised, i.e. weights do not adhere to Dale’s principle [39] during learning and can flip sign.
The teacher neuron parameters determine the target computation and are not altered during the training process. Both the student’s initial values and the teacher’s parameter values are chosen
in the following way: For Integrate-and-fire neurons, the membrane time constant is sampled from τm ∼ U(10ms, 60ms). The synaptic time constant τs is smaller than the membrane time constant and set to τs = τm/4. Reset potentials and target firing rates are sampled from Vr ∼ U(−1.5, 0.9) and rout ∼ U(1Hz, 50Hz) respectively. For Resonate-and-fire neurons, the damping factor−b ∼ U(201/ms, 1201/ms) and the frequency of the subthreshold oscillations ω ∼ U(2Hz, 25Hz) are repeatedly sampled until κ < 4. The two reset parameters Ĩr and Ṽr are both sampled from U(−0.8, 0.8). Target output firing rates are sampled from rout ∼ U(1Hz, 20Hz).

heading: 2.6 Convergence and robustness
text: A simulation’s convergence is determined by the signed relative error between the student’s θ and teacher’s θ′ parameter values:
(θ, θ′) =
∑ i(θi − θ′i)
max(||θ′||2, 0.075) , (7)
where 0.075 in the denominator is included to prevent division by zero. The convergence time is the moment the absolute value of the signed relative error | | falls below and stays below the parameter dependent threshold values Θθ for all parameters. For LIF neurons the target threshold values are Θw = ΘVr = 0.15 and Θτs = Θτm = 0.025 and for LRF neurons Θw̃ = 0.05, Θb = Θω = 0.025 and ΘVr = ΘIr = 0.1.
The robustness of the learning rule is tested by adding temporal noise ξ to the spike time of the teacher signal. The noise is sampled from a Gaussian distribution ξ ∼ bN (µ = 0, σ2ξ )e, where b·e denotes rounding to the nearest integer and σξ is the standard deviation of the noise in ms.

heading: 2.7 Implementation, simulation and code availability
text: All simulations are run for n = 30 different random seeds in discretised time at a temporal resolution of 1ms. Reported times refer to simulated and not real time, i.e. the time required to run a simulation was several orders of magnitude faster. For example, 30 independent simulations of 12,000 minutes of simulated time at 1ms temporal resolution (as in Fig. 3), take less than 30 minutes of real time on a single AMD Ryzen 5950x. C++ code to replicate all simulations and plots is publicly available1 under GPLv3 license and uses the MPL 2.0 licensed Eigen software library v3.3.7 [40].

heading: 3 Results
text: The results are separated into a theoretical and an empirical part. In the theoretical part, we derive the partial derivatives (sec. 3.1) that are required for the parameter update rule (eq. 1) and motivate their event-based versions (sec. 3.2 and derivation in appendix E). In the empirical part, we first show that the EDS rule can recover teacher’s parameters in LIF and LRF neurons (sec. 3.3). Subsequently, we investigate the influence and functionality of the EDS factor and compare it to the SuperSpike surrogate gradient [23] in a lesion study (sec. 3.4). Finally, we test the robustness of the algorithm to temporal noise in the teacher signal (sec. 3.5).

heading: 3.1 Gradients
text: Integrate-and-fire neuron: Without explicit derivation, the partial derivatives of the membrane potential V (eq. 3) with respect to one of the synaptic weights wij , the synaptic time constant τs, membrane time constant τm and the reset potential Vr are
∂Vj(t) ∂wij = ∑ {t̂i} K(∆t̂i), (8)
∂Vj(t) ∂τs = − 1 τ2s N∑ i wij ∑ {t̂i} ∆t̂i exp ( −∆t̂i τs ) , (9)
1https://github.com/lukas-braun/learning-neural-computations
∂Vj(t)
∂τm =
1
τ2m N∑ i wij ∑ {t̂i} ∆t̂i exp ( −∆t̂i τm ) + 1 τ2m (Vr − Vthr) ∑ {t̂j} ∆t̂j exp ( −∆t̂j τm ) , (10)
∂Vj(t) ∂Vr = ∑ {t̂j} exp ( −∆t̂j τm ) . (11)
Resonate-and-fire neuron: Similarly, deriving the partial derivatives for the membrane potential Ṽ (eq. 5) of the leaky resonate-and-fire neuron with respect to a single synaptic weight wij , the damping factor b, the frequency of the subthreshold oscillations ω and the two reset values Ṽr and Ĩr results in
∂Ṽj(t)
∂wij = ∑ {t̂i>t̂j} exp(∆t̂ib) sin(∆t̂iω), (12)
∂Ṽj(t)
∂b = N∑ i wij ∑ {t̂i>t̂j} ∆t̂ie(∆t̂ib) sin(∆t̂iω) + ∆t̂je(∆t̂jb)(Ṽr cos(∆t̂jω) + Ĩr sin(∆t̂jω)),
(13)
∂Ṽj(t)
∂ω = N∑ i wij ∑ {t̂i>t̂j} ∆t̂ie(∆t̂ib) cos(∆t̂iω) + ∆t̂je(∆t̂jb)(Ĩr cos(∆t̂jω)− Ṽr sin(∆t̂jω)),
(14)
∂Ṽj(t)
∂Ṽr = exp(∆t̂jb) sin(∆t̂jω) (15) and
∂Ṽj(t)
∂Ĩr = exp(∆t̂jb) cos(∆t̂jω). (16)

heading: 3.2 Event-based dynamics
text: The closed form solution of the membrane potential dynamics and its partial derivatives comprise sums over input and output spike times. These are computationally expensive and create a large memory footprint, because all past spike times have to be stored and summed over. Since the summands are decaying exponential functions, one approximate solution is to only sum over terms that are still contributing significantly. Here, instead of relying on an approximation, we show that one can replace these sums by factorising the exponentials. As a result, the sums over spike times are replaced by a recurrence relation, which depends on constants that are updated on each spike event only. Using these event-based dynamics, both the computational as well as the memory complexity can be significantly reduced. For all factorised equations and their derivations see appendix E.

heading: 3.3 Joint online learning of synaptic weights and intrinsic parameters
text: We asked whether it is possible to learn a neural computation from sparse temporal feedback. We found that, indeed, the EDS rule can be used to recover teacher’s parameters in LIF and LRF neurons for biologically plausible initial and target parameter values (Fig. 3A-B). The learning dynamics contract but do not converge to the solution i.e., parameter updates fluctuate around the solution without further decay (appendix F). As a consequence, the size of the fluctuations and therefore the final precision of the solution is dependent on the learning rate. Since learning happens online, the student can also recover teachers that vary over time (appendix G).
During training, teacher and student output spike trains become more correlated (Fig. 3C,E). After training, an average of 94%± 3% SE of LIF and 96.3%± 3% SE of LRF output spike times match the teachers’ output exactly. For LIF neurons, 3% of the remaining spikes fall within 1ms before, and 2.3% within 1ms after the target time. For LRF neurons, these values are 1.4% and 1.9%. Which, combined, addresses more than 99% of all spikes. Final paramter fluctuations in LIF neurons, especially in the synaptic weights, are higher than in LRF neurons, possibly because synaptic currents
B A
C D E F
are not reset after output spikes in LIF neurons, making credit assignment more difficult. For both neuron models, the convergence times decrease as the target firing rates increase (Fig. 3D,F). This is expected, since parameter updates occur after each missed or false positive spike and therefore higher firing rates result in more parameter updates in the same amount of training time.
Importantly, synaptic weights and output firing rates deviate from their target values and performance is significantly inclined when restricting training to synaptic weights. For LIF (and LRF) neurons only 18.5% ± 0.7% SE (and 12.9% ± 0.6% SE) of output spikes are direct hits and an additional 25.3% (and 19.8%) fall within the 1ms window around the target spike, addressing for less than 44% (and 33%) of total spikes (appendix H).

heading: 3.4 Influence and necessity of the event-dependent scaling factor
text: If we use the “vanilla” parameter update rule (eq. 1 without λ(∆d±)), synaptic weights vanish, and the intrinsic parameters diverge or settle at a constant but wrong value, i.e., the training fails (appendix I). We designed the EDS factor such that parameter updates are not immediately reversed by the next error signal, i.e., consecutive parameter updates are down-scaled as a function of their temporal proximity. To better understand the effect of the EDS factor, we trained parameters separately, as well as in all possible combinations. To this end, we optimised subsets of the student’s parameters and kept all other parameters fixed to their optimal, target values. In LIF neurons, single parameters, except the reset potential, can be learned reliably without the EDS factor, whereas optimisation of parameter combinations fail (Fig. 4A). Surprisingly, the performance of the SuperSpike surrogate is equal to or worse than the vanilla learning rule (except for the pair w, τm). This is due to local minima, or degenerate solutions in which the firing rate is kept stable while parameters diverge (appendix K), rather than larger than permissible noise levels (section 2.6). When optimising synaptic weights only, the recovery rate increases as, through changing β, the surrogate is getting more similar to the vanilla learning rule (Fig. 4B). A similar trend can be observed for LRF neurons, where the performance of surrogate gradients falls below the vanilla performance (Fig. 4C,D plot with all parameter combinations in appendix J) and the performance of surrogate gradients increases the more similar it becomes to and settles at the optimal performance of vanilla learning. Our learning rule performs well for LRF neurons across all parameter combinations. In fact, including more parameters is helping in situation where a single parameter would get stuck in a local minima (data not shown).
Further insights into why the EDS factor is so effective, can be acquired by analysing expected parameter updates. For example, if the target and output spike trains of a LIF neuron are uncorrelated Poisson spike trains, the expected vanilla parameter updates with respect to the reset potentialE[∆Vr] are accurately positive if the output firing rate is too low but negative if the output firing rate is too high (Fig. 4E). However, as training progresses, output and target spike trains become highly correlated. Accordingly, if we calculate E[∆Vr] with the vanilla learning rule for an output-spike train that is a shifted version of the target spike train, i.e. highly correlated, the average parameter updates are distorted (Fig. 4F, left). They linearly increase for positive offsets but stay around zero for negative offsets, leading to inaccurate parameter updates. In contrast, the EDS rule disentangles E[∆Vr], such that too early spikes lead to positive updates and too late spikes to negative updates of equal size, accurately leading to a respectively reduced or increased reset potential (Fig. 4F, right).

heading: 3.5 Influence of temporal noise in the teacher signal
text: Adding temporal, Gaussian noise with standard deviation σξ to the target spike times decreases the precision of the solution found by the student as noise levels increase, but does generally not brake the algorithm (Fig. 4G,H). However, LRF neurons with target firing rates above 35Hz start to diverge at the highest noise level (data not shown). This is not surprising, when we consider that the average inter-spike-interval at 35Hz is approximately 29ms and temporal noise at σξ = 3 frequently yields perturbations > 5ms, rendering noise indistinguishable from signal. Considering the lower average target firing rates in LRF neurons in comparison to LIF neurons, also explains why higher noise levels yield to, on average, similar precision in LRF comparing to LIF neurons at varying noise levels. More generally, adding temporal noise to the feedback signal reveals a subtle relation between firing rates, noise variance and the resulting permissible signal to noise ratio.

heading: 4 Discussion
text: Learning of input-output functions from a continuous stream of information is difficult. We proposed a gradient-based update rule capable of learning the unique parameters that determine a neural computation from ungraded, unsigned and temporally sparse feedback signals. Importantly, EDS learns online and without backpropagating errors backward in time, i.e., without explicit memory traces of input or target spike times, or past states of the neuron. Instead, we derived an efficient and exact implementation of membrane potential dynamics and their gradients. In contrast to previous theoretical work, but in alignment with observations in biological neurons [3, 4], we regard intrinsic parameters like time constants etc. as crucial components of learning. Therefore, in line with recent theoretical work suggesting that networks benefit from heterogeneous neural computations [41, 42], we optimise the intrinsic parameters along the synaptic weights, and further, derive the learning rule for two neuron models that implement different computations: the frequently used and optimised LIF
neuron and the LRF neuron. Whether EDS can be extended to non-linear models like quadratic [43] or exponential LIF [44] or Izhikevich [1] neurons must be investigated.
Numerical simulations show that our algorithm can learn exact neural computations across a wide range of biologically plausible target values. The EDS rule can adapt to changing target computations mid-flight and is robust against temporal noise in the feedback signal. Optimising different combinations of parameters revealed the crucial role played by the EDS factor λ(∆d±). Surprisingly, the EDS factor cannot be replaced by a surrogate gradient. While surrogate gradients are the current state of the art to train spiking multi-layer neural networks with backpropagation through time, they may not be optimal to learn neural computations that go beyond rate coding. More specifically, the observation that surrogate gradients produce diverging parameter values but stabilised firing rates suggests that mutual online learning of intrinsic parameters and synaptic weights from target spike times requires a time-dependent scaling factor in contrast to the voltage-dependent scaling of surrogate gradients.
It would be interesting, but brake the online learning assumption, to replace surrogate gradients with EDS factors in multi-layer networks and train them with backpropagation through time. An online learning alternative could use EDS to extend difference target propagation [45, 46] from rate-coded neurons to spiking neural networks. In such a network, the parameters of neurons that provide the downstream signals are optimised to predict spike times of neurons in lower layers, and parameters of upstream neurons are updated such that the backpropagated target spike times are generated.

heading: 4.1 Limitations and conclusions
text: While we motivate the EDS factor λ(∆d±) by means of its empirical success, a rigorous theoretical explanation for its function is missing. An analytical derivation of the expected gradient values for correlated output and target spike trains in spirit of figure 4 E and F could further unveil the underlying principles of our algorithm.
In LRF neurons target computations with sustained subthreshold oscillations (−b 201/ms) or spike bursts, i.e., dynamics that reset outside of the separatrix, cannot be learned reliably (data not shown). In addition, target firing rates in LRF neurons remain low, due to the difficulty of achieving high output firing rates from uncorrelated Poisson inputs while keeping a balance between excitation and inhibition. If our algorithm also works for correlated, oscillating or bursting inputs near the neuron’s eigenfrequency, which cause high output firing rates, is left for future studies.
Further, convergence to target values is slow (up to ≈ 5.5 and ≈ 4.3 days of simulated time in LIF and LRF neurons respectively). As a proof of principle, we started each simulation at a random point in parameter space and required convergence to a second, randomly sampled point in parameter space. This could lead to long learning trajectories and slow convergence times (see also [47]), especially when target firing rates were low. In contrast, the parameter subspace that results in good performance in supervised rate-coded networks is large [48], and spiking neural networks may solve a task with a wide range of different architectures or parameters [27]. From a biological perspective, most electrophysiological changes induced by plasticity mechanisms are gradual and a complete reorganisation of a neuron’s electrophysiology (as in our study) implausible. Ultimately, training speed could be increased dramatically (>10x) by increasing learning rates at the cost of losing precision (appendix F) or by using a learning rate schedule to speed up initial learning and subsequently gain precision as the learning rate decays. Convergence time could be further reduced by batch-learning and when modelling biological data, by using educated guesses for initial model parameters (e.g. pick parameters such that the initial firing rate of the model and data are identical).
We conclude that the complex temporal dynamics of neural computations can be learned online with a local gradient-based update rule from feedback signals that indicate target spike times. Our eventdependent scaling rule will help to bring the success of gradient-based optimisation to spiking neural networks performing heterogeneous computations and contribute to the exploration of gradient-based learning in the brain.
Acknowledgments and Disclosure of Funding
We would like to thank Professor Dr. Henning Sprekeler for his valuable suggestions and Dr. Andrew Saxe, Milan Klöwer and Anna Wallis for their constructive feedback on the manuscript. Lukas Braun was supported by the Network of European Neuroscience Schools through their NENS Exchange Grant program, by the European Union through their European Community Action Scheme for the Mobility of University Students, the Woodward Scholarship awarded by Wadham College, Oxford and the Medical Research Council [MR/N013468/1]. Tim P. Vogels was supported by a Wellcome Trust Senior Research Fellowship [214316/Z/18/Z].

Review summaries:
This paper introduces the event-dependent scaling (EDS) rule to use in neural network training that allows neural computations to be learned online - i.e. without error propagation. It presents a teacher-student paradigm for problem formulation, then characterizes the error signal and EDS rule. It then describes its architectures (leaky integrate-and-fire or leaky resonate-and-fire neurons, Poisson-distributed input spike trains, and other standard parameters), and metrics (of convergence and robustness). Next, it discusses resulting learning dynamics and analyses the effectiveness of EDS empirically, and concludes by acknowledging concerns behind lack of theoretical backing and a rarefied testing paradigm.