LETTER

Communicated by Rajesh Rao

Optimal Spike-Timing-Dependent Plasticity for Precise Action
Potential Firing in Supervised Learning
Jean-Pascal Pﬁster
jean-pascal.pﬁster@epﬂ.ch

Taro Toyoizumi
taro@sat.t.u-tokyo.ac.jp

David Barber
david.barber@idiap.ch

Wulfram Gerstner
wulfram.gerstner@epﬂ.chLaboratory of Computational Neuroscience, School of Computer and Communication
Sciences and Brain-Mind Institute, Ecole Polytechnique Fédérale de Lausanne(EPFL), CH-1015 Lausanne, Switzerland

In timing-based neural codes, neurons have to emit action potentials
at precise moments in time. We use a supervised learning paradigm toderive a synaptic update rule that optimizes by gradient ascent the likelihood of postsynaptic ﬁring at one or several desired ﬁring times. We
ﬁnd that the optimal strategy of up- and downregulating synaptic efﬁcacies depends on the relative timing between presynaptic spike arrivaland desired postsynaptic ﬁring. If the presynaptic spike arrives before
the desired postsynaptic spike timing, our optimal learning rule predicts that the synapse should become potentiated. The dependence ofthe potentiation on spike timing directly reﬂects the time course of an excitatory postsynaptic potential. However, our approach gives no uniquereason for synaptic depression under reversed spike timing. In fact, thepresence and amplitude of depression of synaptic efﬁcacies for reversed
spike timing depend on how constraints are implemented in the optimization problem. Two different constraints, control of postsynaptic
rates and control of temporal locality, are studied. The relation of ourresults to spike-timing-dependent plasticity and reinforcement learning
is discussed.
1 Introduction
Experimental evidence suggests that precise timing of spikes is importantin several brain systems. In the barn owl auditory system, for example,
coincidence-detecting neurons receive volleys of temporally precise spikes
from both ears (Carr & Konishi, 1990). In the electrosensory system ofmormyrid electric ﬁsh, medium ganglion cells receive input at precisely
Neural Computation 18, 1318–1348 (2006)


C 2006 Massachusetts Institute of Technology

Optimal STDP in Supervised Learning1319

timed delays after electric pulse emission (Bell, Han, Sugawara, & Grant,
1997). Under the inﬂuence of a common oscillatory drive as present in the
rat hippocampus or olfactory system, the strength of a constant stimulus iscoded in the relative timing of neuronal action potentials (Hopﬁeld, 1995;
Brody & Hopﬁeld, 2003; Mehta, Lee, & Wilson, 2002). In humans, precise
timing of ﬁrst spikes in tactile afferents encodes touch signals at the ﬁngertips (Johansson & Birznieks, 2004). Similar codes have also been suggested
for rapid visual processing (Thorpe, Delorme, & Van Rullen, 2001), and forthe rat’s whisker response (Panzeri, Peterson, Schultz, Lebedev, & Diamond,
2001).
The precise timing of neuronal action potentials also plays an importantrole in spike-timing-dependent plasticity (STDP). If a presynaptic spikearrives at the synapse before the postsynaptic action potential, the synapse
is potentiated; if the timing is reversed, the synapse is depressed (Markram,
Lübke, Frotscher, & Sakmann, 1997; Zhang, Tao, Holt, Harris, & Poo, 1998;
Bi & Poo, 1998, 1999, 2001). This biphasic STDP function is reminiscent ofa temporal contrast or temporal derivative ﬁlter and suggests that STDP
is sensitive to the temporal features of a neural code. Indeed, theoreticalstudies have shown that given a biphasic STDP function, synaptic plasticity
can lead to a stabilization of synaptic weight dynamics (Kempter, Gerstner,& van Hemmen, 1999, 2001; Song, Miller, & Abbott, 2000; van Rossum,
Bi, & Turrigiano, 2000; Rubin, Lee, & Sompolinsky, 2001) while the neuron
remains sensitive to temporal structure in the input (Gerstner, Kempter,van Hemmen, & Wagner, 1996; Roberts, 1999; Kempter et al., 1999; Kistler& van Hemmen, 2000; Rao & Sejnowski, 2001; Gerstner & Kistler, 2002a).
While the relative ﬁring time of pre- and postsynaptic neurons, andhence temporal aspects of a neural code, play a role in STDP, it is less clear
whether STDP is useful to learn a temporal code. In order to elucidate
the computational function of STDP, we ask in this letter the following
question: What is the ideal form of an STDP function in order to generateaction potentials of the postsynaptic neuron with high temporal precision?
This question naturally leads to a supervised learning paradigm: thetask to be learned by the neuron is to ﬁre at a predeﬁned desired ﬁring time t des . Supervised paradigms are common in machine learning in
the context of classiﬁcation and prediction problems (Minsky & Papert,
1969; Haykin, 1994; Bishop, 1995), but have more recently also been studied
for spiking neurons in feedforward and recurrent networks (Legenstein,
Naeger, & Maass, 2005; Rao & Sejnowski, 2001; Barber, 2003; Gerstner, Ritz,& van Hemmen, 1993; Izhikevich, 2003). Compared to unsupervised or
reward-based learning paradigms, supervised paradigms on the level ofsingle spikes are obviously less relevant from a biological point, since it isquestionable what type of signal could tell the neuron about the “desired”
ﬁring time. Nevertheless, we think it is worth addressing the problem ofsupervised learning—ﬁrst, as a problem in its own right, and second, asa starting point of spike-based reinforcement learning (Xie & Seung, 2004;

1320

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstner

Seung, 2003). Reinforcement learning in a temporal coding paradigm implies that certain sequences of ﬁring times are rewarded, whereas others
are not. The “desired ﬁring times” are hence deﬁned indirectly via thepresence or absence of a reward signal. The exact relation of our supervised paradigm to reward-based reinforcement learning will be presented
in section 4. Section 2 introduces the stochastic neuron model and coding
paradigm, which are used to derive the results presented in section 3.
2 Model
2.1 Coding Paradigm. In order to explain our computational paradigm,
we focus on the example of temporal coding of human touch stimuli(Johansson & Birznieks, 2004), but the same ideas would apply analogously
to the other neuronal systems with temporal codes already mentioned (Carr& Konishi, 1990; Bell et al., 1997; Hopﬁeld, 1995; Brody & Hopﬁeld, 2003;
Mehta et al., 2002; Panzeri et al., 2001). For a given touch stimulus, spikes
in an ensemble of N tactile afferents occur in a precise temporal order. If
the same touch stimulus with identical surface properties and force vector
is repeated several times, the relative timing of action potentials is reliably
reproduced, whereas the spike timing in the same ensemble of afferents isdifferent for other stimuli (Johansson & Birznieks, 2004). In our model, we
assume that all input lines, labeled by the index j with 1 ≤ j ≤ N, converge
onto one or several postsynaptic neurons. We think of the postsynaptic
neuron as a detector for a given spatiotemporal spike pattern in the input.
The full spike pattern detection paradigm will be used in section 3.3. As a
preparation and ﬁrst steps toward the full coding paradigm, we also consider the response of a postsynaptic neuron to a single presynaptic spike(section 3.1) or to one given spatiotemporal ﬁring pattern (section 3.2).
2.2 Neuron Model. Let us consider a neuron i that is receiving input
from N presynaptic neurons. Let us denote the ensemble of all spikes ofNneuron j by x j = {t 1j , . . . , t j j }, where t kj denotes the time when neuron
j ﬁred its kth spike. The spatiotemporal spike pattern of all presynaptic
neurons 1 ≤ j ≤ N will be denoted by boldface x = {x1 , . . . , xN }.
f
A presynaptic spike elicited at time t j evokes an excitatory postsynapticfpotential (EPSP) of amplitude wi j and time course (t − t j ). For simplicity,
we approximate the EPSP time course by a double exponential,





s
s
(s) = 0 exp −− exp −
(s),
τm
τs

(2.1)

with a membrane time constant of τm = 10 ms and a synaptic time constantof τs = 0.7 ms, which yields an EPSP rise time of 2 ms. Here (s) denotes
the Heaviside step function with (s) = 1 for s > 0 and (s) = 0 otherwise.

Optimal STDP in Supervised Learning1321

We set 0 = 1.3 mV such that a spike at a synapse with wi j = 1 evokes an
EPSP with amplitude of approximately 1 mV. Since the EPSP amplitude isa measure of the strength of a synapse, we refer to wi j also as the efﬁcacy(or “weight”) of the synapse between neuron j and i.
Let us further suppose that the postsynaptic neuron i receives an additional input I (t) that could arise from either a second group of neurons
or from intracellular current injection. We think of the second input as ateaching input that increases the probability that the neuron ﬁres at or close
to the desired ﬁring time t des . For simplicity, we model the teaching input
as a square current pulse I (t) = I0 (t − t des + 0.5T)(t des + 0.5T − t)of amplitude I0 and duration T. The effect of the teaching current on themembrane potential isuteach (t) = ∞k(s)I (t − s)ds

(2.2)

0

