arXiv:1612.03214v1 [cs.LG] 9 Dec 2016

Towards deep learning with spiking neurons in energy
based models with contrastive Hebbian plasticity

Thomas Mesnard
École Normale Supérieure
Paris, France
thomas.mesnard@ens.fr

Wulfram Gerstner
École Polytechnique Fédéral de Lausanne
1015 Lausanne, Switzerland
wulfram.gerstner@epfl.ch

Johanni Brea
École Polytechnique Fédéral de Lausanne
1015 Lausanne, Switzerland
johanni.brea@epfl.ch

Abstract
In machine learning, error back-propagation in multi-layer neural networks (deep
learning) has been impressively successful in supervised and reinforcement learning
tasks. As a model for learning in the brain, however, deep learning has long been
regarded as implausible, since it relies in its basic form on a non-local plasticity rule.
To overcome this problem, energy-based models with local contrastive Hebbian
learning were proposed and tested on a classification task with networks of rate
neurons. We extended this work by implementing and testing such a model with
networks of leaky integrate-and-fire neurons. Preliminary results indicate that it is
possible to learn a non-linear regression task with hidden layers, spiking neurons
and a local synaptic plasticity rule.

1

Introduction

Error back-propagation is an efficient and vastly used algorithm to train deep neural networks [7].
Although the architecture of deep learning was inspired by biological neural networks, adaptation
of this powerful algorithm for training deep spiking neural networks started to appear only recently
[1, 4].
There are several reasons why we are still waiting for a full “integration of deep learning and
neuroscience” [8]. In standard deep learning with feed-forward networks, rate-based neurons with
real-valued, non-negative output are layer-wise sequentially activated to compute the prediction of
the network given a certain input. After comparing this prediction with a target value, errors are
back-propagated in the reversed order through the layers, using the transposed of the feed-forward
weight matrices, but not the non-linearity of the neurons used in the forward pass to compute the
prediction. The update of weights and biases depends both on the neural activations of the forward
pass and the back-propagated errors. There is no notion of physical time in standard deep learning.
In contrast, biological neural networks consist of spiking neurons with binary outputs that work in
continuous time. It is unclear how such a network composed of real spiking neurons could switch
between non-linear forward passes and linear backward passes to implement error back-propagation.
Equilibrium propagation was recently proposed by Scellier and Bengio [14] as one way to introduce
physical time in deep learning and remove the necessity of different dynamics in forward and
backward pass. Their work combines related ideas like recurrent back-propagation [13], infinitesimal
perturbation of the outputs [10, 5] or contrastive Hebbian learning [9, 16] and moves these concepts
30th Conference on Neural Information Processing Systems (NIPS 2016), Barcelona, Spain.

closer to biology. In a recurrent network of rate-based neurons, whose dynamics is defined by
following the negative gradient of an energy function, Scellier and Bengio propose to relax the
network to a minimum energy state in the forward phase while fixing the rate of the input neurons at a
given value. The rate of the output neurons at the fixed point of the forward phase corresponds to the
prediction in standard deep learning. Moving the rate of the output neurons in direction of the target
value in the backward phase while keeping the input rates fixed, perturbs also the rate of the hidden
neurons, if backward connections exist. The key result of Scellier and Bengio is that supervised
learning with this network is possible with a simple contrastive Hebbian plasticity mechanism that
subtracts the correlation of firing rates at the fixed point of the forward phase from the correlation of
the firing rates after perturbation of the output.
Our contribution is to implement equilibrium propagation in a network of integrate-and-fire neurons
and test it on a regression problem. Our model differs slightly from the rate-based model in that
neither the input rates get explicitly fixed nor the output rates explicitly moved in direction of the
target. Instead, the input neurons receive a constant input current during both phases and the output
neurons are treated as two-compartment neurons that receive an extra somatic input current in the
backward phase [15]. The contrastive Hebbian plasticity mechanism is implemented with an estimate
of the firing rate by low-pass filtering the spike history with a large time constant.

2

Derivation of the learning rule for rate-based neurons

