# Timing-Based Backpropagation in Spiking Neural Networks Without Single-Spike
Restrictions

Kakei Yamamoto, Yusuke Sakemi, Kazuyuki Aihara Kakei Yamamoto is with the
University of Tokyo, Tokyo, Japan (e-mail: kakei@g.ecc.u-tokyo.ac.jp).Yusuke
Sakemi is with Research Center for Mathematical Engineering, Chiba Institute
of Technology, Narashino, Japan.Kazuyuki Aihara is with the International
Research Center for Neurointelligence (WPI-IRCN), The University of Tokyo
Institutes for Advanced Study, The University of Tokyo, Tokyo, Japan.

###### Abstract

We propose a novel backpropagation algorithm for training spiking neural
networks (SNNs) that encodes information in the relative multiple spike timing
of individual neurons without single-spike restrictions. The proposed
algorithm inherits the advantages of conventional timing-based methods in that
it computes accurate gradients with respect to spike timing, which promotes
ideal temporal coding. Unlike conventional methods where each neuron fires at
most once, the proposed algorithm allows each neuron to fire multiple times.
This extension naturally improves the computational capacity of SNNs. Our SNN
model outperformed comparable SNN models and achieved as high accuracy as non-
convolutional artificial neural networks. The spike count property of our
networks was altered depending on the time constant of the postsynaptic
current and the membrane potential. Moreover, we found that there existed the
optimal time constant with the maximum test accuracy. That was not seen in
conventional SNNs with single-spike restrictions on time-to-fast-spike (TTFS)
coding. This result demonstrates the computational properties of SNNs that
biologically encode information into the multi-spike timing of individual
neurons. Our code would be publicly available.

###### Index Terms:

Spiking neural networks (SNNs), supervised learning, backpropagation, temporal
coding, multi-spike.

##  I Introduction

How does the human brain acquire intelligence? Neural systems in the brain
have various functions by forming very complex networks, and much remains a
mystery. One approach to this question is to use mathematical models that
represent the dynamics of biological neurons from a microscopic chemical
viewpoint. You can simulate and analyze spiking phenomena, and investigate the
functions of the biological neural system that the model mimics. These models,
which originated in the Hodgkin-Huxley model [1] and the FitzHugh-Nagumo model
[2, 3], can reproduce various types of neural behavior. Because of their
complexity, however, it is difficult to analyze the neural networks composed
by these neuron models. Therefore, motivated by the desire to understand the
efficient computational architecture underlying biological neural networks,
spiking neural networks (SNNs), which are constructed with neuron models
simplified for computation, have been researched. Among them, the leaky
integrated firing (LIF) neuron model [4] and its variants have been widely
studied due to their mathematical tractability. SNNs are biologically
plausible neural networks and differ significantly from general artificial
neural networks (ANNs) in that they encode information into spikes and
propagate the spikes through synapses. We focus on the dynamic temporal coding
properties specific to neural circuits and SNNs in terms of the encoding
scheme and the learning mechanism of the neural system. In the neural networks
of biological organisms, it was shown experimentally that not only the spike
firing rate but also the individual precise spike timing held information, and
that it was an important factor in achieving more efficient computation [5,
6]. In recent years, gradient methods for end-to-end supervised learning via
direct spiking have been becoming established, but a number of research
challenges remain due to the difficulties specific to SNNs.

One of the fundamental challenges in the supervised learning of SNNs using
gradient methods is the non-differentiability of the discontinuous binary
spike function. Various approaches have been developed and successfully
overcome this difficulty [7, 8, 9]. Most of those used in the recent research
can be divided into two categories depending on a state quantity, through
which the error signals are backpropagated: the rate-based, and the timing-
based method. First, the rate-based method calculates the derivatives between
time-evolving functions such as membrane potentials and spike functions, which
leads to excellent results in a variety of tasks [10, 11]. However, to achieve
the functional derivatives, the rate-based method requires the time
discretization of each function and the approximation of the gradient
computation by replacing discontinuous spike gradients with continuous
surrogate gradients (SG) during error backpropagation. In addition, there are
concerns about how accurately the transfer of information is approximated due
to errors introduced by these techniques of the rate-based method. Note that
the use of SGs, which are continuous real-valued functions, can be interpreted
as approximating the spike rate of each neuron with an activation function,
suggesting that the network is not efficiently dealing with exact spike
timing. Therefore, it is somewhat questionable whether it is appropriate as a
temporal-coding scheme.

In contrast, the timing-based method spearheaded by SpikeProp [12] provides a
solution for computing exact gradients using spike timing. What is more,
Mostafa [13] explicitly derives spike timing gradients for non-leaky IF neuron
models, leading to superior performance. Furthermore, Comsa et al. [14] and
Göltz et al. [15] extended the method to LIF neuron models that also account
for current leakage. The rate-based method backpropagates the gradient from
the postsynaptic current to the presynaptic membrane potential regardless of
the presence of spike firing events, whereas the timing-based method
backpropagates the gradient between layers via the spike timing. Therefore,
the timing-based method can avoid the problem of indifferentiability without
the rate-based coding scheme. This point suggests the potential of the
temporal-coding properties of timing-based methods. In general, timing-based
methods adopt time-to-first-spike (TTFS) as a coding method, which converts
the intensity of input information into absolute spike timing, and also they
impose single-spike restrictions on all neurons [16, 14, 15, 17]. These ideas
are supported by the physiological studies [18, 5] that argue that the first
spike retains more information than the later ones. On the other hand, we
define multi-spike SNNs as the networks in which individual neurons are
allowed to fire multiple times, in contrast to SNNs with single-spike
restrictions. It is clear that the single-spike restriction per neuron limits
the learning efficiency and information processing capacity of SNNs.
Furthermore, since biological neural networks are multi-spike, multi-spike is
essential for investigating the coding properties of SNNs. Some researches
tried on timing-based learning of multi-spike SNNs, but their results do not
indicate that multi-spike SNNs were effective enough to outperform SNNs with
single-spike restrictions [19, 20].

In response to the need to understand the characteristics of multi-spike SNN
models, we propose a new algorithm for timing-based error backpropagation of
multi-spike SNNs. To the best of our knowledge, this is the first study of
supervised learning of SNNs based on spike timing allowing multi-spike. It is
a natural extension of the timing-based method that originated from the
Mostafa [13] in that it adopted the LIF model as the neuron model and did not
involve the intervention of heuristic techniques. The contributions of this
paper are as follows:

  * •

We propose a timing-based backpropagation for multi-spike SNNs with exact
gradients by constructing a learning algorithm via spike timing only.

  * •

Using the MNIST dataset, a fundamental benchmark for image recognition, we
revealed for the first time the fundamental properties of multi-spike
dynamics, the relationship between the leakage time constant and the number of
spikes, and its performance improvement over single-spike restricted timing-
based methods.

  * •

We propose a spike-count loss that can be used in gradient methods based
solely on spike timing, called a dead neuron penalty loss that pulls back
neurons that do not respond to any input pattern into learning.