with k(s) = k0 exp(−s/τm ), where k0 is a constant that is inversely proportional to the capacitance of the neuronal membrane.
In the context of the human touch paradigm discussed in section 2.1,
the teaching input could represent some preprocessed visual information(“object touched by ﬁngers starts to slip now”), feedback from muscle activity (“strong counterforce applied now”), cross-talk from other detector
neurons in the same population (“your colleagues are active now”), or unspeciﬁc modulatory input due to arousal or reward (“be aware—somethinginteresting happening now”).
In the context of training of recurrent networks (e.g., Rao & Sejnowski,2001), the teaching input consists of a short pulse of an amplitude thatguarantees action potential ﬁring.
The membrane potential of the postsynaptic neuron i (spike responsemodel; Gerstner & Kistler, 2002b) is inﬂuenced by the EPSPs evoked by
all afferent spikes of stimulus x, the “teaching” signal, and the refractoryfeffects generated by spikes ti of the postsynaptic neuronui (t|x, yti ) = urest +

N


wi j


f
t j ∈x j

j=1f

(t − t j ) +fη(t − ti ) + uteach (t),

(2.3)

f

ti ∈ytiwhere urest = −70 mV is the resting potential, yti = {ti1 , ti2 , . . . , tiF < t} is theset of postsynaptic spikes that occurred before t, and tiF always denotes thelast postsynaptic spike before t. On the right-hand side of equation 2.3, η(s)denotes the spike afterpotential generated by an action potential. We takesη(s) = η0 exp −
τm


(s),

(2.4)

1322

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstner

A

B
[Hz]

100

ν

post

0.5

prob. density

200

1g [kHz]

C

0

-60

-50

u [mV]

0
0

1

I0

2

0.04

I_0 = 1
I_0 = 1.5
I_0 = 2

0.02
0
0

50

100

150

200

ISI [ms]. (B) Firing rate of the postsynaptic
Figure 1: (A) Escape rate g(u) = ρ0 exp u−ϑuneuron as a function of the amplitude I0 of a constant stimulation current(arbitrary units). (C) Interspike interval (ISI) distribution for different input
currents.

where η0 < 0 is a reset parameter that describes how much the voltage
is reset after each spike (for the relation to integrate-and-ﬁre neurons, see
Gerstner & Kistler, 2002b). The spikes themselves are not modeled explicitly
but reduced to formal ﬁring times. Unless speciﬁed otherwise, we take
η0 = −5 mV.
In a deterministic version of the model, output spikes would be generated whenever the membrane potential ui reaches a threshold ϑ. In order toaccount for intrinsic noise and also for a small amount of synaptic noise generated by stochastic spike arrival from additional excitatory and inhibitory
presynaptic neurons that are not modeled explicitly, we replace the strict
threshold by a stochastic one. More precisely we adopt the following procedure (Gerstner & Kistler, 2002b). Action potentials of the postsynaptic
neuron i are generated by a point process with time-dependent stochastic
intensity ρi (t) = g(ui (t)) that depends nonlinearly on the membrane potential ui . Since the membrane potential in turn depends on both the input andthe ﬁring history of the postsynaptic neuron, we write:ρi (t|x, yti ) = g(ui (t|x, yti )).

(2.5)

We take an exponential
to describe the stochastic escape across threshold:g(u) = ρ0 exp u−ϑ
whereϑ = −50 mV is the formal threshold, u = 3 mVu
is the width of the threshold region and therefore tunes the stochasticityof the neuron, and ρ0 = 1/ms is the stochastic intensity at threshold (seeFigure 1). Other choices of the escape function g are possible with no qualitative change of the results. For u → 0, the model is identical to the deterministic leaky integrate-and-ﬁre model with synaptic current injection(Gerstner & Kistler, 2002b).

Optimal STDP in Supervised Learning1323

We note that the stochastic process, deﬁned in equation 2.5, is similar tobut different from a Poisson process since the stochastic intensity depends
on the set yt of the previous spikes of the postsynaptic neuron. Thus, theneuron model has some memory of previous spikes.
2.3 Stochastic Generative Model. The advantage of the probabilistic
framework introduced above via the noisy threshold is that it is possible
to describe the probability density1 Pi (y|x) of an entire spike train2 Y(t) =fft ∈y δ(t − ti ) (see appendix A for details):
i





  T


ρi (ti |x, yt f ) exp −ρi (s|x, ys )ds

Pi (y|x) = 

f

i

f

ti ∈y

= exp

 T

0


log(ρi (s|x, ys ))Y(s) − ρi (s|x, ys )ds .

(2.6)

0

Thus, we have a generative model that allows us to describe explicitly the
likelihood Pi (y|x) of emitting a set of spikes y for a given input x. Moreover,since the likelihood in equation 2.6 is a smooth function of its parameters,
it is straightforward to differentiate it with respect to any variable. Let us
differentiate Pi (y|x) with respect to the synaptic efﬁcacy wi j , since this is aquantity that we will use later,∂ log Pi (y|x)
=
∂wi j

 T
0


ρi (s|x, ys )f(s − t j )ds,
[Y(s) − ρi (s|x, ys )]ρi (s|x, ys )
f

(2.7)

t j ∈x jdg|u=ui (s|x,ys ) .
where ρi (s|x, ys ) = duIn this letter, we propose three different optimal models: A, B, and C (seeTable 1). The models differ in the stimulation paradigm and the speciﬁc
task of the neuron. In section 3, the task and hence the optimality criteria
are supposed to be given explicitly. However, the task in model C could
also be deﬁned indirectly by the presence or absence of a reward signal, as
discussed in section 4.1. The common idea behind all three approaches isthe notion of optimal performance. Optimality is deﬁned by an objectivefunction L that is directly related to the likelihood formula of equation 2.6
and that can be maximized by changes of the synaptic weights. Throughout

1 For simplicity, we denoted the set of postsynaptic spikes from 0 to T by y insteadof yT .
2 Capital Y is the spike train generated by the ensemble (lowercase) y.

1324

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstner

Table 1: Summary of the Optimality Criterion L for the Unconstrained Scenarios(Au , Bu , Cu ) and the Constrained Scenarios (Ac , Bc , Cc ).
Unconstrained Scenarios

Constrained Scenarios

Au —Postsynaptic spike imposed:
L Au = log(ρ(t des ))Bu —Postsynaptic spike imposed+ spontaneous activity: L Bu = log(ρ̄(t des ))imposed: 
Cu —Postsynaptic
 spike patterns
γ

i i
k M−1
L Cu = log
i Pi (y |x ) k=i Pi (0|x )Ac —No activity:
T
L Ac = L Au − 0 ρ(t)dt
Bc —Stabilized activity:
T
L Bc = L Bu − Tσ1 2 0 (ρ̄(t) − ν0 )2 dt
Cc —Temporal locality constraint:2
L Cc = L Cu , P = a δ  − T̃0

Notes: The constraint for scenario C is not included in the likelihood function L Cc itself,
but rather in the deconvolution with a matrix P that penalizes quadratically the terms
that are nonlocal in time. See appendix C for more details.

the article, this optimization is done by a standard technique of gradientascent,wi j = α

∂L
,
∂wi j

(2.8)

with a learning rate α. Since the three models correspond to three different
tasks, they have a slightly different objective function. Therefore, gradientascent yields slightly different strategies for synaptic update. In the following, we start with the simplest model with the aim of illustrating the basic
principles that generalize to the more complex models.
3 Results
In this section, we present synaptic updates rules derived by optimizing
the likelihood of postsynaptic spike ﬁring at some desired ﬁring time t des .
The essence of the argument is introduced in a particularly simple scenario,
where the neuron is stimulated by one presynaptic spike and the neuron isinactive except at the desired ﬁring time t des . This is the raw scenario thatis developed in several directions.
First, we may ask how the postsynaptic spike at the desired time t des isgenerated. The spike could simply be given by a supervisor. As always inmaximum likelihood approaches, we then optimize the likelihood that thisspike could have been generated by the neuron model (i.e., the generativemodel) given the known input. Or the spike could have been generated bya strong current pulse of short duration applied by the supervisor (teaching
input). In this case, the a priori likelihood that the generative model ﬁres
at or close to the desired ﬁring time is much higher. The two conceptual
paradigms give slightly different results, as discussed in scenario A.

Optimal STDP in Supervised Learning

1325

Second, we may, in addition to the spike at the desired time t des , allow
for other postsynaptic spikes generated spontaneously. The consequencesof spontaneous activity for the STDP function are discussed in scenario B.
Third, instead of imposing a single postsynaptic spike at a desired ﬁringtime t des , we can think of a temporal coding scheme where the postsynaptic
neuron responds to one (out of M) presynaptic spike pattern with a desired output spike train containing several spikes while staying inactive forthe other M − 1 presynaptic spike patterns. This corresponds to a patternclassiﬁcation task, which is the topic of scenario C.
Moreover, optimization can be performed in an unconstrained fashionor under some constraint. As we will see in this section, the speciﬁc form of
the constraint inﬂuences the results on STDP, in particular, the strength ofsynaptic depression for post-before-pre timing. To emphasize this aspect,
we discuss two constraints. The ﬁrst constraint is motivated by the observation that neurons have a preferred working point deﬁned by a typical meanﬁring rate that is stabilized by homeostatic synaptic processes (Turrigiano& Nelson, 2004). Penalizing deviations from a target ﬁring rate is the constraint that we will use in scenario B. For a very low target ﬁring rate, theconstraint reduces to the condition of “no activity,” which is the constraint
implemented in scenario A.
The second type of constraint is motivated by the notion of STDP itself:changes of synaptic plasticity should depend on the relative timing of preand postsynaptic spike ﬁring and not on other factors. If STDP is to be
implemented by some physical or chemical mechanisms with ﬁnite time
constants, we must require the STDP function to be local in time, that is,
the amplitude of the STDP function approaches zero for large time differences. This is the temporal locality constraint used in scenario C. Whilethe unconstrained optimization problems are labeled with the subscript
u (Au , Bu , Cu ), the constrained problems are marked by the subscript c(Ac , Bc , Cc ) (see Table 1).