For completeness, we reproduce here the learning rule derivation described by Scellier and Bengio
[14].
Let us define the dynamics of a neural network by an energy function
X
X
X
2E(s; ŝ, θ, βx , βy ) =
s2i −
wij ρ(si )ρ(sj ) − 2
bi ρ(si ) +
i

+ βx

i

i6=j

X

X
(ŝi − si ) + βy
(ŝi − si )2 ,
2

i∈X

(1)

i∈Y

where si is the state of neuron i, ρ(si ) its firing rate for some non-linear function ρ, parameters
θ = (w, b), with symmetric connection strengths wij = wji and biases bi , Y and X are disjoint
subsets of neurons that may receive external input ŝi , if βx > 0 or βy > 0. Note that the network
does not need to have all-to-all connections. By setting some of the weights wij to zero, a multi-layer
architecture can be achieved.
The neural dynamics is given by
d
E(s; ŝ, θ, βx , βy )
dsi


X
= −si + ρ0 (si ) 
wij ρ(sj ) + bi  + Ii∈X βx (ŝi − si ) + Ii∈Y βy (ŝi − si ) ,

τ ṡi = −

(2)

j

with time constant τ and we used the indicator function Ii∈Y = 1 if i is in set Y and 0 otherwise. Note
that for the rectified-linear function, i.e. ρ(s) = s if s > 0 and ρ(s) = 0 otherwise, the derivative
ρ0 (s) = 1 for s > 0. Since for negative s the derivative ρ0 (s) = 0, s remains non-negative all the
time.
The lowest energy state s∗ given by
dE ∗
(s ; ŝ, θ, ∞, 0) = 0, ∀i
dsi

(3)

defines a map x 7→ y, with x = (s∗i1 , . . . , s∗iNx )ik ∈X and y = (s∗i1 , . . . , s∗iNy )ik ∈Y . We would like to
have a rule for changing the parameters w and b to implement an arbitrary map from x to y.
Let us define a cost function for a single pair of points ŝ and s
X
2C(ŝ, s) =
(ŝi − si )2
i∈Y

2

(4)

and a total cost function C(ŝ1 , s1 , . . . , ŝN , sN ) =

1
N

PN

µ=1

C(ŝµ , sµ ).

To find a learning rule, we look at the constraint optimization problem
dE µ µ
min C(ŝ1 , s1 , . . . , ŝN , sN ) subject to
(s ; ŝ , θ, ∞, 0) = 0, ∀µ, i
θ
dsi
and define the Lagrangian for a single data point ŝ by
X dE
L(s, λ, θ; ŝ) = C(ŝ, s) +
λi
(s; ŝ, θ, ∞, 0)
dsi
i
We minimize this Lagrangian by setting the derivatives with respect to s and λ to zero
dL ∗ ∗
(s , λ , θ; ŝ) = 0, ∀i
dλi
dL ∗ ∗
(s , λ , θ; ŝ) = 0, ∀i
dsi
and performing stochastic descent on the total cost by changing the parameters according to
∂L ∗ ∗
∆θi = −η
(s , λ , θ; ŝ)
∂θi
X
∂dE ∗
(s ; ŝ, θ, ∞, 0)
= −η
λ∗j
∂θ
i dsj
j

(5)

(6)

(7)
(8)

(9)
(10)

where η is a learning rate.
Solving Equation 7, we find that the state s∗ is simply given by the lowest energy state (c.f. Equation 3), i.e. it can be obtained by running the dynamics with input x (βx = ∞) and no target y
(βy = 0).
Equation 8 is a bit harder to solve for λ∗ . Expanding the definition of the Lagrangian it reads
X
dC
d2 E
(ŝ, s∗ ) +
λ∗j
(s∗ ; ŝ, θ, ∞, 0) = 0 .
(11)
dsi
ds
ds
i
j
j
The nice idea of Scellier and Bengio [14] is, to look at the fixed point sβ of the dynamics with weak
target input βy = β > 0 given by
dE β
(s ; ŝ, θ, ∞, β) = 0 , ∀β .
(12)
dsi
Since the left-hand side of this equation is a constant function in β, its derivative with respect to βy is
also zero and evaluated at βy = 0 we get
X dsβj d2 E
d2 E
∂dE
0=
=
+
dβy dsi
∂βy dsi
dβy dsi dsj
j

