Generalized Bienenstock–Cooper–Munro rule
for spiking neurons that maximizes
information transmission
Taro Toyoizumi*†, Jean-Pascal Pfister*, Kazuyuki Aihara†‡, and Wulfram Gerstner*§*Ecole Polytechnique Fédérale de Lausanne, School of Computer and Communication Sciences and Brain-Mind Institute, 1015 Lausanne, Switzerland;†Department of Complexity Science and Engineering, University of Tokyo, 4-6-1, Komaba, Meguro-ku, Tokyo 153-8505, Japan; and ‡Exploratory Research for
Advanced Technology, Aihara Complexity Modelling Project, Japan Science and Technology Agency, 45-18 Oyama-cho, Shibuya-ku, Tokyo 151-0065, Japan

Maximization of information transmission by a spiking-neuronmodel predicts changes of synaptic connections that depend on
timing of pre- and postsynaptic spikes and on the postsynaptic
membrane potential. Under the assumption of Poisson firing statistics, the synaptic update rule exhibits all of the features of the
Bienenstock–Cooper–Munro rule, in particular, regimes of synaptic
potentiation and depression separated by a sliding threshold.
Moreover, the learning rule is also applicable to the more realistic
case of neuron models with refractoriness, and is sensitive tocorrelations between input spikes, even in the absence of presynaptic rate modulation. The learning rule is found by maximizing themutual information between presynaptic and postsynaptic spiketrains under the constraint that the postsynaptic firing rate stays
close to some target firing rate. An interpretation of the synaptic
update rule in terms of homeostatic synaptic processes and spiketiming-dependent plasticity is discussed.

of information maximization to a spiking-neuron model withrefractoriness. Mutual information is maximized under theconstraint that the postsynaptic firing rate stays as close as
possible to the neuron’s typical target firing rate stabilized byhomeostatic synaptic processes (22). In the special case ofvanishing refractoriness, we find that the optimal update rule has
the two BCM properties, i.e., regimes of potentiation anddepression separated by a sliding threshold. In contrast to the
optimality approach of Intrator and Cooper (23), the sliding
threshold follows automatically from our formulation of the
optimality problem. Moreover, our extension of the BCM rule tospiking neurons with refractoriness shows that synaptic changes
should naturally depend on spike timing, spike frequency, andpostsynaptic potential (PSP), which is in agreement with experimental results.

computational neuroscience 兩 information theory 兩 learning 兩spiking-neuron model 兩 synaptic plasticity

Spiking-Neuron Model. We consider a stochastically spiking-The efficacy of synaptic connections between neurons in thebrain is not fixed, but it varies, depending on the firing
frequency of presynaptic neurons (1, 2), the membrane potentialof the postsynaptic neuron (3), spike timing (4–6), and intracellular parameters such as the calcium concentration; for areview see ref. 7. During the last decades, a large number oftheoretical concepts and mathematical models have emergedthat have helped to understand the functional consequences ofsynaptic modifications, in particular, long-term potentiation(LTP) and long-term depression (LTD) during development,
learning, and memory; for reviews see (8–10). Apart from thework of Hebb (11), one of the most influential theoreticalconcepts has been the Bienenstock–Cooper–Munro (BCM)
model originally developed to account for cortical organization
and receptive field properties during development (12). The
model predicted (i) regimes of both LTD and LTP, dependingon the state of the postsynaptic neuron, and (ii) a sliding
threshold that separates the two regimes. Both predictions i andii have subsequently been confirmed experimentally (2, 13, 14).
In this paper, we construct a bridge between the BCM modeland a seemingly unconnected line of research in theoreticalneuroscience centered around the concept of optimality.
There are indications that several components of neuralsystems show close to optimal performance (15–17). Instead oflooking at a specific implementation of synaptic changes, definedby a rule such as in the BCM model, we therefore ask what would
be the optimal synaptic update rule so as to guarantee that aspiking neuron transmits as much information as possible?
Information theoretic concepts have been used by several researchers because they allow to compare performance of neuralsystems with a fundamental theoretical limit (16, 17), but optimal
synaptic update rules have so far been mostly restricted to a purerate description (18–21). In the following, we apply the concept
www.pnas.org兾cgi兾doi兾10.1073兾pnas.0500495102

Methods and Models
neuron model with refractoriness. For simulations, and also forsome parts of the theory, it is convenient to formulate the modelin discrete time with step size ⌬t, i.e., tk ⫽ k⌬t. However, for theease of interpretation and with respect to a comparison withbiological neurons, it is more practical to turn to continuous timeby taking ⌬t 3 0. The continuous time limit is indicated in thefollowing formulas by a right arrow (3). The postsynaptic
neuron receives input at N synapses. A presynaptic spike train
at synapse j is described in discrete time as a sequence xjk (k ⫽1, . . . , K) of zeros (no spike) and ones (spike). The upper index
k denotes time bin k. Thus, xjk ⫽ 1 indicates that a presynaptic
spike arrived at synapse j at a time tj( f ) with tk⫺1 ⱕ tj( f ) ⬍ tk. Eachpresynaptic spike evokes a PSP of amplitude wj and exponential
time course ␧(t ⫺ tj( f )) with time constant ␶m ⫽ 10 ms. The
membrane potential at time step tk is denoted as u(tk) andcalculated as the total PSP冘冘N