We hope that these results will provide a solid foundation for research on
gradient learning methods for timing-based SNNs and serve as a good benchmark
for working with more complex networks and tasks.

##  II Methods

We propose a learning algorithm for feedforward SNNs based solely on spike
timing. In this white paper, we employ the spike timing sequence as the output
of each layer. Note that we use spike timing sequences based on two different
representations of ”spike order” depending on the situation: local spike order
and global spike order. The former refers to the order in which spikes are
output from a single neuron. The local spike order-based spike timing of
neuron $i$ is defined as $t_{i}^{(m)}$, which means the spike timing when
neuron $i$ fires for the $m$-th time. The latter, on the other hand, is the
spike order consistently counted for all output spikes from all neurons in a
layer, which is an integration of the local spike order over all neurons in
the layer. The spike timing based on the global spike order is defined as
$\hat{t}^{(m)}$, which is the spike timing of the $m$-th spike in this layer.
When the total number of spikes in a layer is $M$ and the number of output
spikes for each neuron $i\in\\{1,\cdots,I\\}$ is $M_{i}$, these satisfy
$M=\sum_{i=1}^{I}M_{i}$.

###  II-A Neuron Models

We use a leaky integrate-and-fire (LIF) neuron in continuous time. The neuron
model involves reset dynamics, which is induced as an output spike is
generated. Suppose that the local spike timing sequence of neuron
$i\in\\{1,\cdots,I\\}$ in the previous layer is
$\\{t_{i}^{(m)}\\}_{m=1}^{M_{i}}$. The differential equations satisfied by the
postsynaptic current $I_{j}(t)$ and the membrane potential $V_{j}(t)$ of
postsynaptic neuron $j\in\\{1,\cdots,J\\}$ at time $t\in\mathbb{R}_{\geq 0}$
are defined as

$\displaystyle S_{i}(t)$
$\displaystyle=\sum_{m=1}^{M_{i}}\delta(t-t_{i}^{(m)}),$ (1)
$\displaystyle\frac{\mathrm{d}I_{j}(t)}{\mathrm{d}t}$
$\displaystyle=-\frac{1}{\tau_{I}}I_{j}(t)+\beta_{I}\sum_{i=1}^{I}w_{ij}S_{i}(t),$
(2) $\displaystyle\frac{\mathrm{d}V_{j}(t)}{\mathrm{d}t}$
$\displaystyle=-\frac{1}{\tau_{V}}V_{j}(t)+\beta_{V}I_{j}(t),$ (3)
$\displaystyle V_{j}(t)$ $\displaystyle\leftarrow 0\quad,\text{when}\
V_{j}(t)=V_{th},$ (4)

where $S_{i}(t)$ is the spike sequence of the presynaptic neuron $i$, and
$\delta(t)$ denotes Dirac delta function. $w_{ij}\in\mathbb{R}^{I\times J}$ is
a synaptic weight from presynaptic neuron $i$ to postsynaptic neuron $j$.
$\tau_{I}\in\mathbb{R}_{>0}$ and $\tau_{V}\in\mathbb{R}_{>0}$ denote the decay
time constants for the postsynaptic current and the membrane potential,
respectively. $V_{th}$ is the threshold potential, which is fixed at $1.0$.
$\beta_{I}\in\mathbb{R}_{>0}$ and $\beta_{V}\in\mathbb{R}_{>0}$ are the scale
coefficients, but, for simplicity, we fix both $\beta$s to a $1$, hereafter.

Forward pass based on spike timing In this study, a LIF model with decaying
voltage and current with $\tau_{I},\tau_{V}<\infty$ are adopted, and we let
$\tau_{I}\neq\tau_{V}$ consistently. Assume here that neuron $j$ does not
fire, i.e., this neuron satisfies the initial values $I_{j}(0),V_{j}(0)=0$ at
time $t=0$, continues to receive spikes from the previous layer. In this case,
the postsynaptic potential of neuron $j$ at time $t$ and it’s current,
respectively, are obtained as convolution integral representations:

