# Deep Learning in Spiking Phasor Neural Networks

Connor Bybee
Computational Biology & Redwood Center
UC Berkeley
Berkeley, CA 94720
<EMAIL_ADDRESS>
&E. Paxon Frady
Redwood Center & Intel
UC Berkeley
Berkeley, CA 94720
<EMAIL_ADDRESS>
Friedrich T. Sommer
Redwood Center & Intel
UC Berkeley
Berkeley, CA 94720
<EMAIL_ADDRESS>

###### Abstract

Spiking Neural Networks (SNNs) have attracted the attention of the deep
learning community for use in low-latency, low-power neuromorphic hardware, as
well as models for understanding neuroscience. In this paper, we introduce
Spiking Phasor Neural Networks (SPNNs). SPNNs are based on complex-valued Deep
Neural Networks (DNNs), representing phases by spike times. Our model computes
robustly employing a spike timing code and gradients can be formed using the
complex domain. We train SPNNs on CIFAR-10, and demonstrate that the
performance exceeds that of other timing coded SNNs, approaching results with
comparable real-valued DNNs.

††Work presented at Intel’s Neuromorphic Community Fall 2019 workshop in Graz,
Austria and the UC Berkeley Center for Computational Biology Retreat 2019

##  1 Introduction

The highest classification accuracy on image recognition tasks has been
achieved by real-valued deep neural networks (DNN) with floating point
precision. But recent effort has attempted to reformulate deep learning into
formats more suitable for hardware acceleration and low-energy edge computing.
This work aims to demonstrate the benefits of using spike timing codes in a
way that can advantage deep learning acceleration according to specific
metrics, especially energy efficiency and latency.

Implementation of DNN algorithms on next-generation computing platforms such
as neuromorphic computer architectures and networks of coupled-oscillators is
promising, and may prove to be more energy-efficient and faster than current
hardware based on the Von Neumann architecture [26, 18, 3]. When mapping DNN
neurons to spiking neural networks (SNN), the machine learning and
neuroscience community have largely focused on the use of rate codes [8, 14,
22, 20]. Spiking networks are typically challenged by a non-differentiable
objective function, which prohibits training with backpropagation. Successful
methods include training ANNs with the purpose of converting them to SNNs,
using spiking versions of backpropagation compatible with spike timing-
dependent-plasticity, and use of other learning rules [14, 17]. This paper
presents a method for using spike events to represent complex-valued neural
states, which can easily be differentiated in the complex domain.

In a rate code, a stimulus feature is represented by the firing rate of a
neuron measured by the count of spikes in an extended time window. Conversely,
in a spike timing code a feature is represented by the precise timing of a
single spike. Potential benefits of spike timing codes include (i) efficiency
in terms of the number of spikes per computation, and (ii) computation speed,
since the duration of a step of neural update with a rate code cannot be
shorter than the integration window required for estimating the lowest rate.
On the other hand, spike timing codes are often regarded as brittle and not
suited for implementing robust computations [15].

Recently, Thresholding Phasor Associative Memories (TPAM) were introduced [4]
that illustrate how a complex-valued neural state can be directly mapped to a
spike timing code. Essentially, relative to an ongoing oscillation, the timing
of the spike can be used to indicate the phase of a complex number. This was
shown to provide robust spike timing codes in the context of attractor neural
networks. Here, we extend this timing code to feed forward deep networks, and
show how backpropagation based learning can be applied in the complex domain
and robustly mapped to SNNs. Motivated by that aim, this paper compares and
demonstrates deep learning in spiking phasor networks on benchmarks MNIST [13]
and CIFAR-10 [12].

##  2 Background

###  2.1 Spiking backpropagation

Supervised learning of spike timing codes was explored in [17]. In [17],
neurons are trained to integrate input and fire at a precise time upon
crossing a threshold. Training is performed by first determining the set of
input spikes leading to an output spike, called the causal set. Parameters
updates are then calculated by backpropagating through the causal set. Our
work presents a method to train a spike timing coded SNN without the need to
calculate a casual set.

A method for converting real-valued ANNs to timing coded SNNs was introduced
in [23]. It demonstrates on a simple machine learning task (MNIST) that spike
timing codes are more efficient than rate codes and real-valued DNNs when
considering the number of operations needed to perform a computation. Our
paper extends spike timing coded SNNs to CIFAR-10.