u共tk兲 ⫽ ur ⫹kwj␧共tk ⫺ tn兲xjn 3 ur ⫹j⫽1 n⫽1冘冘jwj␧共t ⫺ tj共f兲兲,f[1]
where ur ⫽ ⫺70 mV is the resting potential. The probability ␳kof firing in time step k is a function of the membrane potentialu and the refractory state R of the neuron,␳k ⫽ 1 ⫺ exp关⫺g共u共t k兲兲R共t k兲⌬t兴 ⬇ g共u共t k兲兲R共t k兲⌬t ,[2]Freely available online through the PNAS open access option.
Abbreviations: BCM, Bienenstock–Cooper–Munro; LTP, long-term potentiation; LTD, longterm depression; PSP, postsynaptic potential.
§To whom correspondence should be addressed. E-mail: wulfram.gerstner@epfl.ch.

© 2005 by The National Academy of Sciences of the USA

PNAS 兩 April 5, 2005 兩 vol. 102 兩 no. 14 兩 5239 –5244

NEUROSCIENCE

Communicated by Leon N. Cooper, Brown University, Providence, RI, January 25, 2005 (received for review November 5, 2004)

where ⌬t is the time step and g is a smooth increasing functionof u. Thus, the larger the membrane potential, the higher thefiring probability. For ⌬t 3 0, we may think of g(u)R(t) as theinstantaneous firing rate, or hazard of firing, given knowledge
about the previous firing history. We focus on nonadapting
neurons where the refractoriness R depends only on the timingof the last postsynaptic spike, but the model can be easily
generalized to include a dependence on earlier spikes as well.
More specifically, we take for the simulations
R共t兲 ⫽

共t ⫺ t̂ ⫺ ␶abs兲 2⌰共t ⫺ t̂ ⫺ ␶ abs兲,2
␶ refr ⫹ 共t ⫺ t̂ ⫺ ␶ abs兲 2[3]

where t̂ denotes the last firing time of the postsynaptic neuron,␶abs ⫽ 3 ms is the absolute refractory time, and ␶refr ⫽ 10 ms is
a parameter characterizing the duration of relative refractoriness. The Heaviside function ⌰(x) takes a value of one forpositive arguments and vanishes otherwise. With a function R(t)such as in Eq. 3 that depends only on the most recent postsynaptic spike, the above neuron model has renewal properties andcan be mapped onto a spike-response model with escape noise(9). Except for Fig. 2, we take throughout the paper g(u) ⫽r0 log{1 ⫹ exp[(u ⫺ u0)兾⌬u]} with u0 ⫽ ⫺65 mV, ⌬u ⫽ 2 mV,
and r0 ⫽ 11 Hz. This set of parameters corresponds to in vivo
conditions with a spontaneous firing rate of ⬇1 Hz. The function
g(u) and the typical firing behavior of the neuron model are
shown in Fig. 1A. For Fig. 2, we consider the case ␶abs ⫽ ␶refr ⫽0 and an instantaneous rate g2(u) ⫽ {10 ms ⫹ [1兾g(u)]}⫺1 withg(u) as above i.e., the neuron model exhibits no refractorinessand is defined by an inhomogeneous Poisson process withmaximum rate of 100 Hz; compare Fig. 1 A. Integration of all
equations is performed in MATLAB on a standard personal
computer by using a time step ⌬t ⫽ 1 ms.
Spike Trains. The output of the postsynaptic neuron at time step

tk is denoted by a variable yk ⫽ 1 if a postsynaptic spike occurred
between tk⫺1 and tk and 0 otherwise. A specific output spike train
up to time bin k is denoted by uppercase letters Yk ⫽ {y1, y2, . . . ,yk}. Because spikes are generated by a random process, we
distinguish the random variable Yk by using boldface charactersfrom a specific realization Yk. Note that the lowercase variableyk refers always to a specific time bin, whereas the uppercase
variable Yk refers to a whole spike train. Similar remarks hold forthe input: X is the random variable characterizing the input at all
synapse 1 ⱕ j ⱕ N; Xk is a specific realization of all input spiketrains up to time tk; Xjk ⫽ {xj1, xj2, . . . , xjk} a specific realizationof an input spike train at synapse j, and xjk, its value in time bin
k. For given presynaptic spike trains Xk and postsynaptic spikehistory Yk⫺1, the probability of emitting a postsynaptic spike is
described by the following binary distribution:
P共yk兩Yk⫺1, Xk兲 ⫽ 共␳k兲y 共1 ⫺ ␳k兲共1⫺y 兲
k

k

[4]

since yk 僆 {0, 1}. Analogously, we find the marginal probabilityof yk given Yk⫺1
P共yk兩Yk⫺1兲 ⫽ 共␳៮ k兲y 共1 ⫺ ␳៮ k兲共1⫺y 兲
k