(13)

X dsβj d2 E
dC
∗
(ŝ, s ) +
(s∗ ; ŝ, θ, ∞, 0) ,
=
dsi
dβ
ds
ds
y
i
j
j
which has the same form as Equation 11 if we identify λ∗j =

(14)

dsβ
j
dβy .

Using this result in Equation 10, we can write
X
X dsβ d∂E

∂dE
d∂E
d
k
λ∗k
=
=
=−
ρ(si )ρ(sj )
(15)
∂wij dsk
dβy dsk ∂wij
dβy ∂wij
dβy
k
k


and approximate this by dβdy ρ(si )ρ(sj ) ≈ β1 ρ(sβi )ρ(sβj ) − ρ(s∗i )ρ(s∗j ) which leads to the contrastive Hebbian update rule

∆wij = η ρ(sβi )ρ(sβj ) − ρ(s∗i )ρ(s∗j ) ,
(16)
where the factor 1/β is absorbed into the learning rate η.
This learning rule can be implemented by iterating through:
3

A

B

Figure 1: Neuron model and network architecture. A Traces of the relevant quantities for a
target neuron. The synaptic weights are updated according to ∆wij = η(ri+ rj+ − ri− rj− ), where
ri+ = ri (1200 ms) and ri− = ri (600 ms), which can be implemented online by decreasing the
weights appropriately at 600 ms and increase them at 1200 ms. B Network architecture. Arrows
indicate all-to-all connectivity.
1. Select a data sample ŝ and relax the system to the lowest energy state with βy = 0 to obtain
y ∗ (Forward phase).
2. Subtract ηρ(s∗i )ρ(s∗j ) from the weights.
3. Set βy = β > 0 and let the system evolve for some time (Backward phase).1
4. Add ηρ(sβi )ρ(sβj ) to the weights.

3

Implementation with leaky integrate-and-fire neurons

To replace rate-based neurons with leaky integrate-and-fire neurons we introduce the somatic membrane potential ui of neuron i that evolves below a threshold θ as
X
τ u̇i (t) = −ui (t) + u0 + 1 − λi (t)
wij sj (t) + λi (t)RIi (t)
(17)
j

with time constant τ , resting potential u0 , post-synaptic response sj given by the dynamics τs ṡj =
P
(f )
(f )
−sj + upsp xj , with presynaptic spike trains xj (t) = t(f ) δ(t − tj ), where tj are the spike times
j

of neuron j, membrane resistance R and additional current input Ii (t), used in the backward phase to
nudge the firing rate of the neuron in direction of the target firing rate. The nudging factor
λi (t) =

RIi (t)
P
,
RIi (t) + j wij sj (t)

(18)

is important at the end of learning, when the predictions by the network match almost the target inputs
and additive instead of convex combination of network input and target input would lead to run-away
dynamics. The nudging factor can be motivated with divisive normalization [2] or an argument
involving conductance-based synapses [15]. The neuron spikes when its membrane potential reaches
threshold θ. The membrane potential is then set to a reset value ur and kept at this value for a
refractory period ∆. Afterwards the dynamics in Eq. 17 determines again the membrane potential.
To implement the contrastive Hebbian learning rule in Equation 16, the pre- and postsynaptic firing
rates are estimated in each synapse with hypothesized processes ri that low-pass filter the spike trains
with a large time constant τr , i.e.
τr ṙi (t) = −ri (t) + si (t) .
1

(19)

Scellier and Bengio [14] observe that relaxation to the fixed-point is not necessary in this second phase.

4

Figure 2: A non-linear regression task. Input
 values θ ∈ [0, 1] and ϕ ∈ [0, 1] are mapped
 to target
values x = cos(πϕ) + cos(π(ϕ + θ)) + 2 /4 and y = sin(πϕ) + sin(π(ϕ + θ)) + 2 /4.

The negative part of the weight update in Equation 16 is applied just before the network receives
target input to the output neurons. The positive part of the weight update is applied when the network
evolved for some time under the influence of the target input.
Figure 1A shows a trace of the relevant quantities for one target neuron.
The firing rate f of a leaky integrate-and-fire neuron with dynamics τ u̇ = −u + v, constant drive v,
reset potential ur , threshold θ and refractory period ∆ is given by
1