###  2.2 Neuroscience

In Theoretical Neuroscience, it is an open question how a network of spiking
biological neurons can precisely coordinate spikes over time and space.
Izhikevich proposed a spiking network that creates “polychronous” spiking
patterns, as a model for describing neural dynamics in brain circuits [9].
Polychronization means that neurons exhibit time-locked spiking patterns, but
with arbitrary spike timing and not just perfectly synchronous spiking. The
spiking patterns we propose are also polychronous patterns, but these are
generated by different mechanisms.

Precise spike timing patterns are observed in multiple brain regions and cell
types. The electrical properties of the neuron membrane give rise to sustained
intrinsic spike oscillations, subthreshold oscillations, and membrane
resonance [10]. Neuron morphology and localized genetic expression patterns of
membrane ion-channels can create compartments with specific electrical
properties [2, 16]. Rate and timing codes have been explored as an method of
mapping ANNs to models of biological plausible neural networks [25, 24, 17,
20]. Inspired by these observations, this work uses such oscillatory
mechanisms to construct and simulate a SNNs with spike timing codes.

##  3 Methods

###  3.1 Spiking Phasor Neural Networks

Complex-valued neural networks have been explored in several contexts. Noest
[19] introduced phasor neural networks (PNN) with applications as an
associative memory and suggested use with backpropagation. Complex domain
backpropagation was introduced in [5] and [7]. Recent developments have been
made with deep neural networks [28]. We provide a concrete example of a feed
forward PNN trained with backpropagation to implement spike timing codes. This
paper refers to PNNs with the TPAM activation function [4] as Spiking PNNs
(SPNNs).

Figure 1: An illustration of a feed forward PNN with normalizing activation
function or Spiking PNN (SPNN).

In PNNs, neurons are represented by complex-valued phasors
$s_{i}\in\mathbb{C}$. Where
$s_{i}=r_{i}e^{\mathrm{i}\theta_{i}}=r_{i}\cos{\theta_{i}}+\mathrm{i}r_{i}\sin{\theta_{i}}=u_{i}+\mathrm{i}v_{i}$.
Here, $\mathrm{i}=\sqrt{-1}$. Phasors are restricted to the unit circle with
the use of a non-linear activation function that preserves phase information.
$\|s_{i}\|_{2}=1,\forall i\in\\{1,..,N_{s}\\}$

The network units typically consist of a set of input, output, and hidden
variables,
$\mathbf{s}=\\{\mathbf{x},\mathbf{y},\mathbf{h}\\}\in\mathbb{C}^{N_{s}}$. The
synaptic weights and neuron biases, $W_{ij}\in\mathbb{C}$ and
$b_{i}\in\mathbb{C}$, act as attenuating or amplifying phase shifters.

We present a two-layer network with layer units
$\mathbf{h}^{(l)}\in\mathbb{C}^{N_{h}},l\in\\{1,2\\}$. Input units can be
real-valued but are transformed to unit length complex-valued vectors whose
phase is a function of the real-valued input.
$x_{i}=e^{\mathrm{i}\theta_{i}}$, where $\theta_{i}=g_{x}(\tilde{x}_{i})$,
$\tilde{x}_{i}\in\mathbb{R}$. Where $g_{x}:\mathbb{R}\rightarrow[0,\pi]$. The
output of the network is mapped to phases similar to binary phase shift
keying. $y_{i}=e^{\mathrm{i}\theta_{i}}$, where
$\theta_{i}=g_{y}(\tilde{y}_{i})$, $\tilde{y}_{i}\in\\{0,1\\}$ and
$g_{y}:\\{0,1\\}\rightarrow\\{0,\pi\\}$.

Inference proceeds by propagating activity from the input layer to the output
layer:

$\mathbf{h}^{(l)}=f(\mathbf{W}^{(l)}\mathbf{h}^{(l-1)}+\mathbf{b}^{(l)},\Theta)$
(1)

The vector $\mathbf{h}^{(l)}$ is the activity of hidden layer $l$, and for
computing the activity in the first hidden layer, equation (1) is applied to
the input vector, i.e., $\mathbf{h}^{(0)}=\mathbf{x}$. The element-wise
thresholding and normalizing function is given by:

$\displaystyle
f(z_{i},\Theta)=\begin{cases}\displaystyle\frac{z_{i}}{|z_{i}|}&\text{if
}|z_{i}|-\Theta>0\\\ 0&\text{otherwise}\end{cases}$ (2)

where $\Theta$ is a trainable thresholding parameter.

The objective is to minimize the phase difference between the network output
and targets. In supervised learning tasks, often the softmax with cross-
entropy loss is used with the target variables are encoded as one-hot binary
vectors. Here, we encode class labels onto the unit circle. Positive class
labels are encoded to lead negative class labels in phase. A classification is
successful when the phase of the positive class leads the phase of all other
outputs. The class targets are binary variables encoded onto the unit circle,
similar to binary phase shift keying. Positive and negative classes are out of
phase by $180\degree$ or $\pi$ radians. The loss function is given by:

$L=\frac{1}{2}\|\mathbf{y}-\hat{\mathbf{y}}\|^{2}_{2}=N_{y}-\sum_{i=1}^{N_{y}}\cos(\theta_{i}-\hat{\theta}_{i})$
(3)

Where $\theta_{i}$ are the target phase angles and $\hat{\theta}_{i}$ are the
estimated phase angles. Minimizing the objective corresponds to phase aligning
the output and targets:

$\frac{\partial
L}{\partial\hat{\theta}_{i}}=\sin(\theta_{i}-\hat{\theta}_{i})$ (4)

As we will show in the results, PNNs with the TPAM activation function can be
trained with the backpropagation algorithm to minimize equation 3.

###  3.2 Complex domain backpropagation

In [5], complex domain backpropagation was derived for the activation function
below.

$\displaystyle y$ $\displaystyle=f(z)$
$\displaystyle=\frac{z}{c+\frac{1}{r}|z|}$ $\displaystyle=u+\mathrm{i}v$

Where, $z=a+\mathrm{i}b\in\mathbb{C}$, $\\{a,b\\}\in\mathbb{R}$ and
$|z|=\sqrt{a^{2}+b^{2}}$. The partial derivatives are

$\displaystyle\frac{\partial u}{\partial
a}=\begin{cases}\frac{r(b^{2}+cr|z|)}{|z|(cr+|z|)^{2}}&\text{if }|z|>0\\\
\frac{1}{c}&\text{if }|z|=0\end{cases}\;\;\;\;\;\;\frac{\partial u}{\partial
b}=\begin{cases}-\frac{rab}{|z|(cr+|z|)^{2}}&\text{if }|z|>0\\\ 0&\text{if
}|z|=0\end{cases}$ $\displaystyle\frac{\partial v}{\partial
a}=\begin{cases}-\frac{rab}{|z|(cr+|z|)^{2}}&\text{if }|z|>0\\\ 0&\text{if
}|z|=0\end{cases}\;\;\;\;\frac{\partial v}{\partial
b}=\begin{cases}\frac{r(a^{2}+cr|z|)}{|z|(cr+|z|)^{2}}&\text{if }|z|>0\\\
\frac{1}{c}&\text{if }|z|=0\end{cases}$

Considering the case where $\Theta=0$ and $|z|>0$, the SPNN activation
function is found by setting $r=1$ and in $\lim_{c\to 0}$. The partial
derivative can be written in the complex domain by adding together the
individual terms.

$\displaystyle(\frac{\partial u}{\partial a}+\frac{\partial v}{\partial
a})+\mathrm{i}(\frac{\partial u}{\partial b}+\frac{\partial v}{\partial
b})=\frac{(b^{2}-ab)+i(a^{2}-ab)}{(a^{2}+b^{2})^{3/2}}$

###  3.3 Spiking neuron model

The continuous-valued phasor network can be mapped to discrete spiking events
in time. For a given choice $T$ for the length of a cycle, the each phase
angle, $\theta_{i}\in[0,2\pi]$ can be mapped to a time
$t_{i}=\theta_{i}*T/2\pi$.

We first consider a network where all layers are phase-locked to the input
layer. The real-valued inputs are mapped to unique phases on the unit circle.
This creates an ambiguity between the largest and smallest values since $2\pi$
is equivalent to $0$. Therefore, we limited the values to $[0,\pi]$. We chose
this option since it is simple and is sufficient for demonstration purposes.
Each input phasor can also be assigned a fixed, random phase shift.