k

[5]

写
k

P共Y 兩X 兲 ⫽k

P共yl兩Yl⫺1, Xl兲,l⫽15240 兩 www.pnas.org兾cgi兾doi兾10.1073兾pnas.0500495102

and an analogous formula for the marginal probability distribution of P(Yk). With Eqs. 4 and 6, we have an expression for the
probabilistic relation between an output spike train and an
ensemble of input spike trains.
Mutual Information Optimization. Transmission of information
between an ensemble of presynaptic spike trains XK of total
duration K⌬t and the output train YK of the postsynaptic neuron
can be quantified by the mutual information (24)

I共YK; XK兲 ⫽冘

Y K,X K

where ␳៮ k ⫽ 具␳k典Xk兩Yk⫺1 and 具䡠典Xk兩Yk⫺1 ⫽ 兺Xk 䡠 P(Xk 兩Yk⫺1).
From probability calculus, we obtain the conditional probability of the output spike train Yk given the presynaptic spiketrains Xk,
kFig. 1. Neuron model. (A Left) Output rate ␯ post of the model neuron (solidline, spiking neuron model used in Figs. 3–5; green dashed line, Poisson model
used in Fig. 2) as a function of presynaptic spike arrival rate at n ⫽ 100 synapses.
All synapses have the same efficacy wj ⫽ 0.5, and are stimulated by independent Poisson trains at the same rate ␯. (Center) Interspike interval distribution
PISI of the spiking neuron model during firing at 10 (blue line), 20 (green line),or 30 Hz (red line). Firing is impossible during the absolute refractory time of␶abs ⫽ 3 ms. (Right) The function g(u) used to generated action potentials (seeMethods and Models for details). (B) From the first row to the fourth row (withthe first row being at the top): The measure Cj that is sensitive to correlations
between the state of the postsynaptic neuron and presynaptic spike arrival atsynapse j, the PSPs caused by spike arrival at the same synapse j, the membranepotential u, and the postsynaptic factor Bpost of Eq. 14 as a function of time.
During postsynaptic action potentials, the postsynaptic factor Bpost has
marked peaks. Their amplitude and sign depend on the membrane potential
at the moment of action potential firing. The coincidence measure Cj exhibits
significant changes only during the duration of PSPs at synapse j.

[6]

P共Y K, X K兲log

P共Y K兩X K兲
.
P共Y K兲

[7]

While it is easier to transmit information if the postsynaptic
neuron increases its firing rate, firing at high rates is costly fromthe point of view of energy consumption and also difficult toimplement by the cells biophysical machinery. We therefore
optimize information transmission under the condition that thefiring statistics P(YK) of the postsynaptic neuron stays as close as
possible to a target distribution P̃(YK). With a parameter ␥ (setto ␥ ⫽ 1 for the simulations), the quantity we maximize istherefore
Toyoizumi et al.

train of the postsynaptic neuron (see Supporting Text and Fig. 6
which are published as supporting information on the PNAS website). The term
Fk ⫽ log

␳k
P共y k兩Y k⫺1, X k兲
1 ⫺ ␳k
k
k
⫽
y
log
⫹
共1
⫺
y
兲log
P共y k兩Y k⫺1兲␳៮ k1 ⫺ ␳៮ k[12]

␳k at time step k[8]where D(P(YK)兩兩P̃(YK)) ⫽ 兺YKP(YK)log(P(Y K)兾P̃(Y K)) denotes
the Kullback–Leibler divergence (24). The target distribution isthat of a neuron with constant instantaneous rate g̃ (set to g̃ ⫽30 Hz throughout the paper except for Fig. 2) modulated by therefractory variable R(t), i.e., that of a renewal process.
The main idea of our approach is as follows. We assume thatsynaptic efficacies wj can change within some bounds 0 ⱕ wj ⱕ
wmax so as to maximize information transmission under theconstraint of a fixed target firing rate. To derive the optimal ruleof synaptic update, we calculate the gradient of Eq. 8. Applying
the chain rule of information theory (24) to both the mutual
information I and the Kullback–Leibler divergence D, we canKwrite L ⫽ 兺 k⫽1⌬Lk where

具P共y k兩Y k⫺1, X k兲P共y k兩Y k⫺1兲⌬L ⫽ log⫺ ␥ logkk⫺1
P共y 兩Y 兲P̃共y k兩Y k⫺1兲

典

,
Yk,Xk

[9]
Assuming slow changes ofwith 具䡠典Xk,Yk ⫽ 兺Xk,Yk 䡠synaptic weights, we apply a gradient ascent algorithm to maximize the objective function and change the synaptic efficacy wj
at each time step by ⌬w jk ⫽ ␣ (⭸⌬Lk兾⭸w j) with an appropriate
learning rate ␣. Evaluation of the gradient (see Supporting Text,
which is published as supporting information on the PNAS website) yields
P(Xk,

Yk).

⌬wjk ⫽ ␣具Cjk共Fk ⫺ ␥Gk兲典Xk,Yk,[10]