f (v) =
τ log



v−ur
v−θ



v > θ,

if

f (v) = 0

otherwise .

(20)

+∆

The rate model “closest”
 P to the spiking model uses this non-linearity, i.e. ρ(vi ) = f (vi ) with
vi = u0 + 1 − λi (t)
j wij sj (t) + λi (t)RIi (t). Even though the derivative of this non-linearity
is not piecewise constant as the one of the rectified-linear function, we did not see a significant
difference in simulations when using ρ0 (v) = 1 for v > 0 and ρ0 (v) = 0 otherwise together with the
leaky integrate-and-fire non-linearity f .
For the experiments described in the next section we used a network architecture with one hidden
layer (see Figure 1). We used only one hidden layer because the task is learnable with one, but
nothing prevents the use of multiple hidden layers. Bias terms are implemented with weights from
bias neurons that are active at a constant firing rate. Instead of clamping the input neurons to the
input values with βx → ∞, we did not allow any feedback from the other layers to the input layer
and provided the input as input currents Ii to the neurons in the input layer.
3.1

Results

In contrast to Scellier and Bengio [14], who tested the model on MNIST classification, we choose a
non-linear regression task (see Figure 2). Since energy-based models, like the Hopfield model [6],
are well known with parameter settings that create multiple local minima, we were curious to see, if
learning is also possible when local minima are undesirable.
3.1.1

Rate-based neurons

In Figure 3, we see that the task can be learned with a single layer of 400 rate-based neurons. Learning
is a bit slower when the firing rate f of the leaky integrate-and-fire neuron (c.f. Equation 20) is used
instead of the rectified linear function. As expected, the task cannot be learned without this hidden
layer.
In these simulations, weights were updated all in once after the backward phase by following
Equation 16. With small learning rates, the same results are expected, if the weights are updated
online, i.e. after the forward phase by subtracting ηρ(s∗i )ρ(s∗j ) and after the backward phase by
adding ηρ(sβi )ρ(sβj ).
5

Figure 3: Learning the regression task with networks of rate-based neurons. The error is the
estimated average Euclidian distance between prediction and target. The differential equations were
integrated with Euler method, time steps 1 ms, τ = 15 ms, τs = 10 ms, τr = 300 ms, u0 = 20,
ur = 0, θ = 20, upsp = 400, R = 40, ∆ = 5 ms,pIi ∈ [0, 1], duration of forward and backward
phase 600 ms, constant learning rates ηi = 0.1/ indegreei , where indegreei is the number of
presynaptic neurons of neuron i. Weights are randomly initialized to small values without symmetry
constraints. Note that the learning rule Equation 16 tends to symmetrize the weight matrix. We tested
two non-linearities: the rectified linear function (green and red) and the Liffi function defined in
Equation 20 (cyan). Preliminary results (not shown) indicate that learning with Liffi is slower, even
when learning rates are optimized for each non-linearity.
3.1.2

Spiking neurons

In the spiking implementation each real-valued input and output dimensions was represented by 20
neurons that received identical current inputs and the hidden layer is composed of 300 neurons.
Figure 4 shows a typical result with spiking neurons. Even though the target is not perfectly learned,
the predictions are much more accurate than without a hidden layer (c.f. Figure 3). In the simulations
with rate-based neurons we noticed the importance of converging close to the fixed point in the
forward phase. The poorer results with spiking neurons could be partially due to not having a good
enough estimate of the fixed point in the forward phase.

Figure 4: Learning the regression task with a network of leaky integrate-and-fire neurons. The
output values x and y are given by the average firing rate of 20 neurons. The differential equations
were integrated with Euler method, time steps 1 ms, τ = 15 ms, τs = 15 ms, τr = 100 ms, u0 = 20,
ur = 0, θ = 20, upsp = 400, R = 40, ∆ = 5 ms, Ip
i ∈ [0, 1], duration of forward and backward
phase 1000 ms, constant learning rates ηi = 5.10−5 / indegreei , where indegreei is the number of
presynaptic neurons of neuron i.