3.1 Scenario A: One Postsynaptic Spike Imposed. Let us start with aparticularly simple model, which consists of one presynaptic neuron and
one postsynaptic neuron (see Figure 2A). Let us suppose that the task ofthe postsynaptic neuron i is to ﬁre a single spike at time t des in responseto the input, which consists of a single presynaptic spike at time t pre , thatis, the input is x = {t pre } and the desired output of the postsynaptic neuron isy = {t des }. Since there is only a single pre- and a single postsynaptic neuron
involved, we drop in this section the indices j and i of the two neurons.
3.1.1 Unconstrained Scenario Au : One Spike at t des . In this section, we
assume that the postsynaptic neuron has not been active in the recent past,that is, refractory effects are negligible. In this case, we have ρ(t|x, yt ) =ρ(t|x) because of the absence of previous spikes. Moreover, since there is

1326

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstner

A

B
0.4

A

∆w u

0.3
0.2
0.1

0

t pre

t des

T

0

-50

t

0

pre

des

-t

50

[ms]

Figure 2: (A) Scenario A: a single presynaptic neuron connected to a postsynaptic neuron with a synapse of weight w. (B) Optimal weight change given byequation 3.2 for scenario Au . This weight change is exactly the mirror image ofan EPSP.

only a single presynaptic spike (i.e., x = {t pre }), we write ρ(t|t pre ) instead of
ρ(t|x).
Since the task of the postsynaptic neuron is to ﬁre at time t des , we candeﬁne the optimality criterion L Au as the log likelihood of the ﬁring intensity
at time t des ,
L Au = log ρ(t des |t pre ) .

(3.1)

The gradient ascent on this function leads to the following STDP function,
w Au = α

∂ L Auρ  (t des |t pre ) des(t − t pre ),=α∂w
ρ(t des |t pre )

(3.2)

dg
|u=u(t|tpre ) . Since this optimal weight change w Au canwhere ρ  (t|t pre ) ≡ dube calculated for any presynaptic ﬁring time t pre , we get an STDP functionthat depends on the time difference t = t pre − t des (see Figure 2B). As we
can see directly from equation 3.2, the shape of the potentiation is exactly
a mirror image of an EPSP. This result is independent of the speciﬁc choiceof the function g(u).
The drawback of this simple model becomes apparent if the STDP function given by equation 3.2 is iterated over several repetitions of the experiment. Ideally, it should converge to an optimal solution given by w Au = 0
in equation 3.2. However, the optimal solution given by w Au = 0 is problematic: for t < 0, the optimal weight tends toward ∞, whereas for t ≥ 0,

Optimal STDP in Supervised Learning

1327

there is no unique optimal weight (w Au = 0, ∀w). The reason for this problem is that the model describes only potentiation and includes no mechanisms for depression.
3.1.2 Constrained Scenario Ac : No Other Spikes Than at t des . In order toget some insight into where the depression could come from, let us consider
a small modiﬁcation of the previous model. In addition to the fact that theneuron has to ﬁre at time t des , let us suppose that it should not ﬁre anywhere
else. This condition can be implemented by an application of equation 2.6
to the case of a single input spike x = {t pre } and a single output spikey = {t des }. In terms of notation, we set P(y|x) = P(t des |t pre ) and similarlyρ(s|x, y) = ρ(s|t pre , t des ) and use equation 2.6 to ﬁnd  TP(t des |t pre ) = ρ(t des |t pre ) exp −ρ(s|t pre , t des )ds .

(3.3)

0

Note that for s ≤ t des , the ﬁring intensity does not depend on t des ; hence,ρ(s|t pre , t des ) = ρ(s|t pre ) for s ≤ t des . We deﬁne the objective function L Ac as
the log likelihood of generating a single output spike at time t des , given asingle input spike at t pre . Hence, with equation 3.3,
L Ac = log(P(t des |t pre ))= log(ρ(t

des

|t

pre

)) −

 T

ρ(s|t pre , t des )ds,

(3.4)

0

and the gradient ascent w Ac = α∂ L Ac /∂w rule yields
w Ac = αρ  (t des |t pre ) des(t − t pre ) − αρ(t des |t pre ) Tρ  (s|t pre , t des )(s − t pre )ds.

0

(3.5)
Since we have a single postsynaptic spike at t des , equation 3.5 can directly beplotted as a STDP function. In Figure 3 we distinguish two different cases.
In Figure 3A we optimize the likelihood L Ac in the absence of any teaching
input. To understand this scenario, we may imagine that a postsynaptic
spike has occurred spontaneously at the desired ﬁring time t des . Applying
the appropriate weight update calculated from equation 3.5 will make sucha timing more likely the next time the presynaptic stimulus is repeated. The
reset amplitude η0 has only a small inﬂuence.
In Figure 3B, we consider a case where ﬁring of the postsynaptic spike at
the appropriate time was made highly likely by a teaching input of duration
T = 1 ms centered around the desired ﬁring t des . The form of the STDP
function depends on the amount η0 of the reset. If there is no reset η0 = 0, the

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstnerwithout teachingη0 = -10 mV
η0 = -5 mV
η0 = 0 mV

0.3

A

∆w c

0.2
0.1

B

with teaching

0.2
0.1

0

0

-0.1

-0.1

-0.2

-50

t

pre

0

des

-t

50

[ms]

η0 = -10 mV
η0 = -5 mV
η0 = 0 mV

0.3

A

A

∆w c

1328

-0.2

-50

t

pre

0

des

-t

50

[ms]

Figure 3: Optimal weight adaptation for scenario Ac given by equation 2.7 in
the absence of a teaching signal (A) and in the presence of a teaching signal(B). The weight change in the post-before-pre is governed by the spike afterpotential u AP (t) = η(t) + uteach (t). The duration of the teaching input is T = 1 ms.
The amplitude of the current I0 is chosen so that maxt uteach (t) = 5 mV. urest is
chosen such that the spontaneous ﬁring rate g(urest ) matches the desired ﬁring
rate 1/T: urest = u log( Tρ1 0 ) + θ −60 mV. The weight strength is w = 1.

STDP function shows strong synaptic depression of synapses that become
active after the postsynaptic spike. This is due to the fact that the teaching
input causes an increase of the membrane potential that decays back torest with the membrane time constant τm . Hence, the window of synaptic
depression is also exponential with the same time constant. Qualitatively
the same is true if we include a weak reset. The form of the depression
window remains the same, but its amplitude is reduced. The inverse of theeffect occurs only for strong reset to or below resting potential. A weakreset is standard in applications of integrate-and-ﬁre models to in vivo dataand is one of the possibilities for explaining the high coefﬁcient of variationof neuronal spike trains in vivo (Bugmann, Christodoulou, & Taylor, 1997;
Troyer & Miller, 1997).
A further property of the STDP functions in Figure 3 is a negative offset
for |t pre − t des | → ∞. The amplitude of ∞the offset can be calculated for w 0
and t > 0, that is, w0 −ρ  (urest ) 0 (s)ds. This offset is due to the factthat we do not want spikes at other times than t des . As a result, the optimal
weight w  (the solution of w Au = 0) should be as negative as possible(w  → −∞ or w  → w min in the presence of a lower bound) for t > 0 ort
0.
3.2 Scenario B: Spontaneous Activity. The constraint in scenario Acof having strictly no other postsynaptic spikes than the one at time t desmay seem artiﬁcial. Moreover, it is this constraint that leads to the negativeOptimal STDP in Supervised Learning

A

1329

B
0.4

j
B

∆w u

0.3
0.2
0.1

i

0

0

pre
tj

t

des

T

-50

t

0

pre

des

-t

50

[ms]

Figure 4: Scenario B. (A) N = 200 presynaptic neurons are ﬁring one after theother at time t j = jδt with δt = 1 ms. (B) The optimal STDP function of scenario
Bu .

offset of the STDP function discussed at the end of the previous paragraph.
In order to relax the constraint of no spiking, we allow in scenario B for areasonable spontaneous activity. As above, we start with an unconstrained
scenario Bu before we turn to the constrained scenario Bc .
3.2.1 Unconstrained Scenario Bu : Maximize the Firing Rate at t des . Let us
start with the simplest model, which includes spontaneous activity. Scenario
Bu is the analog of the model Au , but with two differences. First, we include
spontaneous activity in the model. Since ρ(t|x, yt ) depends on the spiking
history for any given trial, we have to deﬁne a quantity that is independentof the speciﬁc realizations y of the postsynaptic spike train. Second, insteadof considering only one presynaptic neuron, we consider N = 200 presynaptic neurons, each emitting a single spike at time t j = jδt, where δt = 1 ms(see Figure 4A). The input pattern will therefore be described by the set ofdelayed spikes x = {x j = {t j }, j = 1, . . . , N}. As long as we consider onlya single spatiotemporal spike pattern in the input, it is always possible torelabel neurons appropriately so that neuron j + 1 ﬁres after neuron j.
Let us deﬁne the instantaneous ﬁring rate ρ̄(t) that can be calculated byaveraging ρ(t|yt ) over all realizations of postsynaptic spike trains:ρ̄(t|x) = ρ(t|x, yt ) yt |x .