with three functions Cj, F, and G, described in the following.
First, with ␳⬘ denoting the derivative of ␳ with respect to u, thequantity

冘 冘
k

Cjk ⫽

l

l⫽k⫺ka n⫽1

␧共tl ⫺ tn兲xjn冋␳⬘l l 1 ⫺ yl ly ⫺␳␳l
1 ⫺ ␳l册[11]

is a measure that counts coincidences between postsynaptic
spikes (yl ⫽ 1) and the time course of PSPs generated bypresynaptic spikes (xjn ⫽ 1) at synapse j, normalized to an
expected value 具Cjk典Yk兩XK ⫽ 0. The time span ka of the coincidence
window is given by the width of the autocorrelation of the spikeToyoizumi et al.

冋

冉 冊册

g共u共t兲兲 g̃
g៮ 共t兲
g៮ 共t兲

␥⫺R共t兲关g共u共t兲兲 ⫺ 共1 ⫹ ␥ 兲g៮ 共t兲 ⫹ ␥ g̃兴,L ⫽ I共YK; XK兲 ⫺ ␥ D共P共Y K兲储P̃共Y K兲兲,k

Bpost共t兲 ⫽ ␦ 共t ⫺ t̂ ⫺ ␦ 兲log

[13]

where t̂ is the firing time of the last postsynaptic spike. The delay␦ in the Dirac-␦ function reflects the order of updates in a single
time step of the numerical implementation, i.e., we first update
the membrane potential, then the last firing time t̂, then thefactors Cj and B, and finally the synaptic efficacy wj; for themathematical theory, we take ␦ 3 0. The rate g៮ (t) ⫽ 具g(u(t))典X兩Ydenotes an expectation over the input distribution given the
recent firing history of the postsynaptic neuron. For a numerical
implementation, it is convenient to estimate the expected rateg(t) by a running average with exponential time window (a time
constant of 10 s). Similarly, we replace the rectangular coincidence count window in Eq. 11 by an exponential one (with a time
constant of 1 s; see Supporting Text).
Results
We analyzed information transmission for a model neuron thatreceives input from 100 presynaptic neurons. A presynaptic spikethat arrives at time tj( f ) at synapse j evokes an excitatory PSP oftime course ␧(t ⫺ tj( f )). The amplitude wj of the postsynaptic
response is taken as a measure of synaptic efficacy and subject
to synaptic dynamics. Firing of the postsynaptic neuron is more
likely if the total PSP u(t) ⫽ ur ⫹ 兺j 兺f wj␧(t ⫺ tj( f )) is large;
however, because of refractoriness firing is suppressed after apostsynaptic action potential at time t̂ by a factor R, which
depends on the time since the last postsynaptic spike; see
methods for details.
Maximizing the mutual information between several presynaptic spike trains and the output of the postsynaptic neuron canbe achieved by a synaptic update rule that depends on thepresynaptic spike arrival time tj( f ), the postsynaptic membranepotential u, and the last postsynaptic firing time t̂. More precisely,
the synaptic update rule can be written as
dwj⫽ ␣Cj共t兲Bpost共t ⫺ ␦ 兲
dt

[14]

where Cj is a measure sensitive to correlations between pre- andpostsynaptic activity, and B post is a variable that characterizes thestate of the postsynaptic neuron (see Methods and Models). ␣ isa small learning parameter. We note that in standard formulations of Hebbian learning, changes of synaptic efficacies are
driven by correlations between pre- and postsynaptic neurons,
PNAS 兩 April 5, 2005 兩 vol. 102 兩 no. 14 兩 5241

NEUROSCIENCE

Fig. 2. Relation to BCM rule. (A) The function ␾(␯ post) of the BCM learning
rule Eq. 16 derived from our model under assumption of Poisson firing
statistics of the postsynaptic neuron. A value of ␾(␯ post) ⬎ 0 for a given
postsynaptic rate ␯ post means that synapses are potentiated when stimulated
presynaptically. The transition from depression to potentiation occurs at avalue ␪ that depends on the average firing rate ␯៮ post of the postsynaptic
neuron (blue ␯៮ post ⫽ 10 Hz; green ␯៮ post ⫽ 20 Hz; red ␯៮ post ⫽ 30Hz). B. The
threshold ␪ as a function of ␯៮ post for different choices of the parameter ␥, i.e.,␥ ⫽ 0.5 (purple); ␥ ⫽ 1 (black); ␥ ⫽ 2 (orange).

compares the instantaneous firing probability
with the average probability, ␳៮ k and, analogously, the term Gk ⫽log(P(y k兩Y k⫺1)兾P̃(y k兩Y k⫺1)) compares the average probability
with the target value ˜␳k ⫽ g̃R(tk)⌬t. We note that both F and G
are functions of the postsynaptic variables only. We therefore
introduce a postsynaptic factor Bpost by the definition Bpost(tk) ⫽[Fk ⫺ ␥Gk]兾⌬t and take the limit ⌬t 3 0. Under the assumptionof small learning rate ␣ (i.e., ␣ ⫽ 10⫺4 in our simulations), theexpectations 具典X,Y in Eq. 10 can be approximated by averaging
over a single long trial that allows us to define an on-line rule(dw j(t)兾dt) ⫽ ␣Cj(t)Bpost(t ⫺ ␦) with a postsynaptic factor