In the first cycle, the input units spike and each synapse is activated with
magnitude $|W_{ij}|$ after a delay of $\theta_{ij}*T/2\pi$. Each neuron
integrates the information from the current cycle and spikes after the
appropriate delay during the subsequent cycle. Therefore, each cycle
propagates signals to adjacent layers, see figure 2.

###  3.4 Circuit Model

We simulate a circuit model to demonstrate that SPNNs can be implemented by
simple electrical components and are robust to noise induced by numerical
simulation and changing input. A change in soma membrane potential is driven
by a leakage current and combined input from dendrites:

$\displaystyle{\frac{{\mathrm{d}}V_{m}}{{\mathrm{d}}t}}$
$\displaystyle=\frac{g_{l}(V_{l}-V_{m})+g_{c}(V_{d}-V_{m}-\overline{V}_{d})}{C_{m}}$
$\displaystyle{\frac{{\mathrm{d}}\overline{V}_{d}}{{\mathrm{d}}t}}$
$\displaystyle=\frac{(V_{d}-\overline{V}_{d})}{\tau_{d}}$

$V_{m}$ is the membrane potential, $g_{l}$ is the conductance of the leak
channel, $V_{l}$ is the leak reversal potential, $g_{c}$ is the conductance of
between the soma and dendrite, $V_{d}$ is the dendrite potential,
$\overline{V}_{d}$ is the average dendrite potential, $\tau_{d}$ is average
dendrite potential time constant, and $C_{m}$ is the membrane capacitance. To
correct for bias around the threshold voltage, the average dendrite potential
is subtracted from the instantaneous dendrite membrane potential. Dendrites
integrate current from synapses which generate membrane potential oscillations
after a spike arrives:

$\displaystyle V_{d}$ $\displaystyle=\sum_{i}w_{i}V_{s,i}$
$\displaystyle\frac{\mathrm{d}V_{s,i}}{\mathrm{d}t}$
$\displaystyle=\frac{-W_{s,i}}{C_{m}}$
$\displaystyle\frac{\mathrm{d}W_{s,i}}{\mathrm{d}t}$
$\displaystyle=\frac{V_{s,i}}{L}-\frac{W_{s,i}}{\tau_{s}}$

$V_{s,i}$ is the potential at synapse $i$, $w_{i}$ is the synaptic weight.
$W_{s,i}$ are activation parameters which act to generate membrane potential
oscillations at a specific frequency, $L$ is a constant used to tune the
membrane potential frequency, and $\tau_{s}$ is the synaptic time constant
which controls oscillation damping. A spike resets the synapse parameters such
that $V_{s,i}=0$ and $W_{s,i}=W_{spike}$. The neuron spikes when the soma
membrane potential cross a threshold. The refractory period is maintained
until the neuron soma membrane potential becomes negative.

##  4 Results

Figure 2: The first three cycles of inference for a feed forward SPNN trained
on MNIST. The network is phased locked to the input signal. Each cycle allows
information to propagate to the next layer. Each line in the vertical axis is
a neuron. The input, hidden, and output layers are colored red, blue, and
green, respectively. The output neurons are accurately separated in phase. The
input pixel amplitudes are mapped to phase angels inversely proportional to
the magnitude. 2: Similar input amplitudes are phase locked. 2: For
demonstration, each input can have a fixed phase shift while still maintaining
performance.

For visualization purposes, we train a fully connected SPNN on MNIST,
784-512-512-10. We unroll the network activity in time to demonstrate how
signals propagate (Fig. 2 \- left panel). The spike times represent a stable
limit-cycle. For demonstration purposes, we add a fixed, random phase shift to
each component of the input while maintaining correct classification (Fig. 2
\- right panel). A system that has yet to reach the limit-cycle will have a
settling time (transient settling dynamics visible later in Fig. 5).

Convolutional SPNNs are trained on MNIST and CIFAR-10. The convolutional
networks consist of two layers of convolution without padding followed by
three fully connected layers,
conv(6,3x3)$\to$conv(16,3x3)$\to$FC128$\to$FC128$\to$FC10. No normalization or
pooling is used. We train using Pytorch [21]. The Adam optimizer [11] is used
with default parameters and a learning rate of 0.001. We train small networks
compared to state-of-the-ar t for simulation purposes. Future work will
investigate the performance of SPNNs on more difficult tasks and with larger
networks.