4

Discussion

We have implemented Equilibrium Propagation [14] with a multilayer network of leaky integrate-andfire neurons and demonstrated that the network can learn a non-linear regression task. The results
with the spiking networks are not as convincing as with the rate-based networks, which could be due
to non-optimal estimates of the fixed point in the forward phase.
6

The locality of the contrastive Hebbian plasticity rule and the error back-propagation by natural
recurrent dynamics make this approach appealing from the perspective of biological plausibility. The
plasticity mechanism requires quite precisely timed induction of anti-Hebbian and Hebbian plasticity.
Potential processes to implement this involve theta waves [11] or neuromodulators [12, 3], known for
modulating and even sign-reversing synaptic plasticity.
Equilibrium propagation is not the only attempt to implement error back-propagation in a more
biologically plausible way. A recent and compelling approach [4] postulates an implementation of
error back-propagation using the computational richness of multi-compartment neurons.
It is exciting to see the theoretical approaches that try to integrate deep learning and neuroscience.
Ultimately, experiments are needed to support or falsify these hypotheses.

Acknowledgments
References
[1] Yoshua Bengio, Benjamin Scellier, Olexa Bilaniuk, Joao Sacramento, and Walter Senn. Feedforward initialization for fast inference of deep generative networks is biologically plausible.
arXiv preprint arXiv:1606.01651, 2016.
[2] Matteo Carandini and David J Heeger. Normalization as a canonical neural computation. Nature
Reviews Neuroscience, 13(1):51–62, 2012.
[3] Nicolas Frémaux and Wulfram Gerstner. Neuromodulated spike-timing-dependent plasticity,
and theory of three-factor learning rules. Frontiers in neural circuits, 9, 2015.
[4] Jordan Guergiuev, Timothy P Lillicrap, and Blake A Richards. Deep learning with segregated
dendrites. arXiv preprint arXiv:1610.00161, 2016.
[5] John Hertz, Anders Krogh, Benny Lautrup, and Torsten Lehmann. Nonlinear backpropagation:
doing backpropagation without derivatives of the activation function. IEEE Transactions on
neural networks, 8(6):1321–1327, 1997.
[6] John J Hopfield. Neural networks and physical systems with emergent collective computational
abilities. Proceedings of the national academy of sciences, 79(8):2554–2558, 1982.
[7] Yann LeCun, Yoshua Bengio, and Geoffrey Hinton. Deep learning. Nature, 521(7553):436–444,
2015.
[8] Adam H Marblestone, Greg Wayne, and Konrad P Kording. Toward an integration of deep
learning and neuroscience. Frontiers in Computational Neuroscience, 10, 2016.
[9] Javier R Movellan. Contrastive hebbian learning in the continuous hopfield model. In
Connectionist models: Proceedings of the 1990 summer school, pages 10–17, 1991.
[10] Randall C O’Reilly. Biologically plausible error-driven learning using local activation differences: The generalized recirculation algorithm. Neural computation, 8(5):895–938, 1996.
[11] G Orr, G Rao, FP Houston, BL McNaughton, and Carol A Barnes. Hippocampal synaptic
plasticity is modulated by theta rhythm in the fascia dentata of adult and aged freely behaving
rats. Hippocampus, 11(6):647–654, 2001.
[12] Verena Pawlak, Jeffery R Wickens, Alfredo Kirkwood, and Jason ND Kerr. Timing is not
everything: neuromodulation opens the stdp gate. Spike-timing dependent plasticity, page 138,
2010.
[13] Fernando J Pineda. Generalization of back-propagation to recurrent neural networks. Physical
review letters, 59(19):2229, 1987.
[14] Benjamin Scellier and Yoshua Bengio. Towards a biologically plausible backprop. arXiv
preprint arXiv:1602.05179, 2016.
[15] Robert Urbanczik and Walter Senn. Learning by the dendritic prediction of somatic spiking.
Neuron, 81(3):521–528, 2014.
[16] Xiaohui Xie and H Sebastian Seung. Equivalence of backpropagation and contrastive hebbian
learning in a layered network. Neural computation, 15(2):441–454, 2003.

7