similar to the function Cj(t). The above update rule, however,augments these correlations by a further postsynaptic factor
Bpost.
This postsynaptic factor Bpost depends on the firing time of thepostsynaptic neuron, the refractory state of the neuron, its
membrane potential u by means of the instantaneous firing
intensity g(u), and on its past firing history by means of g៮ (t). The
postsynaptic factor can be decomposed into two terms: the firstterm compares the instantaneous firing intensity g(u) with its
running average g៮ (t) and the second term compares the runningaverage with a target rate g̃. Thus, the first term of Bpost measures
momentary significance of the postsynaptic state, whereas its
second term accounts for homeostatic processes; see Methods
and Models for details.
The variable Cj measures correlations between the postsynaptic neuron and its presynaptic input at synapse j
C j共t ⫺ ␦ 兲
dCj共t兲
⫽⫺
dt
␶C
⫹

冘

␧共t ⫺ t j共f兲兲S共t兲关 ␦ 共t ⫺ t̂ ⫺ ␦ 兲 ⫺ g共u共t兲兲R共t兲兴,f[15]
with time constant ␶C ⫽ 1s. Here g(u(t))R(t) is the instantaneous
firing rate modulated by the refractory function R(t), and S(t) ⫽g⬘(u(t))兾g(u(t)) is the sensitivity (the prime denoting the derivatewith respect to u) of the neuron to a change of its membranepotential. The term with the Dirac ␦-function ␦(t ⫺ t̂ ⫺ ␦) induces
a positive jump of Cj immediately (with short delay ␦) after eachpostsynaptic spike. Between postsynaptic spikes, Cj evolves
continuously. Significant changes of Cj are conditioned on thepresence of a PSP ␧(t ⫺ tj( f )) caused by spike arrival at synapse
j. In the absence of presynaptic input, the correlation estimate
decays with time constant ␶C back to 0.
Both the correlation term Cj and the postsynaptic factor Bpost
can be estimated online (Fig. 1) and use only information that
could, directly or indirectly, be available at the site of thesynapse: information about postsynaptic spike firing could be
conveyed by backpropagating action potentials; timing of presynaptic spike arrival is transmitted by neurotransmitter receptors; and the total PSP can be estimated, although not perfectly,
from the local potential at the synapse. The direction of change
of a synapse is determined by a subtle interplay between thecorrelation term Cj and the postsynaptic factor Bpost, which canboth be positive or negative.
To elucidate the balance between potentiation and depressionof synapses, we first considered a simplified neuron model
without refractoriness and firing rate ␯ post ⫽ g2(u). In this specialcase, the synaptic update rule Eq. 10 can be rewritten in thesimpler formdw ⫽ ␣␯j␾共␯ post, ␪ 兲,dt j

冉␯៮ post共t兲
g̃

冊

␥

,

[17]