(3.6)

Here the notation · yt |x means taking the average over all possible conﬁguration of postsynaptic spikes up to t for a given input x. In analogy to a Poisson process, a speciﬁc spike train with ﬁring times yt = {ti1 , ti2 , . . . , tiF < t}is generated with probability P(yt |x) given by equation 2.6. Hence, the1330

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstneraverage · yt |x of equation 3.6 can be written as follows (see appendix B fornumerical evaluation of ρ̄(t)):ρ̄(t|x) =

 t
∞

1
F =0

F!

...

0

 t
0

ρ(t|x, yt )P(yt |x)dtiF , . . . , dti1 .

(3.7)

Analogous to model Au , we can deﬁne the quality criterion as the log
likelihood L Bu of ﬁring at the desired time t des :L Bu = log(ρ̄(t des |x)).

(3.8)

Thus, the optimal weight adaptation of synapse j is given by
w Bj u = α∂ ρ̄(t des |x)/∂w j,ρ̄(t des |x)(3.9)

where ∂ ρ̄(t|x)
is given by∂w j∂ ρ̄(t|x)∂= ρ̄ (t|x)(t − t j ) + ρ(t|x, yt )log P(yt |x)
,
∂w j
∂w j
yt |x
∂
log P(yt |x) is given by equation 2.7 and ρ̄  (t|x) =∂w j(3.10)

dg  .
du u=u(t|x,yt ) yt |x

Figure 4B shows that for our standard set of parameters, the differences
to scenario Au are negligible.
Figure 5A depicts the STDP function for various values of the parameteru at a higher postsynaptic ﬁring rate. We can see a small undershoot in
the pre-before-post region. The presence of this small undershoot can beunderstood as follows: enhancing a synapse of a presynaptic neuron thatﬁres too early would induce a postsynaptic spike that arrives before thedesired ﬁring time and because of refractoriness would therefore preventthe generation of a spike at the desired time. The depth of this undershootdecreases with the stochasticity of the neuron and increases with the amplitude of the refractory period (if η0 = 0, there is no undershoot). In fact,correlations between pre- and postsynaptic ﬁring reﬂect the shape of an
EPSP in the high-noise regime, whereas they show a trough for low noise(Poliakov, Powers, & Binder, 1997; Gerstner, 2001). Our theory shows thatthe pre-before-post region of the optimal plasticity function is a mirror
image of these correlations.
3.2.2 Constrained Scenario Bc : Firing Rate Close to ν0 . In analogy to model
Ac we introduce a constraint. Instead of imposing strictly no spikes attimes t = t des , we can relax the condition and minimize deviations of theOptimal STDP in Supervised Learning

A

1331

B
∆u = 0.5∆u = 1∆u = 3

0.8

B

B

∆w c

∆w u

0.6

σ = 4Hzσ = 6Hzσ = 8Hz

0.2

0.4

0.1

0.2
0

0
-0.2

-50

t

pre

0

des

-t

50

[ms]

-50

t

0

pre

des

-t

50

[ms]

Figure 5: (A) The optimal STDP functions of scenario Bu for different levelsof stochasticity described by the parameter u. The standard value (u =3 mV) is given by the solid line; decreased noise (u = 1 mV and u =0.5 mV) is indicated by dot-dashed and dashed lines, respectively. In the lownoise regime, enhancing a synapse that ﬁres slightly too early can prevent theﬁring at the desired ﬁring time t des due to refractoriness. To increase the ﬁring
rate at t des , it is advantageous to decrease the ﬁring probability time before t des .
Methods: For each value of u, the initial weight w0 is set such that the spontaneous ﬁring rate is ρ̄ = 30 Hz. In all three cases, w has been multiplied by uin order to normalize the amplitude of the STDP function. Reset: η0 = −5 mV.
(B) Scenario Bc . Optimal STDP function for scenario Bc given by equation 3.13 fora teaching signal of duration T = 1 ms. The maximal increase of the membranepotential after 1 ms of stimulation with the teaching input is maxt uteach (t) =5 mV. Synaptic efﬁcacies wi j are initialized such that u0 = −60 mV, which gives
a spontaneous rate of ρ̄ = ν0 = 5 Hz. Standard noise level: u = 3 mV.

instantaneous ﬁring rate ρ̄(t|x, t des ) from a reference ﬁring rate ν0 . This canbe done by introducing into equation 3.8 a penalty term PB given by
1 T (ρ̄(t|x, t des ) − ν0 )2
dt .
PB = exp −
T 0
2σ 2

(3.11)For small σ , deviations from the reference rate yield a large penalty. For
σ → ∞, the penalty term has no inﬂuence. The optimality criterion is acombination of a high ﬁring rate ρ̄ at the desired time under the constraintof small deviations from the reference rate ν0 . If we impose the penalty as amultiplicative factor and take as before the logarithm, we get
L Bc = log ρ̄(t des |x)PB .

(3.12)

1332

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstner

x1

y1

Figure 6: Scenario C. N presynaptic neurons are fully connected to M postsynaptic neurons. Each postsynaptic neuron is trained to respond to a speciﬁc
input pattern and not respond to M − 1 other patterns as described by theobjective function of equation 3.14.

Hence the optimal weight adaptation is given by

w Bj c = α∂ ρ̄(t des |x)/∂w jα−Tσ 2ρ̄(t des |x) T(ρ̄(t|x, t des ) − ν0 )0∂
ρ̄(t|x, t des )dt.
∂w j
(3.13)

Since in scenario B, each presynaptic neuron j ﬁres exactly once at time t j =jδt and the postsynaptic neuron is trained to ﬁre at time t des , we can interpretthe weight adaptation w Bj c of equation 3.13 as an STDP function w Bc
that depends on the time difference t = t pre − t des . Figure 5 shows this
STDP function for different values of the free parameter σ of equation 3.11.
The higher the standard deviation σ , the less effective is the penalty term.
In the limit of σ → ∞, the penalty term can be ignored, and the situation is
identical to that of scenario Bu .
3.3 Scenario C: Pattern Detection
3.3.1 Unconstrained Scenario Cu : Spike Pattern Imposed. This last scenario
is a generalization of scenario Ac . Instead of restricting the study to a single
pre- and postsynaptic neuron, we consider N presynaptic neurons and M
postsynaptic neurons (see Figure 6). The idea is to construct M independent
detector neurons. Each detector neuron i = 1, . . . , M, should respond best
to a speciﬁc prototype stimulus, say xi , by producing a desired spike train
yi , but should not respond to other stimuli, yi = 0, ∀xk , k = i (see Figure 7).
The aim is to ﬁnd a set of synaptic weights that maximizes the probability

Optimal STDP in Supervised Learning

A

B
k

400

200

x

x

i

400

0

1333

0

100

0

200

k

500

0

100

0

200

200

0

100

200

100

200

0.1

ρk

ρi

100

500

0.4
0.2
0

0

1000

y

y

i

1000

0

200

0

100

0.05

200

0

0

time [ms]time [ms]Figure 7: Pattern detection after learning. (Top) The left raster plot represents
the input pattern the ith neuron has to be sensitive to. Each line corresponds toone of the N = 400 presynaptic neurons. Each dot represents an action potential.
The right ﬁgure represents one of the patterns the ith neuron should not respondto. (Middle) The left raster plot corresponds to 1000 repetitions of the output of
neuron i when the corresponding pattern xi is presented. The right plot is the
response of neuron i to one of the patterns it should not respond to. (Bottom)
The left graph represents the probability density of ﬁring when pattern xi ispresented. This plot can be seen as the PSTH of the middle graph. Arrows
indicate the supervised timing neuron i learned. The right graph describes theprobability density of ﬁring when pattern xk is presented. Note the different
scales of vertical axis.

that neuron i produces yi when xi is presented and produces no output
when xk , k = i is presented. Let the likelihood function L Cu beL Cu = log MPi (yi |xi )
i=1



M

Pi (0|xk )

γ
M−1



(3.14)

k=1,k=iwhere Pi (yi |xi ) (see equation 2.6) is the probability that neuron i produces
the spike train yi when the stimulus xi is presented. The parameter γ

1334

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstnercharacterizes the relative importance of the patterns that should not belearned compared to those that should be learned. We get

L Cu =M
log(Pi (yi |xi )) + γ log(Pi (0|xk )) xk =xi ,(3.15)i=1M1where the notation ·xk =xi ≡ M−1
k=i means taking the average over all
patterns other than xi . The optimal weight adaptation yields

wiCj = α