$\displaystyle I_{j}(t)$
$\displaystyle=\int_{0}^{t}\theta(t^{\prime})\mathcal{A}(t^{\prime})\left(\sum_{i=1}^{I}w_{ij}S_{i}(t-t^{\prime})\right)\mathrm{d}t^{\prime},$
(5) $\displaystyle V_{j}(t)$
$\displaystyle=\frac{\tau_{I}\tau_{V}}{\tau_{I}-\tau_{V}}\int_{0}^{t}\theta(t^{\prime})\mathcal{K}(t^{\prime})\left(\sum_{i=1}^{I}w_{ij}S_{i}(t-t^{\prime})\right)\mathrm{d}t^{\prime},$
(6) $\displaystyle\theta(t)$
$\displaystyle=\left\\{\begin{array}[]{ll}1&\text{if}\quad t\geq 0\\\
0&\text{otherwise}\end{array}\right.,$ (9)

where $\theta(t)$ is the step function and the kernel $\mathcal{K}$ and the
function $\mathcal{A},\mathcal{B}$ are defined as

$\displaystyle\mathcal{K}(t)$ $\displaystyle=\mathcal{A}(t)-\mathcal{B}(t),$
(10) $\displaystyle\mathcal{A}(t)$
$\displaystyle=\exp\left(-\frac{t}{\tau_{I}}\right),$ (11)
$\displaystyle\mathcal{B}(t)$
$\displaystyle=\exp\left(-\frac{t}{\tau_{V}}\right).$ (12)

That is, the membrane potential $V_{j}(t)$ is represented as the weighted
convolution integral of this kernel $\mathcal{K}$ with respect to all spikes
in the output of the previous layer.

Then we consider the neuronal dynamics after the postsynaptic neurons fire.
That is, it resets the membrane potential as an output spike is generated.
Note that the synaptic current is not reset at this time. Unlike in the
existing studies of SNNs using timing-based methods [13, 14, 15], where the
single-spike restrictions let the activity of the neuron totally stopped at
its own spike timing, the reset term results in multi-spike in this study.
Thereafter, $\hat{t}^{(m)}$ denotes the presynaptic spike timing at $m$-th in
the global spike order, and the intersynaptic weight from the presynaptic
neuron $i$ which fires at $m$-th in the global spike order to the postsynaptic
neuron $j$ is defined as $w_{ij}^{(m)}$. Assuming that time $t$ satisfies
$t\in[\hat{t}^{(m^{\prime})},\hat{t}^{(m^{\prime}+1)})$ and that a neuron $j$
newly fires during this interval at $t=t_{j}^{(n^{\prime}+1)}$. At any
$t\in[\hat{t}^{(m^{\prime})},\hat{t}^{(m^{\prime}+1)})\land[t_{j}^{(n^{\prime})},t_{j}^{(n^{\prime}+1)})$,
the synaptic current and membrane potential of neuron $j$ can be written
respectively as

$\displaystyle I_{j}(t)$
$\displaystyle=\sum_{m=1}^{m^{\prime}}w_{ij}^{(m)}\mathcal{A}(t-\hat{t}^{(m)}),$
(13)
$\displaystyle\begin{split}V_{j}(t)&=\frac{\tau_{I}\tau_{V}}{\tau_{I}-\tau_{V}}\sum_{m=1}^{m^{\prime}}w_{ij}^{(m)}\left[\mathcal{A}(t-\hat{t}^{(m)})\right.\\\
&-\left.\max\left(\mathcal{A}(t_{j}^{(n^{\prime})}-\hat{t}^{(m)}),\mathcal{B}(t_{j}^{(n^{\prime})}-\hat{t}^{(m)})\right)\mathcal{B}(t-t_{j}^{(n^{\prime})})\right].\end{split}$
(14)

For convenience, we define the leaky factor of the two different time
constants $\tau_{I,V}$ as $p=\tau_{I}/\tau_{V}$. Also, since spike timing only
appears in an exponential form such as $\exp(-t/\tau_{I}),\exp(-t/\tau_{V})$,
we can transform times $t$ to $z$ according to the definition of
$z:=\exp(t/\tau_{I})$, e.g., $\hat{z}^{(m)}:=\exp(\hat{t}^{(m)}/\tau_{I})$ and
$z_{j}^{(n)}:=\exp(t_{j}^{(n)}/\tau_{I})$. Both $I_{j},V_{j}$ at any
$z\in[\hat{z}^{(m^{\prime})},\hat{z}^{(m^{\prime}+1)})\land[z_{j}^{(n^{\prime})},z_{j}^{(n^{\prime}+1)})$
are multiplied by the spike timing sequence using
$\\{\hat{z}^{(m)}\\}_{m=1}^{M}$ as

$\displaystyle I_{j}(z)$
$\displaystyle=\sum_{m=1}^{m^{\prime}}w_{ij}^{(m)}\frac{\hat{z}^{(m)}}{z},$
(15)
$\displaystyle\begin{split}V_{j}(z)&=\frac{\tau_{I}\tau_{V}}{\tau_{I}-\tau_{V}}\\\
&\sum_{m=1}^{m^{\prime}}w_{ij}^{(m)}\left[\frac{\hat{z}^{(m)}}{z}-\frac{\hat{z}^{(m)}(\max\\{z_{j}^{(n^{\prime})},\hat{z}^{(m)}\\})^{p-1}}{(z)^{p}}\right],\end{split}$
(16)
$\displaystyle=\frac{\tilde{\mathcal{A}}_{j}^{(m^{\prime})}}{z}-\frac{\tilde{\mathcal{B}}_{j}^{(m^{\prime},n^{\prime})}}{(z)^{p}},$
(17)

where, for convenience, we use the following definitions:

$\displaystyle\tilde{\mathcal{A}}_{j}^{(m^{\prime})}$
$\displaystyle=\frac{\tau_{I}\tau_{V}}{\tau_{I}-\tau_{V}}\sum_{m=1}^{m^{\prime}}w_{ij}^{(m)}z^{(m)},$
(18) $\displaystyle\tilde{\mathcal{B}}_{j}^{(m^{\prime},n^{\prime})}$
$\displaystyle=\frac{\tau_{I}\tau_{V}}{\tau_{I}-\tau_{V}}\sum_{m=1}^{m^{\prime}}w_{ij}^{(m)}\hat{z}^{(m)}(\max\\{z_{j}^{(n^{\prime})},\hat{z}^{(m)}\\})^{p-1}.$
(19)

For each case, a solution for the spike time $z_{j}^{(n^{\prime}+1)}$, defined
by

$\displaystyle V_{j}(z)=V_{th}$ (20)

has to be found, sequentially. In the whole experiments in the study, let the
leaky factor $p=2$. Therefore, we can solve the next spike timing
$z_{j}^{(n^{\prime}+1)}$ can be solved explicitly and analytically by using
the coefficient parameters
$\tilde{\mathcal{A}}_{j}^{(m^{\prime})},\tilde{\mathcal{B}}_{j}^{(m^{\prime},n^{\prime})}$
based on the historical spike timing sequence of the previous layer
$\\{z^{(m)}\\}_{m=1}^{m^{\prime}}$ and its own previous spike timing
$\hat{z}_{j}^{(n^{\prime})}$ as

$\displaystyle z_{j}^{(n^{\prime}+1)}$
$\displaystyle=\frac{\tilde{\mathcal{A}}_{j}^{(m^{\prime})}-\sqrt{\left(\tilde{\mathcal{A}}_{j}^{(m^{\prime})}\right)^{2}-4V_{th}\tilde{\mathcal{B}}_{j}^{(m^{\prime},n^{\prime})}}}{2V_{th}}.$
(21)

The condition for the existence of $z_{j}^{(n^{\prime}+1)}$ is that each
coefficient parameter
$\tilde{\mathcal{A}}_{j}^{(m^{\prime})},\tilde{\mathcal{B}}_{j}^{(m^{\prime},n^{\prime})}$
and the root content of the equation (21) are non-negative and
$z_{j}^{(n^{\prime}+1)}\in[z^{(m^{\prime})},z^{(m^{\prime}+1)})$. When
actually implementing the algorithm, it can be computed independently for each
$m^{\prime},j$ in accordance with the definition of $z_{j}^{(n^{\prime}+1)}$.
The obtained spike timing sequences $\\{z_{j}^{(n^{\prime})}\\}_{n^{\prime}}$
for each neuron are finally merged into a unified spike timing sequence
$\hat{z}^{(n)}$ in the global spike order, and then subsequent layers can be
calculated in the same way.

Figure 1: Computational graph of a multi-spiking neural network (4-2-3).
Inputs and hidden units are spike units. The SNN has only two neurons in the
hidden layer and, in this case, neuron $1$ and neuron $2$ fire three times and
twice in a given period, respectively. (a) The forward pass is illustrated as
colored arrows: the green solid arrows indicate the propagation of spikes
through the synapse or the reset path, the green dashed arrows indicate
merging the spikes of each neuron in the global spike order. (b)
Backpropagation flow of computational for $\frac{\partial\mathcal{L}}{\partial
z_{1}^{(1)}}$ with orange arrows.

Fig. 1 shows the computational graph where edges start from the input spike
sequence to the loss function and its backward. In the hidden layer, the
computation of the spike timing is done independently for each neuron and
finally merged into a spike timing sequence based on the global spike order
for the entire layer (the dashed arrows). The arrows between spikes within the
same neuron represent the reset paths derived from their own most recent past
spike timing, which are induced by the reset of the membrane potential. When
calculating the hidden layer spike timing sequence from the input spike timing
sequence, the computational path from the input spike received earlier will
continue to affect the calculation of spike timing after that received time
via two types of the path: the postsynaptic currents that remain after resets,
and the latest spike timing of themselves. Note that, considering the
convolution integral representations above, the time range that a neuron is
affected by a previous single spike depends on the leakage time constant. In
the case of a multi-layered SNNs, the spike timing sequence based on the
global spike order of the previous hidden layer can be regarded as the input
spike sequence, and the spike timing sequence of the next layer can be
calculated just as the spike timing is calculated from the input spike
sequence in the figure. The second graph illustrates to computation for
$\frac{\partial\mathcal{L}}{\partial z_{1}^{(1)}}$ by backpropagation from the
bottom $\mathcal{L}$ (the orange arrows). Note that the gradient is propagated
backward without ignoring all the reset paths.

###  II-B Input formats

The input format to the network is latency coding. The latency input consists
of the same number of input neurons as the dimension of the input data, and
each value is mapped to the spike timing of each input neuron in a one-to-one
fashion. Each input neuron fires only once per sample, and its exact spike
timing contains the information. Every input neuron has a one-way weighted,
fully connected, directed path to all neurons in the hidden layer, and spikes
carrying information at the spike timing are input as synaptic currents
through this path to each middle layer.

###  II-C Loss functions

We adopt membrane potential loss as the loss function. The same number of
output neurons as a given number of output dimensions is placed in the final
layer, and the output from one neuron corresponds to one output dimension. The
threshold is set to infinite only for the output neurons, and they never fire.
The value of the membrane potential of each output neuron was observed at the
preset end time $t_{out}$, and this value was taken as the final output of
each neuron. In other words, the output membrane potential
$\textbf{v}_{out}:=\\{v_{o}\\}_{o=1}^{N}$ of this network is defined by the
output layer coefficient parameters
$\mathcal{A}_{o}^{(M)},\mathcal{B}_{o}^{(M,0)}$ as

$\displaystyle
v_{o}=\frac{\mathcal{A}_{o}^{(M)}}{z_{out}}-\frac{\mathcal{B}_{o}^{(M,0)}}{(z_{out})^{2}},$
(22)

where $z_{out}:=\exp(t_{out}/\tau_{I})$. For the multiclass classification
problem, the cross-entropy, which incorporates the correct answer labels into
the softmax of these output values, was used as the error function.

In this paper, we also introduce spike count penalty loss for efficient
learning. By incorporating this into the overall model loss function, the
number of spikes in each neuron can be freely manipulated even in the gradient
method for SNNs based only on spike timing. This penalty loss provides one
solution to the challenge that one of the difficulties in learning timing-
based SNNs is that it is unable to incorporate the spike counts directly into
the model because the timing-based error backpropagation has no state
variables other than the spike timing.

To incorporate the spike count penalty loss, the number of spikes for each
neuron that fires during the forward trial process is counted, and a flag

$\displaystyle\mathcal{P}:=\left\\{\mathcal{P}\in\\{0,\ 1\\}^{J}\ |\
\mathcal{P}_{j}=c(M_{j}),\ \forall j=1,\dots,J\right\\}$ (23)

is the vector for each neuron that meets the specified spike count condition,
where $M_{j}$ is the spike count of neuron $j$ and $J$ is the number of
neurons in the target hidden layer. $c:\mathbb{Z}_{\geq
0}\rightarrow\\{0,1\\}$ is a binary function that defines the spike count
condition. Let $\textbf{v}_{hidden}$ be a $J$-dimensional vector of the
measured values of each output neuron membrane potential at the end time of
the trial $t_{out}$, and then the spike count penalty loss,

$\displaystyle
L_{c}(\textbf{v}_{hidden};\mathcal{P})=\frac{1}{J}\mathcal{P}^{\top}(-\textbf{v}_{hidden}+V_{th}\textbf{1}_{J})$
(24)

is obtained, where $\textbf{1}_{J}$ is an $J$-dimensional 1-vector. By
incorporating this penalty loss, the spike counts of individual neurons can be
conditioned into the SNN timing-based learning algorithm.

In particular, by penalizing dead neurons that do not fire for any input
pattern, it is possible to suppress initial value dependence, which has been
considered difficult in timing-based methods, and to facilitate the
experimental design. In this case, we refer to the loss of spike count penalty
as ’dead neuron loss’. Importantly, by setting the conditional equation to
flag dead neurons, we can avoid the problem of neurons that do not fire
depending on the initial weight variable (since the timing-based method only
handles spike timing, neurons that never fire cannot be included in the
training, and all gradients will be zero, making them untrainable). This
reduces the dependence on the initial weight variables, which is a major
concern with timing-based methods, and makes SNN training design much easier.

The square of the norm of the output multiplied by a small factor is also
added to the overall loss to prevent output values from diverging. Although
learning is possible without this loss term, we added it in our experiments
because we observed that the loss function of the test increased in the middle
of learning, even though the accuracy increased.

In summary, for the obtained network membrane potential output $v_{o}$ and
target class indicator probability $q$, the error function $E$ is

$\displaystyle E(\textbf{v}_{out})$
$\displaystyle=-\sum_{o}q_{o}\log\text{softmax}(v_{o}),$ (25)

and the loss function $\mathfrak{L}(\textbf{v}_{out})$ for the whole model is
given by using the hyperparameters $\lambda,\ \sigma$ as

$\displaystyle\mathfrak{L}(\textbf{v}_{out})$
$\displaystyle=E(\textbf{v}_{out})+\lambda
L_{c}(\textbf{v}_{hidden};\mathcal{P})+\sigma\frac{1}{N_{out}}||\textbf{v}_{out}||_{2}^{2}.$
(26)

###  II-D Memory- and computation-efficient algorithms

Training and inference of SNNs on ordinary Von Neumann type computers are
extremely computation- and memory-intensive due to the nature of the dynamics
in continuous time. Although research on SNNs has been accelerated in recent
years with the rapid development of computers, its computational complexity
still remains an obstacle to research on SNNs. In fact, this limitation has
made it difficult to handle large-scale datasets and models. To realize a
computationally efficient timing-based error backpropagation method that
allows multi-spike, we adopted a computation cycle based on the spike count
independent of the spike timing and introduced the two-stage partitioned
algorithm. First, the spike-count-based computation loop enables parallel
processing of calculations proportional to the input dimension and the number
of hidden neurons, which can be efficiently computed on GPUs and TPUs to
achieve a fast computation rate. Furthermore, though there is a trade-off
between computation time and memory resources, and if the input dimension or
model parameters are large, the burden on memory increases when covering a
large number of spikes, the two-stage cycle algorithm can dramatically reduce
the amount of memory required for computation. The sequence of processes of
the proposed algorithm with two-stage cycle processing is shown in Algorithm
1. In this algorithm, the iterator $j,m^{\prime}$ loops can be computed in
parallel.

Algorithm 1 Forward pass of multi-spiking neural networks

0: $\\{\hat{t}^{(m)}\\}_{m}$: Vector of input spike timing

0: $M$: Number of input spikes

0: $J$: Number of hidden neurons

0: $O$: Number of output neurons

0: $\\{w_{ij}\\}_{i,j},\\{w_{jo}\\}_{j,o}$: Set of weight matrix

0: $N_{1},\ N_{2}$: Maximum number of spikes to break the loop

0: $t_{out}$: Time to observe output membrane potential

0: $\tau_{I}$: Leakage time constant

0: $\\{v_{o}\\}_{o}$: Vector of output membrane potential Initialisation :

1: $\hat{z}^{(m)}\leftarrow\exp(\hat{t}^{(m)}/\tau_{I}),\ \forall
m\in\\{1,\cdots,M\\}$

2: $z_{j}^{(0)}[m^{\prime}]\leftarrow 0,\ \forall j\in\\{1,\cdots,J\\},\forall
m^{\prime}\in\\{1,\cdots,M\\}$

3: $z_{out}\leftarrow\exp(t_{out}/\tau_{I})$ 1st Loop Process :

4: for $n^{\prime}=0$ to $N_{1}-1$ do

5: for $j=1$ to $J$ (Parallel) do

6: for $m^{\prime}=1$ to $M$ (Parallel) do

7:
$\tilde{\mathcal{A}}_{j}^{(m^{\prime})}\leftarrow\tau_{I}\sum_{m=1}^{m^{\prime}}w_{ij}^{(m)}z^{(m)}$

8:
$\tilde{\mathcal{B}}_{j}^{(m^{\prime},n^{\prime})}\leftarrow\tau_{I}\sum_{m=1}^{m^{\prime}}w_{ij}^{(m)}\hat{z}^{(m)}\max\\{z_{j}^{(n^{\prime})},\hat{z}^{(m)}\\}$

9: if
$\tilde{\mathcal{B}}_{j}^{(m^{\prime},n^{\prime})},\tilde{\mathcal{A}}_{j}^{(m^{\prime})},\left(\tilde{\mathcal{A}}_{j}^{(m^{\prime})}\right)^{2}-4\tilde{\mathcal{B}}_{j}^{(m^{\prime},n^{\prime})}V_{th}\geq
0$ then

10:
$z_{j}^{(n^{\prime}+1)}[m^{\prime}]\leftarrow\frac{\tilde{\mathcal{A}}_{j}^{(m^{\prime})}-\sqrt{\left(\tilde{\mathcal{A}}_{j}^{(m^{\prime})}\right)^{2}-4\tilde{\mathcal{B}}_{j}^{(m^{\prime},n^{\prime})}V_{th}}}{2V_{th}}$

11: if $z_{j}^{(n^{\prime}+1)}[m^{\prime}]<\hat{z}^{(m^{\prime})}$ or
$z_{j}^{(n^{\prime}+1)}[m^{\prime}]>\hat{z}^{(m^{\prime}+1)}$ then

12: $z_{j}^{(n^{\prime}+1)}[m^{\prime}]\leftarrow\infty$

13: end if

14: else

15: $z_{j}^{(n^{\prime}+1)}[m^{\prime}]\leftarrow\infty$

16: end if

17: end for

18:
$z_{j}^{(n^{\prime}+1)}\leftarrow\min_{m^{\prime}}z_{j}^{(n^{\prime}+1)}[m^{\prime}]$

19: $\mathcal{A}_{o}\leftarrow\mathcal{A}_{o}+\beta_{I}\beta_{V}\tau_{I}\
w_{jo}\left(z_{j}^{(n^{\prime}+1)}\right)^{2},\ \forall o\in\\{1,\cdots,O\\}$

20: $\mathcal{B}_{o}\leftarrow\mathcal{B}_{o}+\beta_{I}\beta_{V}\tau_{I}\
w_{jo}\ z_{j}^{(n^{\prime}+1)},\ \forall o\in\\{1,\cdots,O\\}$

21: end for

22: end for2nd Loop Process :

23: for $n^{\prime}=N_{1}$ to $N_{2}-1$ do

24: for $j\in\\{j\ |\ z_{j}^{(N_{1})}<\infty\\}$ (Parallel) do

25: the same process as 1st Loop

26: end for

27: end for

28:
$v_{o}\leftarrow\frac{\mathcal{A}_{o}}{z_{out}}-\frac{\mathcal{B}_{o}}{(z_{out})^{2}},\
\forall o$

29: return $\\{v_{o}\\}_{o}$

##  III Results

We use the MNIST benchmark to illustrate the high performance of our method
and the properties of SNNs acquired by multi-spike. Our training model is
characterized by the application of error backpropagation to a feed-forward
SNN with multi-spike, without using surrogate gradient or plasticity. The
weight variables were initialized according to a Gaussian distribution with
expectation and variance of $0.03$ and $0.3$, respectively, for all tasks.
There were $10$ output units since all experiments were solving MNIST
classification problems, and the inputs were composed of the same number of
input neurons and spikes as image dimensionality $784(=28\times 28)$. The
number of neurons in the hidden layer is consistently set to $400$, making a
two-layer network ($784-400-10$ network). Weights were updated using Adam [21]
and trained $100$ epochs with a learning rate of $0.001$. With respect to the
input, the MNIST dataset consists of $70000$ grayscale handwritten digit
images and correct answer labels consisting of any integer from $0-9$. Of
these, $60000$ were used for training and the remaining $10000$ for testing.
The pixel intensity in each dimension was scaled to a value between $[0,1]$
after inverting the black and white of the image and converted to input spike
timing between $t=0.0$ and $t=1.0$. This was split into mini-batches of $100$
each and trained. The output was the membrane potential of the output neuron
at the end time $t_{\text{end}}=1.0$. As a dead neuron loss, we flagged
neurons that only fire on input patterns less than $1/10$ of the mini-batch
size. Hyperparameters were set as $\lambda=0.01,\sigma=0.0001$ after grid-
searched, but their effect on learning performance was small unless taken to
extremes. Bias spiking was not used in any of the experiments. All programs
were implemented on PyTorch [22].

Figure 2: Time evolution of the membrane potential and postsynaptic current of
each neuron when a certain input pattern is given to the learned network. At
the top, a single round symbol indicates a single input spike. The second from
the top shows the membrane potential and postsynaptic current of a neuron
selected in the hidden layer, represented by solid and dotted lines,
respectively. The third figure from the top captures the temporal evolution of
the membrane potential of 20 neurons in the hidden layer. The bottom figure
illustrates the membrane potential of each of the 10 output neurons over time.

Fig. 2 shows the dynamics of neurons in the hidden layer and the output layer
when an input spike sequence was input to a trained SNN with $\tau_{I}=1.0$
under the above parameter settings. The number of spikes for each neuron
varies significantly, and it can be seen that the firing frequency of the same
neuron fluctuates along with the postsynaptic current within a given time
($t_{out}=1.0$ in this case).

First, we investigated how the number of spikes of SNNs can be altered under
multi-spike conditions. In the two-layer feedforward SNNs described above, the
entire computational process is based solely on the spike timing of each
spike. It is worth noting what distribution of the number of spikes per neuron
emerges when multi-spike is allowed.

Figure 3: Distributions of the number of spikes in SNNs (784-400-10) for
various leakage time constants $\tau_{V}=0.2,0.4,6.4$.

Fig. 3 shows histograms of the number of spikes of each hidden neuron, which
is averaged over ten different SNNs trained with different initial weights. As
this figure shows, the distribution depends on the magnitude of the current
and voltage leakage terms, i.e., the leakage time constant. Fig.3 shows that
the smaller $\tau_{I}$ is, the larger the mode of the distribution of the
number of spikes shifts toward. For $\tau_{I}=0.2$, the distribution shifts
unimodal to multimodal and the higher mode of the distribution appears
separated into non-$0$ areas. This indicates a bifurcation of all neurons in
an SNN into two unconsolidated classes, neurons that do not fire at all and
neurons that do fire at least once, which significantly hinders learning.

Figure 4: The relationship between the average number of spikes per neuron and
the leakage constant time $\tau_{I}$. the dashed horizon line denotes $1$,
which the SNN model with single-spike restrictions cannot exceed.

Fig. 4 shows the spike count per hidden neuron in each trained SNN when the
leakage time constant is varied, averaged over all hidden neurons. Fig. 4
indicates that the smaller the time constant of the synaptic leakage term,
i.e., the larger the leakage rate, under conditions in which multi-spike are
allowed, the more average number of neurons fire. By contrast, as the leakage
time constant increases, i.e., as the model asymptotically approaches a non-
leaky neuron model, the average number of spikes per neuron converges to a
certain value $\simeq 1.5$. With the single-spike restrictions, the average
number of spikes per neuron is equal to the density of the entire network, the
proportion of neurons that fire with respect to an input pattern. The figure
with single-spike restrictions indicates that when the leakage time constant
approaches close to zero, almost all neurons will fire. The SNNs with the
single-spike restrictions converge to 0.6 asymptotically (non-leaky models),
but unlike the multi-spike model, the curve is neither smooth nor monotonic.
In terms of sparsity on a per-neuron basis, the multi-spike SNNs are sparser
than the SNNs with single-spike restrictions over a wide range of leakage time
constants.

Figure 5: The relationship between learning performance and the leakage
constant time $\tau_{I}$.

Fig. 5 shows the relationship between leakage time constant and test accuracy
on the MNIST dataset. First, the multi-spike SNNs outperform the SNNs
restricted to single-spike in accuracy under any leakage time constant
condition. Common to both multi-spike and single-spike SNNs is that as the
leakage time constant is very close to zero accuracies deteriorate rapidly as
leakage increases. For both cases, there is a convergence value when the
leakage time constant is increased, converging to $98.23\%$ for multi-spike
and $97.96\%$ for single-spike, respectively. Note that there is a clear
optimum leakage time constant $\tau_{I}=0.8$ (with the test accuracy at
$98.43\pm 0.05\%$) only for the multi-spike case.

TABLE I: Performance of supervised learning for neural networks on MNIST tasks
by neuron type

Network | Coding | BP scheme | Spike type | Leakage | $\tau_{V}$ | $\tau_{I}$ | Test accuracy
---|---|---|---|---|---|---|---
784-800-10 (+ dropout)[23] | ANN |  |  |  |  |  | $98.4\%$
784-800-10[24] | SNN (rate coding) | surrogate-gradient | multi | leaky | $\tau$ | $\tau$ | $98.64\%$
784-800-10[25] | SNN (rate coding) | surrogate-gradient | multi | leaky | $\tau$ | $\tau$ | $98.84\pm 0.02\%$
784-800-10 [13] | SNN (temporal coding) | timing-based | single | non-leaky | $\infty$ | $\tau$ | $97.55\%$
784-340-10 [14] | SNN (temporal coding) | timing-based | single | leaky | $\tau$ | $\tau$ | $97.96\%$
784-500-10 [17] | SNN (temporal coding) | timing-based | single | non-leaky | $\infty$ | $\infty$ | $97.99\pm 0.07\%$
784-350-10 [15] | SNN (temporal coding) | timing-based | single | leaky | $\tau$ | $2\tau$ | $97.2\pm 0.1\%$
256-246-10 [26] | SNN (temporal coding) | surrogate-gradient | multi | leaky | $\tau$ | $\tau$ | $97.5\pm 0.1\%$
784-400-10 [This work] | SNN (temporal coding) | timing-based | single | non-leaky* | $\tau$ | $2\tau$ | $97.96\pm 0.08\%$
784-400-10 [This work] | SNN (temporal coding) | timing-based | single | leaky | $\tau$ | $2\tau$ | $\bf{97.99\pm 0.06\%}$
784-400-10 [This work] | SNN (temporal coding) | timing-based | multi | non-leaky* | $\tau$ | $2\tau$ | $98.23\pm 0.07\%$
784-400-10 [This work] | SNN (temporal coding) | timing-based | multi | leaky | $\tau$ | $2\tau$ | $\bf{98.43\pm 0.05\%}$

*It can be regarded as non-leaky due to extreme conditions: $\tau=3.2$.

A comparison of the results obtained in this experiment with the MNIST test
accuracy of various existing 2-layer feedforward SNNs are summarized in Table
I. The error indications in our study are defined by the standard errors of
the results of 10 experiments performed with different initial weights. The
results show that in supervised learning of simple feedforward SNNs trained
with temporal coding, multi-spike SNNs in the timing-based method achieve the
highest accuracy over all existing methods. In addition, there is already no
significant difference when compared to SNN supervised learning based on rate
coding.

##  IV Discussions

We performed supervised learning on both the existing single-spike-limited
SNNs per neuron and the multi-spike-allowed model in this study and compared
the properties of the error backpropagation method for SNNs based solely on
strict spike timing. In particular, the present study shows that omitting the
single-spike restriction and allowing multiple firings of a neuron improves
the performance of the standard MNIST benchmark under all leakage time
constant conditions. A notable feature of our method is that it uses the exact
spike timing in the backpropagation computation process, thus achieving ideal
temporal coding mediated by spikes on multi-spike SNNs in terms of gradient
computation. In related works on multi-spike SNNs, although types of SG
method, applied an appropriate regularization and achieved sparse, near-
temporal-coding learning [26, 27]. However, since these learning methods using
gradient approximation assume continuous spiking behaviors during
backpropagation, it is open to question whether they are strictly temporal-
coding learning from the viewpoint of gradient propagation. To the best of our
knowledge, our proposed algorithm, for the first time, achieved the supervised
learning of multi-spike SNNs with backpropagation of exact gradients based on
spike timing.

Our proposed model has advantages in precisely dealing with the reset
dynamics. The resetting behavior at the moment of firing is one of the factors
that significantly influence the information processing mechanism. In general,
studies addressing supervised learning of multi-spike SNNs have introduced
double-reset neuron models, which reset the membrane potential and also the
corresponding postsynaptic current at the same time due to its simplicity [10,
20]. However, there are concerns that the adoption of double-reset means that
any information input up to each spike timing is lost on each firing,
completely breaking the dependency between previous and subsequent spikes
within the same neuron, thus lacking the benefits of multi-spike. By contrast,
we reset only the membrane potential that has reached the threshold to the
baseline, $0$, while the postsynaptic current is not reset and retains its
pre-firing values. In this way, the past information received by the neuron
itself is propagated into the future via both its own past spike timing and
the postsynaptic current path. There is a non-negligible difference in the
network structure whether double-reset is employed or not, and we employed the
latter in this study.

Heuristic approximations and abbreviations in the error backpropagation
calculation process for SNNs are also major factors that compose the network
learning paradigm. In general, learning methods using surrogate gradients
ignore the error backpropagation paths led by resets of membrane potential,
what is called reset path, in multi-spike [28, 10, 20]. Backpropagating
through surrogate gradient results in the fact that the reset information is
propagated recursively by a number of times equal to the number of time steps,
even though the actual reset is instantaneous at the spike timing. Even though
the amount of the approximation error induced by a single surrogate derivative
is tolerable, the accumulated error can be huge because the number of time
steps is usually by far larger than the number of spikes [28, 20]. To avoid
this approximation error diverging is the main reason why they ignore the
reset paths. On the other hand, the timing-based method used in this study
does not use any approximation in the whole calculation process of
backpropagation, i.e., it realizes an exact gradient calculation. It should be
noted, however, that the timing-based method cannot capture the occurrence and
disappearance of sudden spikes caused by weight updates and the accompanying
changes in the number of spikes. The sudden spike changes pointed out here are
not spike changes that increase or decrease with continuous changes in
membrane potential in the time direction around the end time $t_{out}$, but
spike changes that occur with discontinuous membrane potential changes in the
middle of a given period $(0,t_{out})$. Although this is a challenge inherent
to the timing-based method as pointed out by SpikeProp [12], it was
experimentally shown that learning works well both for single-spike SNNs and
also for the multi-spike SNNs in this study. The impact of this sudden spike
variation is one of the issues that should be considered in the future.

Through this study, it is clear that the properties of the overall network
spike count are altered depending on the time constants of the membrane
potential and synaptic current, i.e., the synaptic leakage term. The results
basically show that the smaller the leakage time constant (i.e., the more
leakage), the higher the overall network spike count, which leads to higher
learning performance. This change in performance with increasing voltage
leakage is a phenomenon not seen in existing SNNs restricted to single-spike.
In the case of multi-spike SNNs trained through surrogate gradients, through
experiments of rate-coding with feedforward networks, [29] showed that the
leakage of membrane potentials and postsynaptic currents ensures the
robustness of the network, leading to high learning performance. In the
present timing-based method, a similar result may be derived from the
correlation between the increase in the number of spikes and the increase in
performance. In other words, it is possible that the increased rate-coding-
like quality induced by the increase in the number of spikes may contribute to
the superior generalization effect. In contrast, single-spike SNNs are clipped
at a maximum of $1$ in the number of spikes per neuron, so the amount of
leakage of the membrane potential should not affect their performance.
Considering that time constants smaller than the time scale of the experiment
($1\text{ms}$ in this case) would conversely cause a sharp decline in
performance, therefore, there are optimal values for the time constants of
membrane potential and synaptic current for multi-spike SNNs. In particular,
the amount of optimized leakage time constant can be characterized by the
magnitude of the task itself imposed on the SNNs and the local network
structure on which the time scale of the experiment depends. Optimizing the
learning capacity of the entire network through leakage cannot be discussed
only in terms of machine learning aspects. It has been studied and reported
from a physiological perspective that the leakage time constant of a neuron
plays an important role as an indicator to characterize that brain region and
that the diversity of time constants is key to the task-solving ability of
that neural networks [30, 16]. This is in line with the relationship derived
in this study between the optimal leakage time constant and the timescale of
the tasks.

##  V Conclusion

In this study, we proposed a novel efficient backpropagation algorithm for
multi-spike SNN models and pursued the goal of achieving high algorithmic
performance and clarifying its learning characteristics compared to single-
spike SNN models. Moreover, we proposed novel learning techniques such as a
two-stage algorithm based on spike count and the addition of a dead neuron
penalty. The performance of our algorithm on the MNIST dataset outperformed
the state-of-the-art SNNs based on temporal coding and was confirmed to be on
par with the conventional SNNs and ANNs based on rate coding. Focusing on the
leakage time constant of the postsynaptic current and the membrane potential,
we investigated how this affected the characteristics of spike counts of SNNs
with multi-spike neurons. We also showed the existence of an optimal membrane
potential leakage rate that achieved maximum performance, around the time
scale of the task. Moreover, the phenomenon was not observed in training SNNs
with single-spike restrictions. Our multi-spike networks suggested the
possibility of supervised learning based on the spike timing of SNNs with a
biologically plausible recurrent topological structure. In the future, we
would like to pursue timing-based learning of more complex and large-scale
tasks and recurrent networks. Finally, we would like to add that our SNN
model, which learns with information encoded in the spike timing, can also
contribute to the construction of very energy-efficient systems including the
learning process [15, 17].

## Acknowledgments

This work was partially supported by SECOM Science and Technology Foundation,
JST PRESTO Grant Number JPMJPR22C5, the NEC Corporation, Moonshot R&D Grant
Number JPMJMS2021, AMED under Grant Number JP21dm0307009, Institute of AI and
Beyond of UTokyo, the International Research Center for Neurointelligence
(WPI-IRCN) at The University of Tokyo Institutes for Advanced Study (UTIAS),
JSPS KAKENHI Grant Number JP20H05921. Computational resource of AI Bridging
Cloud Infrastructure (ABCI) provided by National Institute of Advanced
Industrial Science and Technology (AIST) was used.

## References

  * [1] A. L. Hodgkin and A. F. Huxley, “A quantitative description of membrane current and its application to conduction and excitation in nerve.” _The Journal of physiology_ , vol. 117, no. 4, pp. 500–544, aug 1952.
  * [2] R. Fitzhugh, “Impulses and Physiological States in Theoretical Models of Nerve Membrane.” _Biophysical journal_ , vol. 1, no. 6, pp. 445–466, jul 1961.
  * [3] J. Nagumo, S. Arimoto, and S. Yoshizawa, “An active pulse transmission line simulating nerve axon,” _Proceedings of the IRE_ , vol. 50, no. 10, pp. 2061–2070, Oct 1962.
  * [4] W. Gerstner and W. M. Kistler, _Spiking Neuron Models: Single Neurons, Populations, Plasticity_. Cambridge University Press, 2002.
  * [5] T. Gollisch and M. Meister, “Rapid neural coding in the retina with relative spike latencies,” _Science_ , vol. 319, no. 5866, pp. 1108–1111, 2008.
  * [6] G. Portelli, J. M. Barrett, G. Hilgen, T. Masquelier, A. Maccione, S. Di Marco, L. Berdondini, P. Kornprobst, and E. Sernagor, “Rank order coding: A retinal information decoding strategy revealed by large-scale multielectrode array retinal recordings,” _eNeuro_ , vol. 3, no. 3, pp. 844–853, 2016.
  * [7] E. O. Neftci, H. Mostafa, and F. Zenke, “Surrogate Gradient Learning in Spiking Neural Networks: Bringing the Power of Gradient-based optimization to spiking neural networks,” _IEEE Signal Processing Magazine_ , vol. 36, no. 6, pp. 51–63, 2019.
  * [8] J. D. Nunes, M. Carvalho, D. Carneiro, and J. S. Cardoso, “Spiking Neural Networks: A Survey,” _IEEE Access_ , vol. 10, pp. 60 738–60 764, 2022\.
  * [9] D. Huh and T. J. Sejnowski, “Gradient descent for spiking neural networks,” _Advances in Neural Information Processing Systems_ , vol. 2018-Decem, no. NeurIPS, pp. 1433–1443, 2018.
  * [10] G. Bellec, F. Scherr, A. Subramoney, E. Hajek, D. Salaj, R. Legenstein, and W. Maass, “A solution to the learning dilemma for recurrent networks of spiking neurons,” _Nature Communications_ , vol. 11, no. 1, dec 2020.
  * [11] C. Lee, S. S. Sarwar, P. Panda, G. Srinivasan, and K. Roy, “Enabling Spike-Based Backpropagation for Training Deep Neural Network Architectures,” _Frontiers in Neuroscience_ , vol. 14, no. February, pp. 1–22, 2020.
  * [12] S. M. Bohte, H. La Poutré, and J. N. Kok, “Error-Backpropagation in Temporally Encoded Networks of Spiking Neurons,” _Neurocomputing_ , vol. 48, pp. 17–37, 2000. [Online]. Available: http://ftp.cwi.nl/CWIreports/SEN/SEN-R0037.pdf
  * [13] H. Mostafa, “Supervised learning based on temporal coding in spiking neural networks,” _IEEE Transactions on Neural Networks and Learning Systems_ , vol. 29, no. 7, pp. 3227–3235, 2018.
  * [14] I. M. Comsa, T. Fischbacher, K. Potempa, A. Gesmundo, L. Versari, and J. Alakuijala, “Temporal Coding in Spiking Neural Networks with Alpha Synaptic Function,” _ICASSP, IEEE International Conference on Acoustics, Speech and Signal Processing - Proceedings_ , vol. 2020-May, pp. 8529–8533, 2020.
  * [15] J. Göltz, L. Kriener, A. Baumbach, S. Billaudelle, O. Breitwieser, B. Cramer, D. Dold, A. F. Kungl, W. Senn, J. Schemmel, K. Meier, and M. A. Petrovici, “Fast and energy-efficient neuromorphic deep learning with first-spike times,” _Nature Machine Intelligence_ , vol. 3, no. 9, pp. 823–835, 2021. [Online]. Available: http://dx.doi.org/10.1038/s42256-021-00388-x
  * [16] A. M. Manea, A. Zilverstand, K. Ugurbil, S. R. Heilbronner, and J. Zimmermann, “Intrinsic timescales as an organizational principle of neural processing across the whole rhesus macaque brain,” _eLife_ , vol. 11, pp. 1–20, 2022\.
  * [17] Y. Sakemi, K. Morino, T. Morie, and K. Aihara, “A Supervised Learning Algorithm for Multilayer Spiking Neural Networks Based on Temporal Coding Toward Energy-Efficient VLSI Processor Design,” _IEEE Transactions on Neural Networks and Learning Systems_ , 2021.
  * [18] E. T. Rolls, L. Franco, N. C. Aggelopoulos, and J. M. Jerez, “Information in the first spike, the order of spikes, and the number of spikes provided by neurons in the inferior temporal visual cortex,” _Vision Research_ , vol. 46, no. 25, pp. 4193–4205, 2006.
  * [19] Y. Xu, X. Zeng, L. Han, and J. Yang, “A supervised multi-spike learning algorithm based on gradient descent for spiking neural networks,” _Neural Networks_ , vol. 43, pp. 99–113, 2013. [Online]. Available: www.elsevier.com/locate/neunet
  * [20] J. Kim, K. Kim, and J. J. Kim, “Unifying activation- And timing-based learning rules for spiking neural networks,” in _Advances in Neural Information Processing Systems_ , vol. 2020-Decem, 2020.
  * [21] D. P. Kingma and J. L. Ba, “Adam: A method for stochastic optimization,” _3rd International Conference on Learning Representations, ICLR 2015 - Conference Track Proceedings_ , pp. 1–15, 2015.
  * [22] A. Paszke, S. Gross, F. Massa, A. Lerer, J. Bradbury, G. Chanan, T. Killeen, Z. Lin, N. Gimelshein, L. Antiga, A. Desmaison, A. Köpf, E. Yang, Z. DeVito, M. Raison, A. Tejani, S. Chilamkurthy, B. Steiner, L. Fang, J. Bai, and S. Chintala, “PyTorch: An imperative style, high-performance deep learning library,” _Advances in Neural Information Processing Systems_ , vol. 32, no. NeurIPS, 2019.
  * [23] N. Srivastava, G. Hinton, A. Krizhevsky, I. Sutskever, and R. Salakhutdinov, “Dropout: A simple way to prevent neural networks from overfitting,” _Journal of Machine Learning Research_ , vol. 15, pp. 1929–1958, 2014.
  * [24] J. H. Lee, T. Delbruck, and M. Pfeiffer, “Training deep spiking neural networks using backpropagation,” _Frontiers in Neuroscience_ , vol. 10, no. NOV, 2016.
  * [25] Y. Jin, W. Zhang, and P. Li, “Hybrid macro/micro level backpropagation for training deep spiking neural networks,” _Advances in Neural Information Processing Systems_ , vol. 2018-December, no. 1, pp. 7005–7015, 2018\.
  * [26] B. Cramer, S. Billaudelle, S. Kanya, A. Leibfried, A. Grübl, V. Karasenko, C. Pehle, K. Schreiber, Y. Stradmann, J. Weis, J. Schemmel, and F. Zenke, “Surrogate gradients for analog neuromorphic computing,” _Proceedings of the National Academy of Sciences of the United States of America_ , vol. 119, no. 4, 2022.
  * [27] Y. Yan, H. Chu, Y. Jin, Y. Huan, Z. Zou, and L. Zheng, “Backpropagation With Sparsity Regularization for Spiking Neural Network Learning,” _Frontiers in Neuroscience_ , vol. 16, no. April, 2022.
  * [28] F. Zenke and T. P. Vogels, “The Remarkable Robustness of Surrogate Gradient Learning for Instilling Complex Function in Spiking Neural Networks,” _Neural computation_ , vol. 33, no. 4, pp. 899–925, 2021.
  * [29] S. S. Chowdhury, C. Lee, and K. Roy, “Towards understanding the effect of leak in Spiking Neural Networks,” _Neurocomputing_ , vol. 464, pp. 83–94, 2021. [Online]. Available: https://doi.org/10.1016/j.neucom.2021.07.091
  * [30] R. Kim and T. J. Sejnowski, “Strong inhibitory signaling underlies stable temporal dynamics and working memory in spiking neural networks,” _Nature Neuroscience_ , vol. 24, no. 1, pp. 129–139, 2021. [Online]. Available: http://dx.doi.org/10.1038/s41593-020-00753-w