Training and test error rates with the convolutional SPNNs were measured for
classifying MNIST data (Fig. 4). Our results were obtained with 32-bit
floating point precision, and therefore should be regarded as an upper-bound
of the performance SPNN hardware implementations are expected to achieve using
parameters from a trained model. The classification precision our model
achieved surpasses the best results previously reported for SNNs with spike
timing codes [23]. Further, SPNNs approach the performance level of comparable
DNNs, our model comes within $0.02\%$ compared to the error rate of the DNN
used for conversion to a SNN. Last, our network outperforms previous methods
for directly training SNNs with spike timing codes [17].

Figure 3: Results for a SPNN on MNIST. The lowest test error rate previously
reported for SNNs with spike timing codes and the ANN used for conversion to
the SNN is $1.43\%$ and $1.04\%$, respectively [23]. The lowest test error
rate previously reported for supervised training of spike timing codes is
$3.03\%$ [17]. SPNNs achieve a test error rate of $1.06\%$.

When trained on CIFAR-10. To our knowledge, there is no previous reports of
performing supervised learning and inference with spike timing codes on
CIFAR-10. However, due to the rather small network size we used, our results
on training and test error rates do not reach the state-of-the-art (Fig. 4).
It is promising that SPNN approach the performance of real-valued DNNs on
CIFAR despite minimizing the mean-squared error opposed to the softmax with
cross-entropy loss.

Figure 4: Training and test accuracy on CIFAR-10. The minimum test error rates
for the relu network and the SPNN network are $30.73\%$ and $34.09\%$,
respectively.

We performed simulations of convolutional SPNNs implemented with Brian2 [6],
an open-source simulator for spiking neural networks. The simulations of SPNNs
in the circuit are modeled at time scales and firing rates relevant to spiking
biological neurons, but the simulation time scale is arbitrary. For the
parameter values, $g_{l}=\pi C_{m}/T$, $g_{c}=60\pi C_{m}/T$, $V_{l}=0$mV,
$C_{m}=10$pF, $L=\frac{1}{(2\pi/T)^{2}C_{m}}$, $W_{\text{spike}}=.3$pA,
$\tau_{d}=0.8T$, and $\tau_{s}=0.0T$. Simulations are for 300ms. The cycle
period, $T$, is 10ms with a frequency of 100Hz. Numerical integration is
performed with forward euler. The time step is 0.025ms. The period to time
step ratio is approximately 400 to 1. The first example is presented for 15
cycles or 150ms (Fig. 5). The second example is then presented for 15 cycles.
The output prediction is calculated as the unit which is furthest out of phase
with respect to the other output units within a 3 cycle window,
$\text{argmax}_{i}\ \mathbb{E}_{t_{spike,i}}\frac{1}{2}(|t_{\text{next
spike},j\neq i}-t_{spike,i}|+|t_{\text{last spike},j\neq i}-t_{spike,i}|)$.
Here, the expectation is taken with respect to the spike times of each
component $i$, within a three cycle window.

At this resolution, we notice a decrease in performance compared to the
original SPNN. Currently, simulating all MNIST and CIFAR test samples is
prohibitively expensive. We report qualitative results, leaving a quantitative
analysis of the SPNN performance versus numerical precision for future work.
The performance gap decreases as the time step decreases relative to the cycle
period and as the confidence in the output prediction increases. It is
promising that even at this resolution SPNN are robust to noise in the spike
times.

The phenomenon of pipelining emerges in SPNNs as an effect of a changing input
stimulus (Fig. 5). Pipelining is a technique used for reducing the latency and
increasing the resource utilization of computing hardware by connecting
processing elements in series and passing the output of one unit to the input
of the next unit in the pipeline. The signals propagate similar to travelling
waves in one direction.

##  5 Discussion

With the advent of powerful hardware solutions (TrueNorth [1], Loihi [3]) for
spiking neurons, the question arises how deep learning can be implemented with
spikes. The current approach is to implement real-valued activation values by
spike rate [8, 14]. Some of these approaches yield near state-of-the-art
results but require a large number of spikes and exhibit rather long response
latency (time to solution) due to the fact that the estimation of small rates
requires large integration windows. A potential remedy to this problem is to
compute with spike times [23], but current solutions are brittle and cannot be
scaled up to challenging machine learning problems.