∂
∂
log Pi (yi |xi ) + αγlog Pi (0|xk )
.
∂wi j
∂wi j
xk =xi(3.16)The learning rule of equation 3.16 gives the optimal weight change foreach synapse and can be evaluated after presentation of all pre- and postsynaptic spike patterns; it is a “batch” update rule. Since each pre- andpostsynaptic neuron emits many spikes in the interval [0, T], we cannot directly interpret the result of equation 3.16 as a function of the time differencet = t pre − t des as we did in scenario A or B.
Ideally, we would like to write the total weight change of the optimal
rule given by equation 3.16 as a sum of contributions
wiCj =



WCu (t pre − t des ),(3.17)t pre ∈xij
t des ∈yi

where WCu (t pre − t des ) is an STDP function and the summation runs overall pairs of pre- and postsynaptic spikes. The number of pairs of pre- andpostsynaptic spikes with a given time shift is given by the correlation function, which is best deﬁned in discrete time. We assume time steps of duration
δt = 0.5 ms. Since the correlation will depend on the presynaptic neuron jand the postsynaptic neuron i under consideration, we introduce a new index, k = N(i − 1) + j. We deﬁne the correlation in discrete time by its matrix
elements Ck that describe the correlation between the presynaptic spiketrain Xij (t) and the postsynaptic spike train Yi (t − T0 + δt). For example,
C3 = 7 implies that seven spike pairs of presynaptic neuron j = 3 withpostsynaptic neuron i = 1 have a relative time shift of T0 − δt. With thisdeﬁnition, we can rewrite equation 3.17 in vector notation (see section C.1
for more details) as!

wC = CWCu ,(3.18)

Optimal STDP in Supervised Learning

A

1335

B
0.4
0.3

∆W c

∆W u

0.2

C

C

a = 0.04
a = 0.4

0.2
0.1

0.1

0

0
-50

t

0

pre

des

-t

50

[ms]

-0.1

-50

t

0

pre

des

-t

50

[ms]

Figure 8: (A) Optimal weight change for scenario Cu . In this case, no locality
constraint is imposed, and the result is similar to the STDP function of scenario
Ac (with η0 = 0 and uteach (t) = 0) represented on Figure 3. (B) Optimal weightchange for scenario Cc as a function of the locality constraint characterized by
a . The stronger the importance of the locality constraint, the narrower is the
spike-spike interaction. For A and B, M = 20, η0 = −5 mV. The initial weights
wi j are chosen so that the spontaneous ﬁring rate matches the imposed ﬁring
rate.
C
C
C
where wC = (w11, . . . , w1N
, w21, . . . , wCMN )T is the vector containing all the optimal weight change given by equation 3.16 and WCu isthe vector containing the discretized STDP function with components
WCu = WCu (−T0 + δt) for 1 ≤  ≤ 2T̃0 with T̃0 = T0 /δ. In particular,
the center of the STDP function (t pre = t des ) corresponds to the index  = T̃0 .
!

The symbol = expresses the fact that we want to ﬁnd WCu such that
wC is as close as possible to CWCu . By taking the pseudo-inverse
C + = (C T C)−1 C T of C, we can invert equation 3.18 and get
WCu = C + wC .

(3.19)The resulting STDP function is plotted in Figure 8A. As it was the case forthe scenario Au , the STDP function exhibits a negative offset. In addition
to the fact the postsynaptic neuron i should not ﬁre at other times thanthe ones given by yi , it should also not ﬁre whenever pattern xk , k = i ispresented. The presence of the negative offset is due to those two factors.
3.3.2 Constrained Scenario Cc : Temporal Locality. In the previous paragraph, we obtained a STDP function with a negative offset. This negative
offset does not seem realistic because it implies that the STDP function is
not localized in time. In order to impose temporal locality (ﬁnite memory

1336

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstner

A

B
0.4
M=20
M=60
M=100

0.2

C

C

∆W c

0.1

w=0.5w0
w=w0
w=1.5w0

0.3

∆W c

0.2

0.1
0

-0.1

0

-0.2
-0.1

-50

t

0

pre

des

-t

50

[ms]

-50

t

pre

0

des

-t

50

[ms]

Figure 9: (A) Optimal STDP function as a function of the number of input
patterns M. (a = 0.04, N = 400). (B) Optimal weight change as a function of
the weight w. If the weights are small (dashed line), potentiation dominates,
whereas if they are big (dotted line), depression dominates.

span of the learning rule), we modify equation 3.19 in the following way(see section C.2 for more details):
WCc = (C T C + P)−1 C T wC ,(3.20)

where P is a diagonal matrix that penalizes nonlocal terms. In this article, wetake a quadratic suppression of terms that are nonlocal in time. With respect
to a postsynaptic spike at t des , the penalty term is proportional to (t − t des )2 .
In matrix notation and using our convention that the postsynaptic spikecorresponds to  = T̃0 , we have:2P = a δ  − T̃0 .

(3.21)The resulting STDP functions for different values of a are plotted in
Figure 8B. The higher the parameter a , the more nonlocal terms are penalized, the narrower is the STDP function.
Figure 9A shows the STDP functions for various number of patterns
M. No signiﬁcant change can be observed for different numbers of input
patterns M. This is due to the appropriately chosen normalization factor
1/(M − 1) in the exponent of equation 3.14.
The target spike trains yi have a certain number of spikes during thetime window T; they set a target value for the mean rate. Let ν post = T1M ×
M T ii=1 0 y (t)dt be the imposed ﬁring rate. Let w0 denote the amplitude of
the synaptic strength such that the ﬁring rate ρ̄w0 given by those weights
is identical to the imposed ﬁring rate: ρ̄w0 = ν post . If the actual weights areOptimal STDP in Supervised Learning

a=0

B

C

0

-0.3

0

-0.3
-0.3

0

opt

∆w

0.3

a = 0.4
0.3

rec

rec

0.3

∆w

∆w

rec

0.3

a = 0.04

∆w

A

1337

0

-0.3
-0.3

0

opt

∆w

0.3

-0.3

0

opt

0.3

∆w

Figure 10: Correlation plot between the optimal synaptic weight changewopt = wCu and the reconstructed weight change wrec = CWCc using
the temporal locality constraint. (A) No locality constraint; a = 0. Deviations
from the diagonal are due to the fact that the optimal weight change given byequation 3.16 cannot be perfectly accounted for the sum of pair effects. The mean
deviations are given by equation C.7. (B) A weak locality constraint (a = 0.04) almost does not change the quality of the weight change reconstruction. (C) Strong
locality constraint (a = 0.4). The horizontal lines arise since most synapses are
subject to a few strong updates induced by pairs of pre- and postsynaptic spiketimes with small time shifts.