where g̃ ⫽ 20 Hz denotes a target value of the postsynaptic rate
implemented by homeostatic processes (22) and ␯៮ post(t) is arunning average of the postsynaptic rate. The function ␾ in Eq.
16 shown in Fig. 2 is characteristic for the BCM learning rule(12). Our approach by information maximization predicts a
specific form of this function that can be plotted either as afunction of the postsynaptic firing rate ␯ post or as a function ofpost), in close agreement with experiments
the total PSP u ⫽ g⫺1
2 (␯
(2, 13). Moreover, because information maximization was performed under the constraint of a fixed target firing rate, ourapproach yields automatically a sliding threshold of the form
postulated in ref. 12, but on different grounds. Thus for neurons
without refractoriness, i.e., a pure rate model, our update rule
for synaptic plasticity reduces exactly to the BCM rule.
An application of rate models to stimulation paradigms thatvary on a time scale of tens of milliseconds or less has often beenquestioned because an interpretation of rate is seen as problematic. Our synaptic update rule is based on a spiking-neuronmodel that includes refractoriness and that captures propertiesof much more detailed neuron models well (25). The learning
rule for spiking neurons has a couple of remarkable properties
that we explore now.
First, we consider a pattern discrimination task in a rate
coding paradigm. Patterns are defined by the firing rate ␯pre of
25 presynaptic neurons (␯pre ⫽ 2, 13, 25, and 40 Hz for patterns,
1–4, respectively) modeled as independent Poisson spike trains.
The remaining 75 synapses received uncorrelated Poisson input
at a constant rate of 20 Hz. Each second, a pattern was chosenstochastically and applied during 1 sec. Those synapses thatreceived pattern-dependent input developed strong efficacies
close to the maximal value wmax ⫽ 1, whereas most of the other75 synapses developed weaker ones; compare Fig. 3. However,[16]

where vj is the instantaneous firing rate of the presynaptic neuron
as estimated from the amplitude of the PSP. vj ⫽ a 兺f ␧(t ⫺ tj( f ))generated at synapse j [If the potential is measured in millivoltsand time in milliseconds, then the proportionality constant a hasunits 1兾(mV ms)]. ␾(␯ post, ␪) ⫽ f(␯ post)log(␯ post兾␪) is a functionthat depends on the instantaneous postsynaptic firing rate␯ post ⫽ g2(u) and a parameter ␪. The function f is proportional
to the derivative of g2, i.e., f(␯ post) ⫽ g⬘2兾a evaluated at u ⫽post). The parameter ␪ denotes the transition from a regimeg⫺12 (␯of potentiation to that of depression. It depends on the recent
firing history of the neuron and is given by5242 兩 www.pnas.org兾cgi兾doi兾10.1073兾pnas.0500495102

␪共t兲 ⫽ ␯៮ post共t兲

Fig. 3. Pattern discrimination. The first 25 synapses 1 ⱕ j ⱕ 25 are stimulated
by Poisson input with a rate ␯ pre ⫽ 2, 13, 25, and 40 Hz that changes eachsecond. The remaining 75 synapses receive Poisson input at a constant rate of
20 Hz. (A Upper) Evolution of all synaptic weights as a function of time (red;strong synapses, wj ⬇ 1; blue: depressed synapses, wj ⬇ 0). All synapses are
initialized at the same value wj ⫽ 0.1. (Lower) The evolution of the average
efficacy of the 25 synapses that receive pattern-dependent input (red line) andthat of the 75 other synapses (blue). Typical examples of individual traces(synapse 1: black and synapses 30: green) are given by the dashed lines. (BUpper) Evolution of the average mutual information I per bin (blue line andleft scale) and of the average Kullback–Leibler distance D per bin as a functionof time. Averages are calculated over segments of 1 min. (Lower) Output rate(spike count during 1 sec) as a function of pattern index before (blue bars) andafter (red bars) learning.

Toyoizumi et al.

Fig. 4. Rate modulation. (A) Distribution of synaptic efficacies of nine
postsynaptic neurons after 60 min of stimulation with identical inputs for all
neurons. Synapses 1 ⱕ j ⱕ 40 (red symbols) received Poisson input with
common rate modulation; the input at synapses 41 ⱕ j ⱕ 80 (blue symbols) was
also rate-modulated but phase-shifted; and the input at the remaining 20
synapses was uncorrelated (green symbols). Four postsynaptic neurons (numbers 1, 3, 5, and 8) develop a spontaneous specialization for the first group of
modulated input (red symbols close to the maximum efficacy of one) and five
(numbers 2, 4, 6, 7, and 9) specialized for the second group. (B) Modulation of
the output rate of the nine postsynaptic neurons before (Left) and after(Right) learning. Red兾blue bars, neurons responding to the first兾second groupof input; Red兾blue lines, modulation of the input of groups 1 and 2.

some of the weakly driven synapses also spontaneously increased
their efficacy. This was necessary for the neuron to achieve amean postsynaptic firing rate close to the target firing rate.
Despite the fact that the mean firing rate approaches on a timescale of tens of seconds the target rate, the spike count in each
1-sec segment is strongly modulated by the input pattern and canbe used to distinguish the patterns (Fig. 3B) with a misclassification error of ⬍20%. (See Supporting Text and Fig. 7, which is
published as supporting information on the PNAS web site.)
In a second set of simulation experiments, we again studied
rate coding, but considered two groups of input defined byperiodic modulation of the presynaptic firing rates. More precisely, 40 neurons received input spike trains generated by aninhomogeneous Poisson process with common rate modulationv(t) ⫽ v0 ⫹ A sin (2␲ t兾T) with amplitude A ⫽ 10 Hz and periodT ⫽ 100 ms. Another group of 40 neurons received modulated
input of the same amplitude and period, however, with a phaseshift ␲. The remaining 20 synapses received Poisson input at a
fixed rate v0 ⫽ 20 Hz. All 100 inputs projected onto nine
postsynaptic neurons. Synaptic weights were initialized randomly(between 0.10 and 0.12) and evolved according to the update
rule Eq. 14. Of the nine postsynaptic neurons, four developed
strong connections to the first group of correlated inputs,
whereas five developed strong connections to the second group
Toyoizumi et al.

of correlated inputs (Fig. 4A). The rate of the output neurons
reflects the modulation of their respective inputs (Fig. 4B). To
summarize, the synaptic update rule derived from the principle
of information maximization drives neurons to spontaneously
detect and specialize for groups of coherent inputs. Just as in thestandard BCM rule, several output neurons (with different
specialization) are needed to account for the different featuresof the input.
Whereas the two preceding paradigms focused on rate modulation, we now show that, even if all presynaptic neurons fire at
the same mean rate, the presence of weak spike-spike correlations in the input is sufficient to bias the synaptic selection
mechanism; compare Fig. 5. All synapses received spike input at
the same rate of 20 Hz, but the spike trains of 50 synapses showed
weak correlations (c ⫽ 0.1), whereas the remaining 50 synapses
received uncorrelated input. Most of the 50 synapses thatreceived a weakly correlated input increased their weights underapplication of the synaptic update rule (Fig. 5). If at a later stage
the group of correlated inputs changes, the newly included
synapses will be strengthened as well, whereas those that are nolonger correlated decay. Moreover, the final distribution ofsynaptic efficacies persists, even if the input is switched torandom spike arrival. Thus the synaptic update rule is sensitive
to spike–spike correlations on a millisecond scale, which would
be difficult to account for in a pure rate model. A biophysicalsignature of spike–spike correlations are systematic and large
fluctuations of the membrane potential. If several postsynaptic
neurons receive the same input, their outputs are again correlated and generate large fluctuations in the membrane potential
PNAS 兩 April 5, 2005 兩 vol. 102 兩 no. 14 兩 5243

NEUROSCIENCE

Fig. 5. Spike–spike correlations. The n ⫽ 100 synapses have been separated
into four groups of 25 neurons each (group A, 1 ⱕ j ⱕ 25; group B, 26 ⱕ j ⱕ 50;
group C, 51 ⱕ j ⱕ 75; group D, 76 ⱕ j ⱕ 100). All synapses were stimulated at
the same rate of 20 Hz. However, during the first 15 min of simulated time,neurons in groups C and D were uncorrelated, whereas the spike trains of theremaining 50 neurons (groups A and B) had correlations of amplitude c ⫽ 0.1,
i.e., 10% of the spike arrival times were identical between each pair ofsynapses. After 15 min, correlations changed so that group A became correlated with C, whereas B and D were uncorrelated. After 45 min of simulated
time, correlations stopped, but stimulation continued at the same rate. (AUpper) Evolution of all 100 weights (red, potentiated; blue, depressed).
(Lower) Average mutual information per bin as a function of time. In theabsence of correlations (t ⬎ 45 min), mutual information is lower than before,
but the distribution of synaptic weights remains stable. (B) Nine postsynaptic
neurons 1 ⱕ i ⱕ 9 with membrane potential ui(t) are stimulated as discussed
in A and project to a readout unit with potential h(t) ⫽ 兺i9⫽ 1 兺m␧(t ⫺ tim)
where the sum runs over all output spikes m of all nine neurons. Mean
membrane potentials are u៮ and h៮ , respectively. The fluctuations ␴u ⫽ 具 (ui(t) ⫺u៮ )2典0.5 of the PSPs (blue line, top graph) and those of the readout potentials(␴h ⫽ 具(h(t) ⫺ h៮ )2典0.5, green line) are correlated (Lower) with the mutual
information and can serve as neuronal signal.

of a readout neuron further down the processing chain (Fig. 5B).
Thus, information that is potentially encoded in millisecond
correlations in the input can be detected, enhanced, transmitted,
and read out by other neurons.

ingredients of the spiking neuron model considered here were (i)PSPs generated by presynaptic spike arrival, (ii) a heuristic
spiking probability that depends on the total PSP, and (iii) a
phenomenological account of absolute or relative refractoriness.
The synaptic update rule depends on all three of these quantities.
While we do not imply that synaptic potentiation and depressionof real neurons are implemented the way it is suggested by our
update rule, the rule shows nevertheless some interesting
features.
First, in contrast to pure Hebbian correlation driven learning,
the update rule uses a correlation term modulated by anadditional postsynaptic factor. Thus, presynaptic stimulation is
combined with a highly nonlinear function of the postsynaptic
state to determine the direction and amplitude of synaptic
changes. The essence of the BCM rule (presynaptic gating
combined with nonlinear postsynaptic term) is hence translated
into a spike-based formulation.
Second, the spike-based formulation of a synaptic update ruleshould allow a connection to spike-timing-dependent plasticity(5, 6) and allow its interpretation in terms of optimal informationtransmission (27–29). Given the highly nonlinear involvement of
postsynaptic spike times and PSP in the optimal synaptic update
rule, a simple interpretation in terms of pairs or pre- andpostsynaptic spikes as in many standard models of synaptic
plasticity (30, 31) can only capture a small portion of synaptic
plasticity phenomena. The optimal learning rule suggests that
nonlinear phenomena (32–35) are potentially highly relevant.

Discussion
The synaptic update rule discussed in this paper relies on the
maximization of the mutual information between an ensemble ofpresynaptic spike trains and the output of the postsynaptic
neuron. As in all optimization approaches, optimization has tobe performed under some constraint. Because information
scales with the postsynaptic firing rate, but high firing rates
cannot be sustained by the biophysical machinery of the cell overlong times, we imposed that, on average, the postsynaptic firing
rate should stay close to a desired firing rate. This idea isconsistent with the widespread finding of homeostatic processes
that tend to push a neuron always back into its preferred firing
state (22). The implementation of this idea in our formalism gavenaturally rise to a control mechanism that corresponds exactly tothe sliding threshold in the original BCM model (12). Whilederivations of the BCM model from optimality concepts (23, 19)
or statistical approaches (26) are not new, our approach gives
another perspective on the concept of a sliding threshold.
Our derivation extends the BCM model, which was originallydesigned for rate models of neuronal activity to the case of
spiking-neuron models with refractoriness. Spiking-neuronmodels of the integrate-and-fire type, such as the one presented
here, can be used to account for a broad spectrum of neuronalfiring behavior, including the role of spike–spike correlations,interspike interval distributions, coefficient of variations, andeven timing of single spikes; for a review, see ref. 9. The essential

This work was supported by the Japan Society for the Promotion of
Science and a Grant-in-Aid for Japan Society for the Promotion of
Science Fellows (to T.T.) and by the Swiss National Science Foundation(to J.-P.P.).

1. Bliss, T. V. P. & Gardner-Medwin, A. R. (1973) J. Physiol. 232, 357–374.
2. Dudek, S. M. & Bear, M. F. (1993) J. Neurosci. 13, 2910–2918.
3. Artola, A., Bröcher, S. & Singer, W. (1990) Nature 347, 69–72.
4. Levy, W. B. & Stewart, D. (1983) Neuroscience 8, 791–797.
5. Markram, H., Lübke, J., Frotscher, M. & Sakmann, B. (1997) Science 275,
213–215.
6. Bi, G.-Q. & Poo, M.-M. (2001) Annu. Rev. Neurosci. 24, 139–166.
7. Malenka, R. C. & Nicoll, R. A. (1999) Science 285, 1870–1874.
8. Dayan, P. & Abbott, L. F. (2001) Theoretical Neuroscience (MIT Press,
Cambridge, MA).
9. Gerstner, W. & Kistler, W. K. (2002) Spiking Neuron Models (Cambridge Univ.
Press, Cambridge, U.K.).
10. Cooper, L. N., Intrator, N., Blais, B. S. & Shouval, H. Z. (2004) Theory ofCortical Plasticity (World Scientific, Singapore).
11. Hebb, D. O. (1949) The Organization of Behavior (Wiley, New York).
12. Bienenstock, E. L., Cooper, L. N. & Munro, P. W. (1982) J. Neurosci. 2, 32–48;
13. Artola, A. & Singer, W. (1993) Trends Neurosci. 16, 480–487.
14. Kirkwood, A., Rioult, M. G. & Bear, M. F. (1996) Nature 381, 526–528.
15. Barlow, H. B. & Levick, W. R. (1969) J. Physiol. (London) 200, 1–24.
16. Laughlin, S. (1981) Z. Naturforsch. 36, 910–912.
17. Fairhall, A. L., Lewen, G. D, Bialek, W. & vanSteveninck, R. R. D. (2001)
Nature 412, 787–792.
18. Linsker, R. (1989) Neural Comput. 1, 402–411.
19. Nadal, J.-P. & Parga, N. (1997) Neural Comput. 9, 1421–1456.
20. Bell, A. J. & Sejnowski, T. J. (1995) Neural Comput. 7, 1129–1159.

21. Chechik, G. (2003) Neural Comput. 15, 1481–1510.
22. Turrigiano, G. G. & Nelson, S. B. (2004) Nat. Rev. Neurosci. 5, 97–107.
23. Intrator, N. & Cooper, L. N. (1992) Neural Networks 5, 3–17.
24. Cover, T. M. & Thomas, J. A. (1991) Elements of Information Theory (Wiley,
New York).
25. Jolivet, R., Lewis, T. J. & Gerstner, W. (2004) J. Neurophsyiol. 92,
959 –976.
26. Beggs, J. M. (2000) Neural Comput. 13, 87–111.
27. Toyoizumi, T., Pfister, J.-P., Aihara, K. & Gerstner, W. (2005) in Advances in
Neural Information Processing Systems, eds. Saul, L. K., Weiss, Y. & Bottou, L.
(MIT Press, Cambridge, MA), Vol. 17, pp. 1409–1416.
28. Bell, A. J. & Parra, L. C. (2005) in Advances in Neural Information Processing
Systems, eds. Saul, L. K., Weiss, Y. & Bottou, L. (MIT Press, Cambridge, MA),
Vol. 17, pp. 121–128.
29. Bohte, S. M. & Mozer, M. C. (2005) in Advances in Neural Information
Processing Systems, eds. Saul, L. K., Weiss, Y. & Bottou, L. (MIT Press,
Cambridge, MA), Vol. 17, pp. 201–208.
30. Gerstner, W., Kempter, R., vanHemmen, J. L. & Wagner, H. (1996) Nature 383,
76–78.
31. Song, S., Miller, K. D. & Abbott, L. F. (2000) Nat. Neurosci. 3, 919–926.
32. Froemke, R. & Dan, Y. (2002) Nature 416, 433–438.
33. Bi, G.-Q. & Wang, H.-X. (2002) Physiol. Behav. 77, 551–555.
34. Senn, W. (2002) Biol. Cybern. 87, 344–355.
35. Sjöström, P. J., Turrigiano, G. G. & Nelson, S. B. (2001) Neuron 32,1149 –1164.

5244 兩 www.pnas.org兾cgi兾doi兾10.1073兾pnas.0500495102

Toyoizumi et al.