Here, we present SPNNs employing spike timing codes and show them to be
successful at solving challenging machine learning tasks. In our SPNN model,
spike times represent phase angles of a complex number. Upon receiving input,
the network computes robustly by relaxing to a stable limit-cycle of periodic
spike trains, which is a fixed-point in the complex domain. Our work is
probably the first demonstration how a SNN employing a spike timing code can
learn the CIFAR-10 dataset.

Another problem that our model can address is how to implement gradient
learning in spiking neural networks using local signals. If deep learning is
implemented by a rate code, pre- and post-synaptic spikes need to be
integrated to estimate the activity values required in the gradient-based
learning rule. In contrast, in the SPNN the timing of spikes exactly encode
the continuous-valued phase variables, which are required in the learning
rule.

Our model is based on a phasor network with a complex-valued normalizing
activation function, first proposed in [4], constraining the phasor variables
to binary magnitude values of 0 or 1. It may be possible to extend our model
to combine rate and timing codes by representing continuous-valued magnitudes
of phasor variables in spike rates. In addition, there is other recent work on
deep learning in complex networks [28], proposing activation functions
different from the one used here. Future work should explore these other
activation functions in the context of SPNNs.

Our work has also interesting implications for neuroscience. As demonstrated
in the results, SPNNs exhibit time-locked, but not necessarily synchronous
spiking. This phenomenon is similar to polychronous spike patterns [9].
Therefore, SPNNs may serve as a model for how to compute with such spike
patterns, which also have been observed in neural activity [27].

Figure 5: Simulation of the circuit model SPNN on MNIST (left) and CIFAR
(right). Left: An example with label 5 is presented for 15 cycles or 150ms.
Immediately following, an example with label 8 is presented for 15 cycles. 5:
The membrane potential for a unit in the output layer. The unit is unstable
until cycle 7 when it beings to lock to the input signal. 5: Spike times for
the output layer. 5: Predicted class label over time. Right: An example with
label 0 is presented for 15 cycles or 150ms. Immediately following, an example
with label 3 is presented for 15 cycles. 5: The membrane potential for a unit
in the output layer. The unit is unstable until cycle 9 when it beings to lock
to the input signal. 5: Spike times for the output layer. 5: Predicted class
label over time.