smaller than w0 , almost all the weights should increase, whereas if they arebigger than w0 , depression should dominate (see Figure 9B). Thus, the exact
form of the optimal STDP function depends on the initial weight value w0 .
Alternatively, homeostatic process could ensure that the mean weight value
is always in the appropriate regime.
In equations 3.17 and 3.18, we imposed that the total weight change
should be generated as a sum over pairs of pre- and postsynaptic spikes.
This is an assumption that has been made in order to establish a link tostandard STDP theory and experiments where spike pairs have been in thecenter of interest (Gerstner et al., 1996; Kempter et al., 1999 ; Kistler & vanHemmen, 2000; Markram et al., 1997; Bi & Poo, 1998; Zhang et al., 1998). It
is, however, clear by now that the timing of spike pairs is only one of several
factors contributing to synaptic plasticity. We therefore asked how much wemiss if we attribute the optimal weight changes calculated in equation 3.16
to spike pair effects only. To answer this question, we compared the optimal
weight change wiCj from equation 3.16 with that derived from the pairCc pre− t des ) with or without
based STDP rule wirec
j =t pre ∈xij
t des ∈yi W (tlocality constraint, that is, for different values of the locality parameter(a = 0, 0.04, 0.4) (see Figure 10). More precisely, we simulate M = 20 detector neurons, each having N = 400 presynaptic inputs, so each subplotof Figure 10 contains 8000 points. Each point in a graph corresponds tothe optimal change of one weight for one detector neuron (x-axis) compared to the weight change of the same weight due to pair-based STDP

1338

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstner(y-axis). We found that in the absence of a locality constraint, the pair-wise
contributions are well correlated with the optimal weight changes. Withstrong locality constraints, the quality of the correlation drops signiﬁcantly.
However, for a weak locality constraint that corresponds to an STDP function with reasonable potentiation and depression regimes, the correlationof the pair-based STDP rule with the optimal update is still good. Thissuggests that synaptic updates with an STDP function based on pairs of
pre- and postsynaptic spikes is close to optimal in the pattern detection
paradigm.

4 Discussion
4.1 Supervised versus Unsupervised and Reinforcement Learning.
Our approach is based on the maximization of the probability of ﬁring
at desired times t des with or without constraints. From the point of viewof machine learning, this is a supervised learning paradigm implemented
as a maximum likelihood approach using the spike response model withescape noise as a generative model. Our work can be seen as a continuoustime extension of the maximum likelihood approach proposed in Barber(2003).
The starting point of all supervised paradigms is the comparison of a desired output with the actual output a neuron has, or would have, generated.
The difference between the desired and actual output is then used as thedriving signal for synaptic updates in typical model approaches (Minsky& Papert, 1969; Haykin, 1994; Bishop, 1995). How does this compare to experimental approaches? Experiments focusing on STDP have been mostly
performed in vitro (Markram et al., 1997; Magee & Johnston, 1997; Bi &Poo, 1998). Since in typical experimental paradigms, ﬁring of the postsynaptic neuron is enforced by strong pulses of current injection, the neuron
is not in a natural unsupervised setting; but the situation is also not fully
supervised, since there is never a conﬂict between the desired and actual
output of a neuron. In one of the rare in vivo experiments to STDP (Frégnac,
Shulz, Thorpe, & Bienenstock, 1988, 1992), the spikes of the postsynaptic
neuron are also imposed by current injection. Thus, a classiﬁcation of STDP
experiments in terms of supervised, unsupervised, or reward based is not
as clear-cut as it may seem at a ﬁrst glance.
From the point of view of neuroscience, paradigms of unsupervised orreinforcement are probably much more relevant than the supervised scenario discussed here. However, most of our results from the supervised
scenario analyzed in this article can be reinterpreted in the context of reinforcement learning following the approach proposed by Xie and Seung(2004). To illustrate the link between reinforcement learning and supervised
learning, we deﬁne a global reinforcement signal R(x, y) that depends on the
spike timing of the presynaptic neurons x and the postsynaptic neuron y.

Optimal STDP in Supervised Learning1339

The quantity optimized in reinforcement learning is the expected rewardRx,y averaged over all pre- and postsynaptic spike trains:Rx,y =

R(x, y)P(y|x)P(x).

(4.1)

x,y

If the goal of learning is to maximize the expected reward, we can deﬁne alearning rule that achieves this goal by changing synaptic efﬁcacies in the
direction of the gradient of the expected reward Rx,y :∂ log P(y|x)
wx,y = α R(x, y)
,
∂w
x,y

(4.2)

∂ log P(y|x)

is the quantity we discussed
where α is a learning parameter and∂w
in this article. Thus, the quantities optimized in our supervised paradigm
re-appear naturally in a reinforcement learning paradigm.
For an intuitive interpretation of the link between reinforcement learning, and supervised learning, consider a postsynaptic spike that (spontaneously) occurred at time t0 . If no reward is given, no synaptic change takes
place. However, if the postsynaptic spike at t0 is linked to a rewarding situation, the synapse will try to recreate in the next trial a spike at the same time,that is, t0 has the role of the desired ﬁring time t des introduced in this article.
Thus, the STDP function with respect to a postsynaptic spike at t des derived
in this article can be seen as the spike timing dependence that maximizes
the expected reward in a spike-based reinforcement learning paradigm.
4.2 Interpretation of STDP Function. Let us now summarize and discuss our results in a broader context. In all three scenarios, we found an
STDP function with potentiation for pre-before-post timing. Thus, this result
is structurally stable and independent of model details. However, depression for post-before-pre timing does depend on model details.
In scenario A, we saw that the behavior of the post-before-pre region isdetermined by the spike afterpotential (see Table 2 for a result summaryof the three models). In the presence of a teaching input and ﬁring rate
constraints, a weak reset of the membrane potential after the spike meansthat the neuron effectively has a depolarizing spike after potential (DAP).
In experiments, DAPs have been observed by Feldman (2000), Markram
et al. (1997), and Bi and Poo (1998) for strong presynaptic input. Other
studies have shown that the level of depression does not depend on thepostsynaptic membrane potential (Sjöström, Turrigiano, & Nelson, 2001).
In any case, a weak reset (i.e., to a value below threshold rather than to
the resting potential) is consistent with the ﬁndings of other researchersthat used integrate-and-ﬁre models to account for the high coefﬁcient of

1340

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstner

Table 2: Main Results for Each Scenario.
Unconstrained Scenarios
Au —pre-before-post: LTP ∼ EPSP
Bu —pre-before-post: LTP/LTD ∼ reverse
correlation
Cu —pre-before-post: LTP ∼ EPSP,LTD ∼ background patterns

Constrained Scenarios
Ac —post-before-pre: LTD (or LTP) ∼ spikeafterpotential
Bc —post-before-pre: LTD ∼ increased
ﬁring rate
Cc —post-before-pre: LTD ∼ backgroundpatterns ∼ temporal locality

variation of spike trains in vivo (Bugmann et al., 1997; Troyer & Miller,
1997).
In the presence of spontaneous activity (scenario B), a constraint on the
spontaneous ﬁring rate causes the optimal weight change to elicit a depression of presynaptic spikes that arrive immediately after the postsynaptic
one. In fact, the reason for the presence of the depression in scenario Bc is
directly related to the presence of a DAP caused by the strong teaching stimulus. In both scenarios A and B, depression occurs in order to compensate
the increased ﬁring probability due to the DAP.
In scenario C, it has been shown that the best way to adapt the weights(in a task where the postsynaptic neuron has to detect a speciﬁc input
pattern among others) can be described as an STDP function. This task issimilar to the one in Izhikevich (2003) in the sense that a neuron is designed
to be sensitive to a speciﬁc input pattern, but different since our work does
not assume any axonal delays. The depression part in this scenario arises
from a locality constraint. We impose that weight changes are explained bya sum of pair-based STDP functions.
There are various ways of deﬁning objective functions, and we haveused three different objective functions in this article. The formulation ofan objective function gives a mathematical expression of the functionalrole we assign to a neuron. The functional role depends on the type ofcoding (temporal coding or rate coding) and hence on the information thepostsynaptic neurons will read out. The functional role also depends on
the task or context in which a neuron is embedded. It might seem thatdifferent tasks and coding schemes could thus give rise to a huge numberof objective functions. However, the reinterpretation of our approach in thecontext of reinforcement learning provides a unifying viewpoint: even if
the functional role of some neurons in a speciﬁc region of the brain canbe different from other neurons of a different region, it is still possible tosee the different objective functions as different instantiations of the sameunderlying concept: the maximization of the reward, where the reward istask speciﬁc.
More speciﬁcally, all objective functions used in this letter maximized
the ﬁring probability at a desired ﬁring time t des , reﬂecting the fact that

Optimal STDP in Supervised Learning1341

in the framework of timing-based codes, the task of a neuron is to ﬁreat precise moments in time. With a different assumption on the neuron’srole on signal processing, different objective functions need to be used.
An extreme case is a situation where the neuron’s task is to avoid ﬁring
at time t des . A good illustration is given by the experiments done in theelectrosensory lobe (ELL) of the electric ﬁsh (Bell et al., 1997). These cells
receive two sets of input: the ﬁrst one contains the pulses coming from theelectric organ, and the second input conveys information about the sensory
stimulus. Since a large fraction of the sensory stimulus can be predictedby the information coming from the electric organ, it is computationallyinteresting to subtract the predictable contribution and focus on only theunpredictable part of the sensory stimulus. In this context, a reasonabletask would be to ask the neuron not to ﬁre at time t des where t des is the time
where the predictable simulation arrives, and this task could be deﬁned
indirectly by an appropriate reward signal. An objective function of thistype would, in the end, reverse the sign of the weight change of the causal
part (LTD for the pre-before-post region), and this is precisely what is seen
experimentally (Bell et al., 1997).
In our framework, the deﬁnition of the objective function is closely related to the neuronal coding. In scenario C, we postulate that neurons emita precise spike train whenever the “correct” input is presented and aresilent otherwise. This coding scheme is clearly not the most efﬁcient one.
Another possibility is to impose postsynaptic neurons to produce a speciﬁc but different spike train for each input pattern, and not only for the“correct” input. Such a modiﬁcation of the scenario does not dramaticallychange the results. The only effect is to reduce the amount of depression
and increase the amount of potentiation.
4.3 Optimality Approaches versus Mechanistic Models. Theoretical
approaches to neurophysiological phenomena in general, and to synaptic plasticity in particular, can be roughly grouped into three categories:
biophysical models that aim at explaining the STDP function from principles of ion channel dynamics and intracelluar processes (Senn, Tsodyks,& Markram, 2001; Shouval, Bear, & Cooper, 2002; Abarbanel, Huerta, &Rabinovich, 2002; Karmarkar & Buonomano, 2002); mathematical models
that start from a given STDP function and analyze computational principles such as intrinsic normalization of summed efﬁcacies or sensitivity tocorrelations in the input (Kempter et al., 1999; Roberts, 1999; Roberts & Bell,2000; van Rossum et al., 2000; Kistler & van Hemmen, 2000; Song et al., 2000;
Song & Abbott, 2001; Kempter et al., 2001; Gütig, Aharonov, Rotter, & Sompolinsky, 2003); and models that derive “optimal” STDP properties for agiven computational task (Chechik, 2003; Dayan & Häusser, 2004; Hopﬁeld& Brody, 2004; Bohte & Mozer, 2005; Bell & Parra, 2005; Toyoizumi, Pﬁster,
Aihara, & Gerstner, 2005a, 2005b). Optimizing the likelihood of postsynaptic
ﬁring in a predeﬁned interval, as we did in this letter, is only one possibility

1342

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstneramong others of introducing concepts of optimality (Barlow, 1961; Atick &Redlich, 1990; Bell & Sejnowski, 1995) into the ﬁeld of STDP. Chechik (2003)
uses concepts from information theory but restricts his study to the classiﬁcation of stationary patterns. The paradigm considered in Bohte and Mozer(2005) is similar to our scenario Bc , in that they use a fairly strong teaching
input to make the postsynaptic neuron ﬁre. Bell and Parra (2005) and Toyoizumi et al. (2005a) are also using concepts from information theory, butthey are applying them to the pre- and postsynaptic spike trains. The workof Toyoizumi et al. (2005a) is a clearcut unsupervised learning paradigm
and hence distinct from our approach. Dayan and Häusser (2004) use concepts of optimal ﬁlter theory but are not interested in precise ﬁring of thepostsynaptic neuron. The work of Hopﬁeld and Brody (2004) is similar to
our approach in that it focuses on recognition of temporal input patterns,
but we are also interested in triggering postsynaptic ﬁring with precise
timing. Hopﬁeld and Brody emphasize the repair of disrupted synapses
in a network that has previously acquired its function of temporal patterndetector.
Optimality approaches such as ours will never be able to make strictpredictions about the properties of neurons or synapses. Optimality criteria
may, however, help to elucidate computational principles and provide insights into potential tasks of electrophysiological phenomena such as STDP.
Appendix A: Probability Density of a Spike Train
The probability density of generating a spike train yt = {ti1 , ti2 , . . . , tiF < t}with the stochastic process deﬁned by equation 2.5 can be expressed asfollows,
P(yt ) = P(ti1 , . . . , tiF )R(t|yt ),(A.1)

where P(ti1 , . . . , tiF ) is the probability density of having F spikes at times
tti1 , . . . , tiF and R(t|yt ) = exp(− t F ρ(t  |yt )dt  ) corresponds to the probabilityiof having no spikes from tiF to t. Since the joint probability P(ti1 , . . . , tiF ) canbe expressed as a product of conditional probabilities,
F

P(ti1 , . . . , tiF ) = P(ti1 ) f f −1
P ti |ti , . . . , ti1 ,f =2Equation A.1 becomesP(yt ) = ρ(ti1 |yti1 ) exp

  1

ti


−
ρ(t |yt )dt
0

(A.2)

Optimal STDP in Supervised Learning

·


 F



f =2

1343

  f

 


ti
t
f




ρ(ti |yt f ) exp −ρ(t |yt )dt
ρ(t |yt )dt
exp −
i
f −1

tiF
ti


  t

f
ρ(ti |yt f ) exp −ρ(t  |yt )dt  .

=

i

(A.3)

0

f

ti ∈ytAppendix B: Numerical Evaluation of ρ̄(t)Since it is impossible to numerically evaluate the instantaneous ﬁring rateρ̄(t) with the analytical expression given by equation 3.6, we have to do
it in a different way. In fact, there are two ways to evaluate ρ̄(t). Before
going into the details, let us ﬁrst recall that from the law of large numbers,
the instantaneous ﬁring rate is equal to the empirical density of spikes attime t,ρ(t|yt ) yt = Y(t)Y(t) ,

(B.1)
f

where Y(t) = t f ∈yt δ(t − ti ) is one realization of the postsynaptic spikei
train. Thus, the ﬁrst and simpler method based on the right-hand side ofequation B.1 is to build a PSTH by counting spikes in small time bins [t, t +δt] over, say, K = 10,000 repetitions of an experiment. The second, and moreadvanced, method consists in evaluating the left-hand side of equation B.1
by Monte Carlo sampling. Instead of averaging over all possible spike trains
yt , we generate K = 10,000 spike trains by repetition of the same stimulus.
A speciﬁc spike train yt = {ti1 , ti2 , . . . , tiF < t} will automatically appear withappropriate probability given by equation 2.6. The Monte Carlo estimationρ̃(t) of ρ̄(t) can be written as
1 ρ(t|ytm ),PPρ̃(t) =

(B.2)

m=1

where ytm is the mth spike train generated by the stochastic process given
by equation 2.5. Since we use the analytical expression of ρ(t|ytm ), we willcall equation B.2 a semianalytical estimation. Let us note that the semianalytical estimation ρ̃(t) converges more rapidly to the true value ρ̄(t) thanthe empirical estimation based on the PSTH.
In the limit of a Poisson process, η0 = 0, the semianalytical estimation ρ̃(t)given by equation B.2 is equal to the analytical expression of equation 3.6,
since the instantaneous ﬁring rate ρ of a Poisson process is independent of
the ﬁring history yt = {ti1 , ti2 , . . . , tiF < t} of the postsynaptic neuron.

1344

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. GerstnerAppendix C: Deconvolution
C.1 Deconvolution for Spike Pairs. With a learning rule such as equation 3.16, we know the optimal weight change wi j for each synapse, but
we still do not know the corresponding STDP function.
Let us ﬁrst deﬁne the correlation function c k (τ ), k = N(i − 1) + j between
the presynaptic spike train Xij (t) = tpre ∈xij δ(t − t pre ) and the postsynaptic
spike train Yi (t) =c k (τ ) = T
0t des ∈yi δ(t − tdesXij (s)Yi (s + τ )ds,),
k = 1, . . . , NM,(C.1)

T. Since the sum of thewhere we allow a range −T0 ≤ τ ≤ T0 , with T0
pair-based weight change W should be equal to the total adaptation ofweights wk , we can write
 T0
−T0

!

c k (s)W(s)ds = wk

k = 1, . . . , NM.

(C.2)

If we want to express equation C.1 in a matrix form, we need to descretize
time in small bins δt and deﬁne the matrix element,Ck = (+1)δt−T0
δt−T0c k (s)ds.

(C.3)

Now equation C.2 becomes!

w = CW,(C.4)

where w = (w11 , . . . , w1N , w21 , . . . , w MN )T is the vector containing
all the optimal weight change and W is the vector containing the discretized STDP function: W = W(−T0 + δt), for  = 1, . . . , 2T̃0 withT̃0 = T0 /δt.
In order to solve the last matrix equation, we have to compute the inverseof the nonsquare NM × 2T̃0 matrix C, which is known as the Moore-Penrose
inverse (or the pseudo-inverse),
C + = (C T C)−1 C T ,(C.5)

which exists only if (C T C)−1 exists. In fact, the solution given byW = C + w

(C.6)

Optimal STDP in Supervised Learning

1345

minimizes the square distanceD=

1(CW − w)2 .
2

(C.7)

C.2 Temporal Locality Constraint. If we want to impose a constraintof locality, we can add a term in the minimization process of equation C.7
and deﬁne the following,
1
E = D + WT PW,
2

(C.8)

where P is a diagonal matrix that penalizes nonlocal terms. In this article,
we take a quadratic suppression of terms that are nonlocal in time:2P = a δ  − T̃0 .

(C.9)

T̃0 corresponds to the index of the vector W in equations C.4 and C.8
for which t pre − t des = 0. Calculating the gradient of E given by equation
C.8 with respect to W yields∇W E = C T (CW − w) + PW.

(C.10)

By looking at the minimal value of E, that is, ∇W E = 0, we have
W = (C T C + P)−1 C T w.

(C.11)

By setting a = 0, we recover the previous case.
Acknowledgments
This work was supported by the Swiss National Science Foundation(200020-103530/1 and 200020-108093/1). T.T was supported by the Research Fellowships of the Japan Society for the Promotion of Science for
Young Scientists and a Grant-in-Aid for JSPS Fellows.
References
Abarbanel, H., Huerta R., & Rabinovich, M. (2002). Dynamical model of long-term
synaptic plasticity. Proc. Natl. Academy of Sci. USA, 59, 10137–10143.
Atick, J., & Redlich, A. (1990). Towards a theory of early visual processing. NeuralComputation, 4, 559–572.
Barber, D. (2003). Learning in spiking neural assemblies. In S. Becker, S. Thrun, &K. Obermayer (Eds.), Advances in neural information processing systems, 15 (pp.
149–156). Cambridge, MA: MIT Press.

1346

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstner

Barlow, H. B. (1961). Possible principles underlying the transformation of sensory
messages. In W. A. Rosenblith (Ed.), Sensory communication (pp. 217–234). Cambridge, MA: MIT Press.
Bell, A., & Sejnowski, T. (1995). An information maximization approach to blindseparation and blind deconvolution. Neural Computation, 7, 1129–1159.
Bell, A. J., & Parra, L. C. (2005). Maximising sensitivity in a spiking network. In
L. K. Saul, Y. Weiss, & L. Bottou (Eds.), Advances in neural information processing
systems, 17 (pp. 121–128). Cambridge, MA: MIT Press.
Bell, C., Han, V., Sugawara, Y., & Grant, K. (1997). Synaptic plasticity in a cerebellumlike structure depends on temporal order. Nature, 387, 278–281.
Bi, G., & Poo, M. (1998). Synaptic modiﬁcations in cultured hippocampal neurons:
Dependence on spike timing, synaptic strength, and postsynaptic cell type. J.
Neurosci., 18, 10464–10472.
Bi, G., & Poo, M. (1999). Distributed synaptic modiﬁcation in neural networks induced by patterned stimulation. Nature, 401, 792–796.
Bi, G., & Poo, M. (2001). Synaptic modiﬁcation of correlated activity: Hebb’s postulate
revisited. Ann. Rev. Neurosci., 24, 139–166.
Bishop, C. M. (1995). Neural networks for pattern recognition. Oxford: Clarendon Press.
Bohte, S. M., & Mozer, M. C. (2005). Reducing spike train variability: A computational theory of spike-timing dependent plasticity. In L. K. Saul, Y. Weiss, & L.
Bottou (Eds.), Advances in neural information processing systems, 17 (pp. 201–208).
Cambridge, MA: MIT Press.
Brody, C., & Hopﬁeld, J. (2003). Simple networks for spike-timing-based computation, with application to olfactory processing. Neuron, 37, 843–852.
Bugmann, G., Christodoulou, C., & Taylor, J. G. (1997). Role of temporal integrationand ﬂuctuation detection in the highly irregular ﬁring of leaky integrator neuronmodel with partial reset. Neural Computation, 9, 985–1000.
Carr, C. E., & Konishi, M. (1990). A circuit for detection of interaural time differencesin the brain stem of the barn owl. J. Neurosci., 10, 3227–3246.
Chechik, G. (2003). Spike-timing-dependent plasticity and relevant mututal information maximization. Neural Computation, 15, 1481–1510.
Dayan, P., & Häusser, M. (2004). Plasticity kernels and temporal statistics. In S. Thrun,
L. K. Saul, & B. Schölkopf (Eds.), Advances in neural information processing systems,
16. Cambridge, MA: MIT Press.
Feldman, D. (2000). Timing-based LTP and LTD and vertical inputs to layer II/III
pyramidal cells in rat barrel cortex. Neuron, 27, 45–56.
Frégnac, Y., Shulz, D. E., Thorpe, S., & Bienenstock, E. (1988). A cellular analogue ofvisual cortical plasticity. Nature, 333(6171), 367–370.
Frégnac, Y., Shulz, D. E., Thorpe, S., & Bienenstock, E. (1992). Cellular analogs ofvisual cortical epigenesis. I: Plasticity of orientation selectivity. Journal of Neuroscience, 12(4), 1280–1300.
Gerstner, W. (2001). Coding properties of spiking neurons: Reverse- and crosscorrelations. Neural Networks, 14, 599–610.
Gerstner, W., Kempter, R., van Hemmen, J. L., & Wagner, H. (1996). A neuronallearning rule for sub-millisecond temporal coding. Nature, 383, 76–78.
Gerstner, W., & Kistler, W. K. (2002a). Mathematical formulations of Hebbian learning. Biological Cybernetics, 87, 404–415.

Optimal STDP in Supervised Learning

1347

Gerstner, W., & Kistler, W. K. (2002b). Spiking neuron models. Cambridge: Cambridge
University Press.
Gerstner, W., Ritz, R., & van Hemmen, J. L. (1993). Why spikes? Hebbian learningand retrieval of time–resolved excitation patterns. Biol. Cybern., 69, 503–515.
Gütig, R., Aharonov, R., Rotter, S., & Sompolinsky, H. (2003). Learning input correlations through non-linear temporally asymmetry Hebbian plasticity. J. Neuroscience, 23, 3697–3714.
Haykin, S. (1994). Neural networks. Upper Saddle River, NJ: Prentice Hall.
Hopﬁeld, J. J. (1995). Pattern recognition computation using action potential timing
for stimulus representation. Nature, 376, 33–36.
Hopﬁeld, J. J., & Brody, C. D. (2004). Learning rules and network repair in spiketiming-based computation networks. Proc. Natl. Acad. Sci. USA, 101, 337–342.
Izhikevich, E. (2003). Simple model of spiking neurons. IEEE Transactions on Neural
Networks, 14, 1569–1572.
Johansson, R., & Birznieks, I. (2004). First spikes in ensembles of human tactile
afferents code complex spatial ﬁngertip events. Nature Neuroscience, 7, 170–177.
Karmarkar, U., & Buonomano, D. (2002). A model of spike-timing dependent plasticity: One or two coincidence detectors. J. Neurophysiology, 88, 507–513.
Kempter, R., Gerstner, W., & van Hemmen, J. L. (1999). Hebbian learning and spiking
neurons. Phys. Rev. E, 59, 4498–4514.
Kempter, R., Gerstner, W., & van Hemmen, J. L. (2001). Intrinsic stabilization of
output rates by spike-based Hebbian learning. Neural Computation, 13, 2709–2741.
Kistler, W. M., & van Hemmen, J. L. (2000). Modeling synaptic plasticity in conjunction with the timing of pre- and postsynaptic potentials. Neural Comput., 12,385–405.
Legenstein, R., Naeger, C., & Maass, W. (2005). What target functions can be learnt
with spike-timing-dependent plasticity? Neural Computation, 17, 2337–2382.
Magee, J. C., & Johnston, D. (1997). A synaptically controlled associative signal for
Hebbian plastiticy in hippocampal neurons. Science, 275, 209–213.
Markram, H., Lübke, J., Frotscher, M., & Sakmann, B. (1997). Regulation of synaptic
efﬁcacy by coincidence of postysnaptic AP and EPSP. Science, 275, 213–215.
Mehta, M. R., Lee, A. K., & Wilson, M. A. (2002). Role of experience of oscillations intransforming a rate code into a temporal code. Nature, 417, 741–746.
Minsky, M. L., & Papert, S. A. (1969). Perceptrons. Cambridge MA: MIT Press.
Panzeri, S., Peterson, R., Schultz, S., Lebedev, & Diamond, M. (2001). The role of spiketiming in the coding of stimulus location in rat somatosensory cortex. Neuron, 29,
769–777.
Poliakov, A. V., Powers, R. K., & Binder, M. C. (1997). Functional identiﬁcation of the
input-output transforms of motoneurones in the rat and cat. J. Physiology, 504,
401–424.
Rao, R. P. N., & Sejnowski, T. J. (2001). Spike-timing-dependent Hebbian plasticity
as temporal difference learning. Neural Computation, 13, 2221–2237.
Roberts, P. (1999). Computational consequences of temporally asymmetric learning
rules: I. Differential Hebbian learning. J. Computational Neuroscience, 7, 235–246.
Roberts, P., & Bell, C. (2000). Computational consequences of temporally asymmetriclearning rules: II. Sensory image cancellation. Computational Neuroscience, 9, 67–
83.

1348

J.-P. Pﬁster, T. Toyoizumi, D. Barber, and W. Gerstner

Rubin, J., Lee, D. D., & Sompolinsky, H. (2001). Equilibrium properties of temporally
asymmetric Hebbian plasticity. Physical Review Letters, 86, 364–367.
Senn, W., Tsodyks, M., & Markram, H. (2001). An algorithm for modifying neurotransmitter release probability based on pre- and postsynaptic spike timing.
Neural Computation, 13, 35–67.
Seung, S. (2003). Learning in spiking neural networks by reinforcement of stochastic
synaptic transmission. Neuron, 40, 1063–1073.
Shouval, H. Z., Bear, M. F., & Cooper, L. N. (2002). A uniﬁed model of NMDA
receptor dependent bidirectional synaptic plasticity. Proc. Natl. Acad. Sci. USA,99, 10831–10836.
Sjöström, P., Turrigiano, G., & Nelson, S. (2001). Rate, timing, and cooperativity
jointly determine cortical synaptic plasticity. Neuron, 32, 1149–1164.
Song, S., & Abbott, L. (2001). Column and map development and cortical re-mapping
through spike-timing dependent plasticity. Neuron, 32, 339–350.
Song, S., Miller, K., & Abbott, L. (2000). Competitive Hebbian learning throughspike-time-dependent synaptic plasticity. Nature Neuroscience, 3, 919–926.
Thorpe, S., Delorme, A., & Van Rullen, R. (2001). Spike-based strategies for rapid
processing. Neural Networks, 14, 715–725.
Toyoizumi, T., Pﬁster, J.-P., Aihara, K., & Gerstner, W. (2005a). Spike-timing dependent plasticity and mutual information maximization for a spiking neuronmodel. In L. K. Saul, Y. Weiss, and L. Bottou (Eds.), Advances in neural information
processing systems, 17 (pp. 1409–1416). Cambridge, MA: MIT Press.
Toyoizumi, T., Pﬁster, J.-P., Aihara, K., & Gerstner, W. (2005b). Generalized
Bienenstock-Cooper-Munro rule for spiking neurons that maximizes information transmission. Proc. National Academy Sciences (USA), 102, 5239–5244.
Troyer, T. W., & Miller, K. (1997). Physiological gain leads to high ISI variability in a
simple model of a cortical regular spiking cell. Neural Computation, 9, 971–983.
Turrigiano, G., & Nelson, S. (2004). Homeostatic plasticity in the developing nervous
system. Nature Reviews Neuroscience, 5, 97–107.
van Rossum, M. C. W., Bi, G. Q., & Turrigiano, G. G. (2000). Stable Hebbian learningfrom spike timing-dependent plasticity. J. Neuroscience, 20, 8812–8821.
Xie, X., & Seung, S. (2004). Learning in neural networks by reinforcement of irregularspiking. Phys. Rev. E, 69, 041909.
Zhang, L., Tao, H., Holt, C., Harris, W. A., & Poo, M.-M. (1998). A critical window forcooperation and competition among developing retinotectal synapses. Nature,395, 37–44.

Received August 10, 2004; accepted September 8, 2005.