## References

  * [1] F. Akopyan, J. Sawada, A. Cassidy, R. Alvarez-Icaza, J. Arthur, P. Merolla, N. Imam, Y. Nakamura, P. Datta, G.-J. Nam, et al.  Truenorth: Design and tool flow of a 65 mw 1 million neuron programmable neurosynaptic chip.  IEEE Transactions on Computer-Aided Design of Integrated Circuits and Systems, 34(10):1537–1557, 2015.
  * [2] L. Beaulieu-Laroche, E. H. Toloza, M.-S. van der Goes, M. Lafourcade, D. Barnagian, Z. M. Williams, E. N. Eskandar, M. P. Frosch, S. S. Cash, and M. T. Harnett.  Enhanced dendritic compartmentalization in human cortical neurons.  Cell, 175(3):643–651, 2018.
  * [3] M. Davies, N. Srinivasa, T.-H. Lin, G. Chinya, Y. Cao, S. H. Choday, G. Dimou, P. Joshi, N. Imam, S. Jain, et al.  Loihi: A neuromorphic manycore processor with on-chip learning.  IEEE Micro, 38(1):82–99, 2018.
  * [4] E. P. Frady and F. T. Sommer.  Robust computation with rhythmic spike patterns.  Proceedings of the National Academy of Sciences, 116(36):18050–18059, 2019.
  * [5] G. M. Georgiou and C. Koutsougeras.  Complex domain backpropagation.  IEEE transactions on Circuits and systems II: analog and digital signal processing, 39(5):330–334, 1992.
  * [6] D. F. Goodman and R. Brette.  The brian simulator.  Frontiers in neuroscience, 3:26, 2009.
  * [7] A. Hirose.  Continuous complex-valued back-propagation learning.  Electronics Letters, 28(20):1854–1855, 1992.
  * [8] E. Hunsberger and C. Eliasmith.  Spiking deep networks with lif neurons.  arXiv preprint arXiv:1510.08829, 2015.
  * [9] E. M. Izhikevich.  Polychronization: computation with spikes.  Neural computation, 18(2):245–282, 2006.
  * [10] E. M. Izhikevich.  Dynamical systems in neuroscience.  MIT press, 2007.
  * [11] D. P. Kingma and J. Ba.  Adam: A method for stochastic optimization.  arXiv preprint arXiv:1412.6980, 2014.
  * [12] A. Krizhevsky and G. Hinton.  Learning multiple layers of features from tiny images.  Technical report, CIFAR, 2009.
  * [13] Y. LeCun, L. Bottou, Y. Bengio, P. Haffner, et al.  Gradient-based learning applied to document recognition.  Proceedings of the IEEE, 86(11):2278–2324, 1998.
  * [14] J. H. Lee, T. Delbruck, and M. Pfeiffer.  Training deep spiking neural networks using backpropagation.  Frontiers in neuroscience, 10:508, 2016.
  * [15] M. London, A. Roth, L. Beeren, M. Häusser, and P. E. Latham.  Sensitivity to perturbations in vivo implies high noise and suggests rate coding in cortex.  Nature, 466(7302):123, 2010.
  * [16] Y. Matsumoto-Makidono, H. Nakayama, M. Yamasaki, T. Miyazaki, K. Kobayashi, M. Watanabe, M. Kano, K. Sakimura, and K. Hashimoto.  Ionic basis for membrane potential resonance in neurons of the inferior olive.  Cell reports, 16(4):994–1004, 2016.
  * [17] H. Mostafa.  Supervised learning based on temporal coding in spiking neural networks.  IEEE transactions on neural networks and learning systems, 29(7):3227–3235, 2018.
  * [18] D. E. Nikonov, G. Csaba, W. Porod, T. Shibata, D. Voils, D. Hammerstrom, I. A. Young, and G. I. Bourianoff.  Coupled-oscillator associative memory array operation for pattern recognition.  IEEE Journal on Exploratory Solid-State Computational Devices and Circuits, 1:85–93, 2015.
  * [19] A. J. Noest.  Phasor neural networks.  In Neural information processing systems, pages 584–591, 1988.
  * [20] P. O’Connor, E. Gavves, and M. Welling.  Training a network of spiking neurons with equilibrium propagation.  arxiv, 2018.
  * [21] A. Paszke, S. Gross, S. Chintala, G. Chanan, E. Yang, Z. DeVito, Z. Lin, A. Desmaison, L. Antiga, and A. Lerer.  Automatic differentiation in pytorch.  NeurIPS 2017, 2017.
  * [22] D. Roy, P. Panda, and K. Roy.  Learning spatio-temporal representations using spike-based backpropagation.  arxiv, 2018.
  * [23] B. Rueckauer and S.-C. Liu.  Conversion of analog to spiking neural networks using sparse temporal coding.  In 2018 IEEE International Symposium on Circuits and Systems (ISCAS), pages 1–5. IEEE, 2018.
  * [24] J. Sacramento, R. P. Costa, Y. Bengio, and W. Senn.  Dendritic cortical microcircuits approximate the backpropagation algorithm.  In Advances in Neural Information Processing Systems, pages 8735–8746, 2018.
  * [25] B. Scellier and Y. Bengio.  Equilibrium propagation: Bridging the gap between energy-based models and backpropagation.  Frontiers in computational neuroscience, 11:24, 2017.
  * [26] C. D. Schuman, T. E. Potok, R. M. Patton, J. D. Birdwell, M. E. Dean, G. S. Rose, and J. S. Plank.  A survey of neuromorphic computing and neural networks in hardware.  arXiv preprint arXiv:1705.06963, 2017.
  * [27] D. Tingley and G. Buzsáki.  Transformation of a spatial map across the hippocampal-lateral septal circuit.  Neuron, 98(6):1229–1242, 2018.
  * [28] C. Trabelsi, O. Bilaniuk, Y. Zhang, D. Serdyuk, S. Subramanian, J. F. Santos, S. Mehri, N. Rostamzadeh, Y. Bengio, and C. J. Pal.  Deep complex networks.  arXiv preprint arXiv:1705.09792, 2017.