David Sussillo, Taro Toyoizumi and Wolfgang Maass
J Neurophysiol 97:4079-4095, 2007. First published Apr 4, 2007; doi:10.1152/jn.01357.2006
You might find this additional information useful...
Supplemental material for this article can be found at:
http://jn.physiology.org/cgi/content/full/01357.2006/DC1
This article cites 34 articles, 14 of which you can access free at:
http://jn.physiology.org/cgi/content/full/97/6/4079#BIBL
Updated information and services including high-resolution figures, can be found at:
http://jn.physiology.org/cgi/content/full/97/6/4079
Additional material and information about Journal of Neurophysiology can be found at:
http://www.the-aps.org/publications/jn

This information is current as of June 12, 2007 .
Downloaded from jn.physiology.org on June 12, 2007

Journal of Neurophysiology publishes original articles on the function of the nervous system. It is published 12 times a year
(monthly) by the American Physiological Society, 9650 Rockville Pike, Bethesda MD 20814-3991. Copyright © 2005 by theAmerican Physiological Society. ISSN: 0022-3077, ESSN: 1522-1598. Visit our website at http://www.the-aps.org/.

J Neurophysiol 97: 4079 – 4095, 2007.
First published April 4, 2007; doi:10.1152/jn.01357.2006.

Self-Tuning of Neural Circuits Through Short-Term Synaptic Plasticity
David Sussillo,1 Taro Toyoizumi,1 and Wolfgang Maass2
1

Center for Theoretical Neuroscience in the Center for Neurobiology and Behavior, Columbia University, New York, New York;and 2Institute for Theoretical Computer Science, Technische Universitaet Graz, Graz, Austria
Submitted 24 December 2006; accepted in final form 28 March 2007

INTRODUCTION

Cortical neurons fire not only in response to sensory inputs,
but also spontaneously, as part of an ongoing network activityof the neocortex. According to Steriade (2001), the rates of
spontaneous firing recorded extracellularly from motor andassociation areas of awake monkeys and cats are typically inthe range of 10 –15 Hz. Baddeley et al. (1997) found that cells
in area IT (recorded extracellularly in awake macaque) fire at
a rate of 14 ⫾ 8.3 Hz during blank-screen viewing. Intracellular recordings from regular-spiking neurons (presumed to be
pyramidal cells) in motor, association, primary somatosensory,and visual association areas of awake adult cats yielded aspontaneous firing rate of 9.4 ⫾ 1.7 Hz (Steriade et al. 2001).
The functional role of this spontaneous firing of neocorticalneurons is unknown. It was previously conjectured that persistent neuronal activity, which is “conspicuously absent in cerebellar and basal ganglia circuits,” is an essential feature of the
neocortex, which enables it “to incorporate the past into thesystem’s present state” (Buzsáki 2006). It was also previously
conjectured to support “flexible cooperation among local anddistant cell assembles,” which is “believed to underlie theefficacy of cortical performance and is an essential ingredientfor cognition” (Buzsáki 2006). On a more technical level, it
was pointed out that not only recurrent circuits of neurons in
the neocortex, but dynamical systems in general support higher
computational performance if they operate not in the ordered(or dissipative) regime, but closer to a critical state—whichamounts to a substantial spontaneous firing activity in a neuralcircuit—so that also small external inputs can cause significant
Address for reprint requests and other correspondence: W. Maass, Institutefor Theoretical Computer Science, Technische Universitaet Graz, A-8010
Graz, Austria (E-mail: maass@igi.tugraz.at).
www.jn.org

system responses (Buzsáki 2006; Legenstein and Maass 2007).
Spontaneous brain activity is commonly assumed to begenerated by a combination of intrinsic electrophysiologicalproperties of single neurons and synaptic interactions in the
network (for a recent review see Destexhe and Contreras2006). It was already pointed out in Abeles (1991) that “itseems to be very difficult to stabilize a network of inhibitoryand excitatory neurons at very low levels of activity.” Throughtheoretical analysis and extensive computer simulations weshow that the short-term dynamics of synapses—more precisely, the empirically found diversity of mixtures of pairedpulse depression and paired-pulse facilitation that depend on
the type (excitatory or inhibitory) of pre- and postsynaptic
neurons— endows networks of leaky integrate-and-fire (LIF)
neuron models with either conductance-based or current-based
synapses, with a remarkable stability of spontaneous network
activity, even for rather low average firing rates, such as 10 Hz.
It is well known that synapses in the neocortex do notrespond in the same way to each spike in a train of spikes.
Rather, the amplitude of the postsynaptic response [excitatory
postsynaptic potential (EPSP) or excitatory postsynaptic current (EPSC)] to a presynaptic spike depends on the history ofpresynaptic firing activity. If the response amplitude for thesecond spike in a pair of spikes is smaller than that for the firstone, one calls this effect paired-pulse depression. If the response amplitude for the second spike is larger, one calls thispaired-pulse facilitation. Actually, however, most neocorticalsynapses exhibit a mixture of depression and facilitation andthe amplitude of the postsynaptic response depends in a complex manner on the temporal pattern of several preceding
presynaptic spikes. This effect is commonly referred to assynaptic dynamics or short-term plasticity. In contrast to longterm synaptic plasticity this short-term plasticity does not
usually depend on the pattern of postsynaptic firing and its
effect is not long-lasting (more precisely: it depends only on
the temporal pattern of presynaptic spikes during the last few
seconds). We refer to Abbott and Nelson (2003), Thomson(2003), and Abbott and Regehr (2004) for recent reviews of
biological data and proposed functional roles of short-time
synaptic plasticity. The dynamics of cortical synapses can be
fitted quite well by the phenomenological model fromMarkram et al. (1998), where the short-term dynamics ofsynapses is characterized by three parameters: U (whichroughly models the release probability of a synaptic vesicle forthe first spike in a train of spikes), D (time constant forrecovery from depression), and F (time constant for recovery
from facilitation). Hereafter we call this set of parameters
UDF. An essentially equivalent model was previously proposed in Varela et al. (1997).

0022-3077/07 $8.00 Copyright © 2007 The American Physiological Society

4079

Downloaded from jn.physiology.org on June 12, 2007

Sussillo D, Toyoizumi T, Maass W. Self-tuning of neural circuits
through short-term synaptic plasticity. J Neurophysiol 97: 4079 – 4095,
2007. First published April 4, 2007; doi:10.1152/jn.01357.2006. Numerous experimental data show that cortical networks of neurons are notsilent in the absence of external inputs, but rather maintain a low
spontaneous firing activity. This aspect of cortical networks is likely
to be important for their computational function, but is hard to
reproduce in models of cortical circuits of neurons because thelow-activity regime is inherently unstable. Here we show—throughtheoretical analysis and extensive computer simulations—that shortterm synaptic plasticity endows models of cortical circuits with aremarkable stability in the low-activity regime. This short-term plasticity works as a homeostatic mechanism that stabilizes the overall
activity level in spite of drastic changes in external inputs and internal
circuit properties, while preserving reliable transient responses tosignals. The contribution of synaptic dynamics to this stability can bepredicted on the basis of general principles from control theory.

4080

D. SUSSILLO, T. TOYOIZUMI, AND W. MAASSẋ ⫽ ⫺␮0x⫹␮x ⫹ v共t兲(1)␮˙ ⫽ f共x兲 ⫺ g共␮兲(2)

The parameter ␮0 is here assumed to be unknown and conditions1 on the functions f and g are given under which thevariable x moves (in the absence of external input) to a steadystate x* with ẋ ⫽ 0 (i.e., ␮ ⫽ ␮0) regardless of the value of thisparameter ␮0. The system expressed by Eqs. 1 and 2 is
1

It is shown in Moreau and Sontag (2003) that it suffices to assume that 1)
f is decreasing, 2) g is strictly increasing, and 3) there exists some x* such thatf (x*) ⫽ g(␮0).

obviously quite far removed from models for recurrent neuralnetworks consisting of excitatory and inhibitory neurons with
nonnegative firing rates. However, it is remarkable that its
self-tuning property relies on the assumption that the strength
␮ of the positive feedback in Eq. 2 varies as a function of x,
unlike the usual “static” synaptic weights in artificial neuralnetworks. Furthermore, the required assumptions on the functions f and g in Eq. 2 imply that ␮ decreases when x grows
beyond its steady-state value x*. Thus the resulting dynamics
of the “synaptic weight” ␮ has some properties in common
with synapses between excitatory neurons in the cortex (whichhave been found to be usually depressing), although apparently
no relation to synaptic dynamics had been intended (or evendiscussed) in Moreau and Sontag (2003). Herein we show (seethe end of RESULTS and the Supplementary material2 that thiscontrol principle can be generalized and adapted to more
complex dynamical systems that consist of two interacting
modules (e.g., a population of excitatory neurons and a population of inhibitory neurons) that have profoundly different
functional roles. It turns out that an application of the control
principle to such a distributed system requires not only an
activity-dependent modulation of “weights” ␮(x), but different
rules for the modulations ␮mn(x) of different weights, that
depend on the type n of the source of the feedback ␮mn(xn)䡠xn
and the type m of the target of the synaptic connection (wherem, n range over excitatory and inhibitory neurons). Thus inapplications to models for cortical circuits, this principle requires differential synaptic dynamics for different types ofsynaptic connections.
More precisely, our theoretical analysis predicts that theUDF parameters that characterize the dynamics of synapses
belong to a specific region within the UDF space (which wecall the “N-volume”) if a strengthening of this type of synaptic
connection tends to increase the firing rate of excitatory neurons in the circuit, or into another region (“P-volume”) if astrengthening of this type of synaptic connection tends todecrease the firing rate of excitatory neurons. A comparison ofthese predicted ranges of UDF values with empirically foundUDF values gives a quite satisfactory (although not perfect)
match. This is remarkable because many aspects that are
relevant for the dynamics of cortical circuits of neurons (different input– output behaviors of excitatory and inhibitory
neurons), further differentiation of neuronal dynamics according to subclasses of inhibitory (Gupta et al. 2000) and excitatory neurons, subnetworks formed by particular subtypes ofneurons (Yoshimura and Callaway 2005; Yoshimura et al.
2005), and modulation of neuronal dynamics through neuromodulators (Gulledge and Jaffe 1998, 2001) are not reflected in
the network models consisting of LIF neurons that we examine
herein.
Altogether this article provides a novel explanation for thefunctional role of short-term synaptic dynamics, in particularalso for the empirically found diversity of synaptic dynamics independency on the type of pre- and postsynaptic neurons. In
addition it contributes new ideas and principles for the investigation of a research topic, obviously of primary importance—
the remarkable stability of cortical circuits of neurons against
numerous perturbations and parameter drifts. This stability has
so far not been reproduced with either artificial neural networks
2

The online version of this article contains supplemental data.

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

In contrast to the static synapse of an artificial neuralnetwork, a dynamic synapse does not just provide a fixed
amplification of presynaptic activity, but rather implements a
nonlinear filter that could potentially serve a number of different purposes, such as cortical gain control (Abbott et al. 1997),
burst detection (Lisman 1997), or temporal integration of
information (Haeusler and Maass 2007; Maass et al. 2004). In
Tsodyks et al. (1998) it was shown that depressing synapses
can create in a mean-field model for a population of excitatoryneurons a stable fixed point for its activity and that dynamic
synapses can create a variety of rhythmic and irregular activity
patterns if one adds a population of inhibitory neurons.
One intriguing aspect of biolocial data on synaptic dynamics is that this short-term dynamics—more precisely,
the parameters UDF— differ from synapse to synapse.
Moreover, the UDF values form clusters within the threedimensional (3D) UDF parameter space and the cluster towhich the UDF value of a particular synapse belongs depends on the type (excitatory or inhibitory) of its presynaptic neuron and in some cases also on the type of thepostsynaptic neuron (Galarreta and Hestrin 1998; Gupta etal. 2000; Markram et al. 1998; Thomson 1997). It has beendifficult to identify a possible functional role for the rulesthat assign a particular region in the 3D UDF parameter
space to different types of synaptic connections. Herein weshow that the assignment of different ranges of UDF values
to different types of synapses makes it possible to implement a self-tuning principle for the firing rate of a corticalcircuit. This self-tuning principle enables neural circuits to
respond to external perturbations with a characteristic transient response in the firing rate of excitatory neurons andthen returns to its previous firing rate within a few hundredsof milliseconds back into a given target range. The firing
rate of inhibitory neurons is automatically adjusted by thesynaptic dynamics so that it compensates the external perturbation. A similar self-tuning property can be demonstrated for changes within the circuit, such as those caused
by long-term synaptic plasticity or changes in the concentration of neurotransmitters.
We show that this self-tuning principle in models of cortical
networks of neurons is related to an abstract self-tuning principle recently proposed for much simpler dynamical systems incontrol theory (Moreau and Sontag 2003a,b). Their abstract
self-tuning principle enables some simple two-dimensional(2D) dynamical system to return in a reliable manner to aparticular operating regime (e.g., the vicinity of a bifurcation),in spite of external perturbations. They considered the following dynamical system with external input v(t) that has just twovariables: x (firing rate) and ␮ (synaptic strength)

SELF-TUNING OF NEURAL CIRCUITSMETHODS

Neuron models
We used two leaky integrate-and-fire neuron models (LIF): one
with current-based synapses and the other with conductance-based
synapses. The current-based LIF neuron is defined by
␶ mem

dV共t兲
⫽ ⫺ 关V共t兲 ⫺ Vrest兴 ⫹ Rm
dt

冋冘
KE

冘
KI

jE, j共t兲 ⫹

j⫽1

j⫽1

册

jI, j共t兲 ⫹ Iext共t兲

(3)where V is the membrane voltage; ␶mem is the membrane time
constant; Vrest is the resting membrane voltage; Rm is the membrane
resistance; JE,j (t) and JI,j (t) are the KE and KI synaptic currents fromthe excitatory and inhibitory synapses, respectively; and Iext(t) is anexternal input composed of mean input current Iinject and a white noise
current Inoise(t), with zero mean and SD ␴noise. The LIF neuron with
conductance-based synapses is defined by
Cm

dV共t兲
⫽ ⫺ gleak关V共t兲 ⫺ Vrest兴 ⫺
dt

冘
KE

gE, j共t兲关V共t兲 ⫺ EE兴

j⫽1

冘
KI

⫺

gI,j共t兲关V共t兲 ⫺ EI兴 ⫹ Iext共t兲

(4)

j⫽1

where Cm is the membrane capacitance; gleak is the leak conductanceof the neuron; and gE,j(t) and gI,j(t) are the KE and KI synaptic
conductances from the excitatory and inhibitory synapses, respectively. The constant EE is the reversal potential for the excitatory
synapses and the constant EI is the reversal potential for the inhibitory
synapses. The dynamics of each synapse over time are defined by
j 共t兲
dj 共t兲
⫽ ⫺⫹ W␦共t ⫺ tsp兲
dt
␶syn

(5)

where j(t) is the instantaneous synaptic current or conductance,
depending on the model, and W is the synaptic weight. The
currents or conductances decrease exponentially with time constant␶syn and increase instantaneously by adding W to the running valueof j(t) whenever a spike occurs in the presynaptic neuron at time␶sp. Herein we compare the use of static synapses to dynamic
synapses and so the synaptic weight W was constant in somesimulations and in other simulations allowed to vary according toknown biological rules of synaptic plasticity [i.e. W(t)]. In the caseof dynamic conductance-based synapses we varied the maximal
conductance and in the case of dynamic current-based synapses wevaried the maximal current. When static synapses were used (i.e.,
when W is constant) we use the notation J to denote the synaptic
weight and always state whether the simulation used conductanceor current-based synapses. When dynamic synapses were used
[i.e., when W(t) is implied] we use the notation ␮ to denote thesynaptic weight.
We chose values for these LIF neurons that qualitatively reflect thehigh conductance, in vivo UP state measured in in vivo studies ofcortical neurons (Destexhe et al. 2003). The membrane resistance was
Rm ⫽ 10 M⍀, the membrane time constant was ␶mem ⫽ 10 ms, and,in the case of neurons with conductance-based synapses, the capacitance was set to Cm ⫽ 10 ms/10 M⍀ ⫽ 1 nF. The neurons’ resting
potential was Vrest ⫽ ⫺60 mV and the firing threshold was Vthresh ⫽⫺50 mV. When the neuron spiked, the membrane voltage was reset to
Vreset ⫽ ⫺60 mV, where it remained for a refractory period of ␶refr ⫽3 ms. The excitatory synaptic time constant ␶syn was ␶e ⫽ 4 ms; theinhibitory value was ␶i ⫽ 8 ms. In the case of neurons with conductance-based synapses gleak ⫽ 1/Rm ⫽ 100 nS and the reversal
potentials for the excitatory and inhibitory neurons were EE ⫽ 0 mV
and EI ⫽ ⫺80 mV. We modeled background synaptic inputs (fromneurons outside of the simulated circuit) through Iext ⫽ Iinject ⫹ Inoisewith a constant Iinject ⫽ 2.455 ⫻ 10⫺9 A and values of Inoise drawn
from a 0-mean Gaussian with ␴noise ⫽ 6 ⫻ 10⫺9 A SD. Iext caused anaverage membrane potential of ⫺55.4 mV and SD of 4.3 mV in theabsence of recurrent synaptic feedback (and when the threshold Vthresh
was temporarily raised to 0 mV to prevent spikes). This resulted in a
firing rate of approximately 20 Hz for each neuron in the absence ofrecurrent synaptic input Irec and was chosen because it is the medianvalue of the firing rates of irregular and tonic discharges reported in invivo conditions in waking animals (Destexhe et al. 2003; Steriade etal. 2001) and is thus suitable for studying stability involving both upand downregulation of firing rates.

Spiking network modelIn our simulations we created sparse, random-spiking networks of
5,000 neurons, with 4,000 excitatory neurons and 1,000 inhibitory
neurons. Neurons were either excitatory or inhibitory. The networkmodels contained either purely LIF neurons with conductance-based
synapses or purely LIF neurons with current-based synapses; we referto the former as “conductance-based networks” and the latter as“current-based networks.” We use the term sparse to denote that the
number of excitatory connections Ke and inhibitory connections Ki toa given neuron are much less than both the total number of excitatoryneurons Ne or inhibitory neurons Ni in the network. We use the term
random to denote that the connections between neurons are chosenrandomly and follow no other design principle. Following Vogels and
Abbott (2005), the neurons were connected randomly to each otherwith a 2% connection probability. The network connectivity between
neuron types is shown in Fig. 1.

Dynamic synapse model
We used two equivalent forms of dynamic synapse models for
networks with spiking and nonspiking neurons. The first set of
expressions constituting the Eq. 6 system is taken from Markram et al.

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

or with biologically more realistic models for cortical circuits
and thus an understanding of its prerequisites may provide
important insight into design principles of cortical networks.
The general program to investigate fundamental aspects of thedesign of neural circuits by focusing on their stability properties has already produced remarkable insight into another typeof neural circuit: the pattern-generating pyloric circuit in thecrustacean stomatogastric ganglion (for a recent review see
Marder and Goaillard 2006). It turns out that this circuit isendowed with intricate homeostatic mechanisms that move ahigh-dimensional vector of internal parameters (conductances,etc.) to a different setting if one method for producing thedesired firing pattern becomes unavailable, such as because ofchanges in the concentration of neuromodulators. Thus it isessential for the functioning of this pattern-generating circuitthat it can produce “multiple solutions to the production ofsimilar behavior” (Marder and Goaillard 2006). We proposethat analogous design principles are implemented in corticalcircuits of neurons, which have to maintain a stable computational function in spite of drastic changes in the intensity of
external inputs and in the concentrations of neuromodulators(that affect the excitability of neurons within the network;
Gulledge and Jaffe 1998, 2001). More specifically, we proposethat the short-term dynamics of synapses is essentially involved in design principles that produce a stable function ofcortical networks.

4081

4082

D. SUSSILLO, T. TOYOIZUMI, AND W. MAASS

on the current presynaptic firing rate x(t), but through the auxiliary
variables u and R also on the history of this presynaptic firing rate.
The steady-state equations for this continuous dynamic synapse
model yield in the case of a constant presynaptic firing rate x(t) ⫽ x*
␮ * ⫽ AR*U1*
u* ⫽

FUx*
1 ⫹ FUx*

R* ⫽1
1 ⫹ DU1x*

U 1* ⫽ u*共1 ⫺ U兲 ⫹ U(8)where the * notation denotes steady-state values of u, U , R, and ␮. In
particular for m, n 僆 {e, i} the value ␮*mn (xn) defined according to Eq.
8 denotes the steady-state synaptic weight from any neuron in population n to any neuron in population m as a function of the presynaptic
firing rate x* ⫽ xn.
1(1998)3 and describes for the case of spiking neurons the magnitudeof the synaptic weight ␮k (a current or conductance, depending on the
network model) for the kth spike in a spike train with interspikeintervals ⌬1, ⌬2, . . . , ⌬k⫺1
␮ k ⫽ AR kuk
u k ⫽ U ⫹ uk⫺1共1 ⫺ U兲 exp共 ⫺ ⌬k⫺1/F兲
R k ⫽ 1 ⫹ 共Rk⫺1 ⫺ uk⫺1Rk⫺1 ⫺ 1兲 exp(⫺⌬k⫺1/D)
u1 ⫽ U
R1 ⫽ 1(6)where u is the running variable for synaptic utilization and R is therunning variable for synaptic availability, with u, R 僆 [0, 1]. The
constants U, D, and F (for curtness further denoted as UDF or AUDF
when the strength scaling factor A is included) represent the releaseprobability for the first spike, the depression time constant, andfacilitation time constant, respectively. The synaptic weight for the kth
spike is then defined by ␮k (i.e., W :⫽ ␮k in Eq. 5) and has an arbitrarystrength scale set by A.
An equivalent, continuous mean-field model for dynamic synapses
was developed in Tsodyks et al. (1998), which can be used innonspiking mean-field models with continuous time. The mean-field
equations for dynamic synapses are given by␮ ⫽ ARU 1
u̇ ⫽ ⫺

Linking UDF parameters to network stability
We show that the fundamental property of dynamic synapses forcreating a network with a stable firing rate is the relationship between
the slopes of the steady-state synaptic weight curves ␮*mn (xn) for
different synapse types. We searched for UDF parameters that generate a negative derivative for the connections E 3 E and I 3 I
because for our parameters an increase in the strength of either leads
to a higher firing rate xe of the excitatory neurons and a positive
derivative for the steady-state curves for the connections E 3 I andI 3 E because an increase in their strengths leads to a lower firing ratexe (see Supplemental materials for further explanation). More specifically, we searched for UDF settings that generated either a strictly
positive or strictly negative derivative between the frequencies 10 and
100 Hz of the presynaptic firing rate and that additionally were in aphysiological range (see Table 1). We uniformly sampled the UDF
space in the range U 僆 (0, 1], D 僆 (0, 1] s, and F 僆 (0, 1] s (atTABLE 1.

Parameter

UDF parameters
Synapsis
Type

U

D, s

F, s

E3E
E3I
I3E
I3I
E3E
E3I
I3E
I3I
E3E
E3I
I3E
I3I
E3E
E3I
I3E
I3I

0.59390.40280.00070.50890.63190.15170.07460.30290.57620.10100.08650.5521
0.59
0.049
0.016
0.25

0.53330.00160.11530.1744
0.94680.00630.00010.44290.61870.01050.00040.6139
0.813
0.399
0.045
0.706

0.18280.08480.17950.49730.99490.27010.90430.99630.79890.10030.57790.4220
0.001
1.79
0.376
0.021

R1

R2

R3

u
⫹ U共1 ⫺ u兲x共t兲
F

Experimental

1⫺R⫺ U1Rx共t兲
Ṙ ⫽D
U 1 ⫽ u共1 ⫺ U兲 ⫹ U(7)where x(t) represents the instantaneous firing rate of the presynaptic
neuron; U1 is simply an algebraic term introduced for clarity; and ␮,u, R, U, D, and F have the same meaning as in the Eq. 6 system. Notethe magnitude ␮(t) of the synaptic weight at time t depends not only3

The variable uk⫺1 in the third expression in the Eq. 6 system waserroneously replaced by uk in Markram et al. (1998); see the discussion in
Maass and Markram (2002).

UDF parameters—R1, R2, and R3—were chosen randomly based on theoretical criteria (see METHODS). The UDF parameters for synapse types E 3 I
and I 3 E were randomly chosen from the P-volume (Fig. 4) and the UDF
parameters for synapses E 3 E and I 3 I were chosen randomly from the
N-volume. Experimental dynamic synapse parameters—U, D, and F values—
were from the experimental literature. The E 3 E and E 3 I synapse types had
UDF values taken from the means given in Markram et al. (1998). The I 3 E
and I 3 I had UDF values taken from the means given in Gupta et al. (2000).
The F1 (facilitating) parameters were used for the I 3 E synapses and the F2(depressing) parameters were used for I 3 I synapses.

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

FIG. 1. Basic network architecture used in mean-field simulations andspiking network simulations. Two populations of neurons connect to eachother, one excitatory Exc with firing rate xe and one inhibitory Inh with firing
rate xi. Blue denotes excitatory populations or recurrent excitatory connections
and red denotes inhibitory populations or inhibitory recurrent connections. If
the synaptic connections have a static synaptic weight they are denoted Jmn
from presynaptic pool n to postsynaptic pool m; if the synaptic connections
have dynamic synaptic weights then they are denoted as ␮mn. External
excitatory input, labeled Iinject ⫹ Inoise, is denoted by green. External input inthis paper is a constant noisy current with a stationary mean.

SELF-TUNING OF NEURAL CIRCUITS

1,000 ms. The networks with dynamic synapses ran for 2,000 ms andthe steady-state firing rates were measured for the last 1,000 ms,
which allowed extra time for the dynamic synapses to equilibrate. The
measured firing rates for the theoretically derived UDF values werethen averaged across all three dynamic synapse parameters sets R1,
R2, and R3. Network simulations were performed using the CSIM
neural simulator (Natschläger et al. 2003) and MatLab (The MathWorks, Natick, MA). The simulation time step was 0.1 ms. The
mean-field simulations were also performed using MatLab.

Correlation function
Network synchrony plays a large role in the behavior of networksof sparsely connected, spiking neurons (Brunel 2000). We wereinterested in separating asynchronous irregular activity from othertypes of network activity. We used the voltage traces of randomly
sampled neurons to analyze the network correlation and then applied
the standard cross-covariance function. The covariance coefficient forthe sequences y and z, with means y៮ and z៮ , defined at the zeroth lag,is given by

冘
冑冘
N⫺1

C yz共0兲 ⫽共yn ⫺ y៮ 兲共zn ⫺ z៮兲

n⫽0

N⫺1

n⫽0

冘

N⫺1

共 yn ⫺ y៮ 兲2

共zn ⫺ z៮ 兲2n⫺0

To measure the network correlation, we created two random
groups of neurons from the excitatory neuron population (excepting differences in firing rate, the inhibitory neurons behave in samemanner as the excitatory neurons). The two groups each consisted
of 5% or 250 neurons. The cross-covariance coefficient was computed for every pair of neurons, one from each group (excluding
identical neurons, which occurred by chance) and the average crosscovariance coefficient was taken as a measure of the network
correlation.

Self-tuning simulation setup

Mean-field models for sparsely connected neural circuits

The self-tuning simulations we performed involved measuring andcomparing the equilibrated (steady-state) average firing rates of excitatory neurons with static synapses and then with dynamic synapses
for three network models: conductance-based networks, current-based
networks, and a rate-based mean-field network (see following text forthe definition). There were three types of perturbation simulations run,with two parameters perturbed in each case:1) mean input current strength and SD of current injection2) excitatory synaptic strength and inhibitory synaptic strength3) percentage excitatory neuron inactivation (a temporary cessation of function potentially caused by a variety of circumstances suchas neuromodulators) and percentage inhibitory neuron inactivation.
This resulted in 841 networks for the input perturbation experiments
and 900 networks for the synaptic weight perturbation and neuronal
inactivation experiments. Each network was generated with random
connections and first run in the control case with static synapses. The
network was then regenerated with the exact same parameters, but adifferent random number seed, and run with dynamic synapses. To
test the robustness of the ideas presented, each network was tested
with four different dynamical synapse parameters sets, three based ontheoretical considerations called R1, R2, and R3 and the parameters
based on experimental data, thus requiring 841 ⫹ (841 ⫻ 4 ⫻ 3) ⫽10,933 separate network simulations for the input perturbation experiments and 900 ⫹ (900 ⫻ 4 ⫻ 3) ⫽ 11,700 separate network
simulations for experiments with synapse strength and neuronal inactivation perturbations. To measure the firing rates we let the control
networks with static synapses run for 1,500 ms and the average firing
rates of excitatory and inhibitory neurons were measured for the last

For our analysis of network stability we developed neuronal population models—that is, mean-field firing rate models—for sparse,
random networks (van Vreeswijk and Sompolinsky 1996, 1998). Theymodel populations of current-based neurons with static and dynamic
synapses and are based on a balance of excitation and inhibition thatcreates temporally irregular spiking. We assume that there are Ne
excitatory neurons and Ni inhibitory neurons. Kmn is the number ofsynaptic connections from neurons in pool n to neurons in pool m (weassume that Kmn ⬍⬍ Ne, Ni). We are interested in networks withrandomly chosen sparse synaptic connections because they provide
reasonable first approximation models for cortical circuits.
Because the firing rate of a neuron is dependent on both the mean input
current and the variance of the input current, we refer to its input– output
function as the FMS (firing rate, mean input, SD) surface, which wedenote by F. A sparse, random network has a parameter range that allows
for asynchronous irregular activity (Brunel 2000; for a review, see Vogels
et al. 2005) and the synaptic currents to a given cell are uncorrelated intime for this regime, aside from a potentially fluctuating population firingrate. In the steady-state condition, which is when all the dynamics of thesystem have been allowed to equilibrate, this allows one to use the time
average of the synaptic input for estimating the synaptic input at eachpoint in time. Thus in the asynchronous condition, and assuming Poisson
spiking statistics, one can write down the average recurrent synaptic
current Irec and recurrent synaptic input variance ␴rec2 for each neuron in
their respective pools m (Shriki et al. 2003) as

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

I rec
m ⫽

冘
n

Kmn␶nxn Jmn

Downloaded from jn.physiology.org on June 12, 2007

increments of 0.014 in each dimension) and found two topologicallysimple volumes of parameter space that satisfied these requirements.
These were designated “P-volume” for the positive derivative and the“N-volume” for the negative derivative. We randomly chose threeassignments R1, R2, and R3 to the UDF values of the four synapse
types to create three separate dynamic synapse UDF parameters sets.
Thus each parameter set Ri (i ⫽ 1, 2, 3) consisted of two points fromthe N-volume, which were randomly assigned as the mean for the
UDF-parameters of the synaptic connections E 3 E and I 3 I, andtwo points from the P-volume, which were also randomly assigned as
the mean to the synaptic connections E 3 I and I 3 E. Eachindividual Ri parameter set was used separately from the others in the
network simulations. Additionally, we used experimental data on
UDF parameters from Gupta et al. (2000) and Markram et al. (1998)
for all four synapse types (see Table 1). All self-tuning simulations
were carried out for these parameters based on experimental data, as
well as all three parameter assignments R1, R2, and R3 based on atheoretical analysis of desirable network performance. Performancedifferences between the R1–R3 parameter settings were small and weaveraged the responses across the parameters sets. These UDF parameters were fixed for all self-tuning experiments and only thescaling factors Amn of the dynamic synapses were modified on a
per-simulation or per-network basis.
Altogether synaptic parameters were chosen as follows. Allsynapses were given a uniform delay of 0.1 ms. The means of thescaling factor parameters (the A terms in Eqs. 6 and 7), mean UDF
parameters, and mean values for the initial conditions for dynamic
synapses u0 and R0 were set as described earlier. After selection ofa mean value, each actual value was drawn for each synapse in the
network from a Gaussian distribution with SD of 10% of the meanvalue (e.g., all UDF parameters for every synapse were chosenfrom a Gaussian distribution with a SD equal to 10% of the mean
UDF parameters, which depended on the synapse type). If theactual value would have been set below zero, then the value wasredrawn from a uniform distribution between zero and two times
the mean.

4083

4084

D. SUSSILLO, T. TOYOIZUMI, AND W. MAASS

␴ rec2⫽m

1
2

冘K mn␶nxn J 2mn(9)nwhere m, n 僆 {e, i} represent indices into the post- and presynaptic
pools, respectively. So the total average current and total input current
variance into a neuron in a pool m isrecinjectI tot⫹ I noisem ⫽ Im ⫹ I
rec2
␴ tot2⫹ ␴ inject2 ⫹ ␴ noise2m ⫽ ␴m

(10)

␶ mem

dx e
tot
⫽ ⫺ xe ⫹ F៮ 共I tot
e , ␴e 兲
dt

␶ mem

dx i
tot
⫽ ⫺ xi ⫹ F៮ 共I tot
i , ␴i 兲
dt(11)

where xe is the average firing rate of the excitatory neurons and xi isthe average firing rate of the inhibitory neurons. For a current-based
LIF network with ␶e, ␶i ⬍ ␶mem, the correct time constant of the modelis the neuronal membrane time constant ␶mem.
This model can also be used to analyze system dynamics, eventhough the expressions in the Eq. 11 system were developed from an

I rec
m ⫽

冘
冘

Kmn␶nxn␮mn共xn, t兲n␴ rec2⫽m

1
2

K mn␶nxn␮mn共xn, t兲2

(12)

n

where have replaced the static synaptic weight value Jmn with the
dynamic synaptic weight value ␮mn(xn, t).
RESULTS

Conditions on dynamic synapses to create network stability
The goal of introducing dynamic synapses into the sparse,
random model is to set up a stable fixed point at a desired target
network firing rate xe ⫽ xi ⫽ x*, so that xe is robust to various
types of perturbations (i.e., the circuit is self-tuning). We view
self-tuning as the ability of a system to ultimately return asystem value to given range, despite a perturbation in systemparameters or input. We define the output of the network to be
the average excitatory firing rate xe, which is the firing rate thatis actively controlled, perhaps by increasing or decreasing xi.
Although sparse, random networks also have steady-state firing
rates as a function of static synaptic weight values Jmn, we willshow that introducing dynamic synapses with relatively general constraints endows the network with much greater robustness under a variety of perturbations. This improved robustness

FIG. 2. FMS (firing rate, mean current injection, SD of
current noise injection) surface surface for leaky integrateand-fire (LIF) neuron. A: approximated FMS surface F៮ of LIF
neuron with ␶m ⫽ 10 ms, ␶rp ⫽ 3 ms, Vrest ⫽ ⫺60 mV, andVreset ⫽ ⫺60 mV. Approximated F៮ gives a very good fit withmean absolute error 0.2 Hz. Dashed black line shows thefiring rate as a function of the external current noise levelassuming a constant current injection of 2.455 ⫻ 10⫺9 A.
Solid black line shows the firing rate as a function of mean
input, with injected noise value of 6 ⫻ 10⫺9 A. Intersection ofthese 2 lines denotes the external input (mean current andnoise level) that every neuron in all simulations of spiking
neurons received and resulted in a steady-state firing rate of
20 Hz in the absence of synaptic feedback.

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

where Iinject and ␴noise2 represent the mean and variance of external
tot2
input to the neuron, and Itotm and ␴m are the mean and variance of thetotal current injection to each cell (see METHODS).
To use these formulas to create a working neuron model one caneasily sample the FMS surface of a single current-based neuron by
injecting current with an arbitrary mean and Gaussian noise component and then use the Eq. 9 system to translate the parameters (Jee, Jie,
xe, Jei, Jii, xi) into (Itot, ␴tot). Using this technique we sampled the
FMS surface F of a current-based LIF neuron with the same neuronalparameters as the neurons used in our current-based spiking networks.
Using the samples generated from this neuron we used smoothing
splines to create a nonparametric approximation F៮ to F that could be
evaluated at any given point to yield a good approximation of thepostsynaptic firing rate (Fig. 2). Because the excitatory and inhibitory
neurons in our current-based networks have the same neuronal parameters, we used the same neuron model for both types. Finally, themean-field model for a network of excitatory and inhibitory neurons
in an asynchronous network, which also represents a sparse, random
network under the same assumptions, is defined byanalysis of the steady state, because the neuronal dynamics of acurrent-based neuron acts merely as a low-pass filter. A major assumption of this approach is that all combinations of (Jee, Jie, xe, Jei,
Jii, xi) lead to asynchronous, irregular network firing, which is notalways the case for sparse, random-spiking networks. Thus for themean-field model to correctly predict the behavior of the spiking
network, the range of parameters must be limited to ensure thisassumption and, consequently, the network correlation and coefficientof variation (CV) of interspike intervals were measured.
Finally, we introduce dynamic synapses (Eq. 7 system) into thesparse, random mean-field model of the Eq. 11 system. The model isdefined by the expressions in Eq. 11 with slightly different definitions
rec2
for I recm and ␴ m

SELF-TUNING OF NEURAL CIRCUITS

0.25

0.2

0.2
µ* (x )eeµ* (x )
ei

*
ie

e

0.15

µ*ii(x )

*
ii

0.15

µ , µ (nA)*ee

Input to Inhibitory Population

0.25

*
ei

µ ,µ (nA)

B

Input to Excitatory Population

A

0.1

4085

0.1

0.05

0.05

0

0

µ* (x )
ie

0

5

10
Frequency (Hz)

15

20

0

5

10
Frequency (Hz)1520

is achieved if the UDF parameters (the parameters of the
dynamic synapse model; see METHODS) of the dynamic synapses
are chosen to satisfy the following three conditions.4
1) Maintain the synaptic weight values if the presynaptic
firing rate is equal to the desired firing rate.
2) Increase the strength of the E 3 E and I 3 I synapses
and decrease the strength of the E 3 I and I 3 E synapses if
the presynaptic firing rate is below the target firing rate.
3) Decrease the strength of the E 3 E and I 3 I synapses
and increase the strength of the E 3 I and I 3 E synapses if
the presynaptic firing rate is above the target firing rate.
The first constraint involves creating the fixed point at atarget firing rate by setting the scaling factor A of the dynamic
synapse model according to the formula
A mn ⫽Jmn
R*mn共x*兲U*mn1 共x*兲(13)where Jmn represents specific synaptic weight values of acorresponding network with static synapses and R*mn (x*) andU1*
mn (x*) are the steady-state values for the internal dynamic
synapse variables as a function of presynaptic firing rate, here
set to the desired steady-state firing rate (see Eqs. 6 and 7 in
METHODS). Equations 6 and 13 ensure that the steady-state
dynamic synapse weight is equal to the static synaptic weight[␮*mn (x*) ⫽ Jmn] for a presynaptic firing rate x*. We refer to␮*mn (x*) as the steady-state value for ␮mn (see METHODS, Eq. 8)
and designate the resulting function of x* as the steady-state
synaptic weight curve (see Fig. 3). In the asynchronous range,
when static synaptic weight values (Jee, Jei, Jie, Jii) are set fora sparse, random network, the network has a specific steadystate firing rate associated with those synaptic strengths, x*(van Vreeswijk and Sompolinsky 1996, 1998). Because the
network already has a fixed point with firing rate x* for thesestatic weights, we simply scale the dynamic synapses to takethese static weights at x*. We do this by setting the value ofAmn according to Eq.13, which gives ␮*mn (x*) ⫽ AmnR*mn
4
See the Supplemental materials for a theoretical justification of thesechoices.

(x*)U1*(x*) ⫽ Jmn[R*mn (x*)U1*(x*)/R*mn (x*)U1*(x*)] ⫽ Jmn,
thus returning the firing rate x*.
Satisfying the second and third constraints involve setting
the UDF parameters of the dynamic synapses such that if there
is a perturbation in the network that affects the steady-statefiring rate then the dynamic synapses attempt to compensate
and bring xe back to x*. The UDF parameters are responsible
for the shape of the steady-state synaptic weight curve. Thus tocreate an attractor at x*, which produces firing rate stability, wechoose UDF parameters such that the synapse types E 3 E andI 3 I have negative derivatives in their steady-state synaptic
weight curves and the synapse types E 3 I and I 3 E havepositive derivatives in their steady-state synaptic weight curves(see Fig. 3). To test the validity of this analysis, we ran aparameter search over the range of physiological UDF values
to determine whether there were simple volumes of UDF spacethat fulfill these two requirements through the range of firing
rates 10 –100 Hz. (See METHODS for a description of how andwhy these values were chosen.) Based on these theoreticalconsiderations we chose three sets of UDF parameters (R1, R2,
and R3) for each of the four synapse types randomly from theappropriate subvolume of UDF space: for synaptic connections
E 3 E and I 3 I six points were chosen from the blue
N-volume, for E 3 I and I 3 E six points were chosen fromthe red P-volume (see Fig. 4, A and B). All resulting UDF
parameters are listed in Table 1. The three conditions explained
earlier create a stable fixed point for the firing rates of theresulting dynamical system. At the target firing rate xe ⫽ xi ⫽x*, the dynamic synapses equilibrate to (Jee, Jei, Jie, Jii), which
maintains the firing rate. At a firing rate less than x* thesynapses E 3 E and I 3 I increase their absolute synaptic
weight values to values greater than Jee and Jii and the synapses
E 3 I and I 3 E decrease their absolute strength to less thanJei and Jie, leading to an increased excitatory firing rate. At afiring rate greater than x* the synapses E 3 E and I 3 I
decrease their absolute synaptic weight values to less than Jee
and Jii and the synapses E 3 I and I 3 E increase theirabsolute synaptic weight values to greater than Jei and Jie,
leading to a decreased excitatory firing rate.

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

FIG. 3. Steady-state synaptic strengths.␮*
mn are used to control the excitatory firing rate in a network with dynamic synapses. A: synaptic strengths to theexcitatory population as a function of presynaptic firing rate. Dashed, straight lines show static synapse values such that xe ⫽ xi ⫽ 10 Hz (blue, Jee; red, Jei).
Thicker, curved lines show the dynamical synapse steady-state strengths tuned such that exactly at 10 Hz (black “✕”), the static and dynamic synapse values matchup [blue, ␮*ee(xe); red, ␮*ei(xi)] and therefore a network with dynamic synapses also supports the 10-Hz excitatory and inhibitory firing rate of the network withstatic synapses. B: same as A but for the inhibitory population. Note the inversion of the slopes of the steady-state synaptic strengths. Dashed, straight lines: blue,
Jie; red, Jii; thicker, curved lines: blue, ␮*ie(xe); red, ␮*ii(xi). Arrows in A show how the dynamic synapses will push the excitatory firing rate toward 10 Hz in
the event of a system perturbation. If a perturbation causes xe to increase, then ␮*ee(xe) decreases (A) and ␮*ie(xe) increases (B), thereby causing a decrease in xe
and an increase in xi; the latter induces an increase in ␮*ei(xi), causing an additional decrease in xe by xi. If a perturbation causes xe to decrease then the same
curves explain an increase in xe.

4086

D. SUSSILLO, T. TOYOIZUMI, AND W. MAASS

In addition, we took experimentally found UDF values
reported in Markram et al. (1998) and Gupta et al. (2000) andcompared their contribution to network stability with that of
the theoretically deduced parameters that were selected aspreviously described. We sampled random points from Gaussian distributions with means and SDs as reported in thesearticles and plotted the resulting UDF values in Fig. 4, C–F.
The UDF values for E 3 E synapses from Markram et al.
(1998) (blue circles in Fig. 4, C and D) lie clearly in the
N-volume. The UDF values for E 3 I synapses (green circles)from that study lie much closer to the P-volume, but are notcontained in it. The UDF values of GABAergic synapses oftype F1 (facilitating) , F2 (depressing), and F3 (recovering)
reported in Gupta et al. (2000) are contained in the N-volume(see green, cyan, and black pluses in Fig. 4, E and F), but thosefor types F1 and F3 are close to the P-volume. This ambiguous
position of UDF parameters for GABAergic synapses could be
related to the fact that, from a theoretical perspective for
improved homeostatic plasticity, the I 3 I and I 3 E connections may have either positive or negative derivatives depending on other network parameters (see Supplemental materials).
In Wang et al. (2006) several additional clusters of UDF values
for E 3 E synapses in prefrontal cortex are reported. These

values (plotted in Fig. 4, G and H) all fall clearly into the
N-volume (into which they should fall according to the proposed theoretical analysis).
Network parameter ranges
To find ranges of synapse strengths for static synapses
that were compatible with low firing rates and low network
correlation in current-based or conductance-based networks,
we varied the excitatory and inhibitory synaptic weights ofconductance-based and current-based networks (see METHODS and Fig. 5). Specifically, we varied two synaptic weightvalues, one excitatory synaptic weight Je and one inhibitory
synaptic weight Ji, and set the all synaptic weight values
such that Jee ⫽ Jie ⫽ Je and Jei ⫽ Jii ⫽ Ji. So that results
between conductance- and current-based networks could becompared we scaled the conductance-based synapses so thatthe conductance- and current-based synapses had equal
strength when the conductance-based synapses were multiplied by their respective driving forces at the neurons’
resting potential. We found a broad parameter range ofsynaptic weights that could potentially support dynamic
synapses as a network feedback mechanism in the asynchronous, irregular regime. However, because we measured the

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

FIG. 4. Two volumes in the UDF [designation for collective set of parameters: U (which roughly models the release probability of a synaptic vesicle
for the first spike in a train of spikes), D (time constant for recovery from depression), and F (time constant for recovery from facilitation)] parameter
space where the steady-state synaptic weights.␮*mn (xmn) have either positive or negative derivatives with respect to presynaptic firing rate, which are essential
for the theoretical prediction of resulting self-tuning properties of the resulting network. A and B: randomly generated UDF parameters (E 3 E and I 3 I, denotedby blue “*”; and I 3 E and E 3 I denoted by red “*”), which are used to test the theoretically predicted stability properties of UDF values from these 2 volumes.
A: “N-Volume” of UDF space (blue) denotes that parameter range in which ␮*(xpre) has a negative derivative between the firing rates of 10 and 100 Hz (alsoin C, E, and G). B: “P-Volume” of UDF space (red) denotes that parameter range in which ␮*(xpre) has a positive derivative between the firing rates of 10 and
100 Hz (also in D, F, and H). C and D: sampled UDF values for pyramid to pyramid (blue circles) and pyramid to interneuron (green circles) connections basedon means and SDs from Markram et al. (1998). E and F: sampled UDF values for F1 facilitating synapses (green “⫹”), F2 depressing synapses (cyan “⫹”), andF3 recovering synapses (black “⫹”) based on means and SDs from Gupta et al.(2000). These are all UDF values of GABAergic synapses whose targets areeither excitatory or inhibitory neurons; thus no clear prediction for the ranges of their values can be made through the proposed theoretical analysis. UDF valuesof F1 and F2 synapses lie close to the border region between the P- and N-volume. G and H: sampled UDF values for medial prefrontal cortex
excitatory-to-excitatory synapses (PFCe1, PFCe1a, and PFCe1b samples as yellow“✕”; PFCe2, PFCe2a, and PFCe2b samples as black “✕”; and PFCe3, PFCe3a,and PFCe3b samples as magenta “✕”) based on means and SDs from Wang et al. (2006). All these UDF values lie in the N-volume, as predicted by the theoreticalanalysis.

SELF-TUNING OF NEURAL CIRCUITS

200

1

100

0.5
2

4

6

8

1.5
1

0

2

100
0

Inh. syn. strength (nA)

6

8

1.5

0.03

1

0.02

0.5

0.01

0

2

CV

1

0.08
0.06
0.04
0.02

0.5

0.05 0.1 0.15

4

6

8

0

Inh. syn. strength (nS)

F 0.1

Network Correlation

0.03

0.08
0.06
0.04
0.02

0

Inh. syn. strength (nA)

0.02
0.01
0.05 0.1 0.15
Inh. syn. strength (nA)Firing Rate (Hz)

200

0.08
0.06
0.04
0.02

100

0.05 0.1 0.15

00FIG. 5. Parameter search over synaptic
weights for conductance-based spiking networks, current-based spiking networks, andmean-field firing-rate networks with static synapses. y-axis is Je ⫽ Jee ⫽ Jie and the x-axis is
Ji ⫽ Jei ⫽ Jii. A–C: network statistics for
conductance-based spiking networks withstatic synapses. A: network steady-state firingrate. B: coefficient of variation (CV) for conductance-based spiking networks over therange of synaptic weights. C: network correlation of excitatory population as measured as
described in METHODS. There was no significant network correlation. D–F: network statistics for current-based spiking networks withstatic synapses. D: network steady-state firingrate. Five colored circles represent the currentbased spiking network settings used in thetarget firing-rate simulation. E: CV for currentbased spiking networks over the range of synaptic weights. F: network correlation of excitatory population. There was no significant
correlation over the range of synaptic weights
tested. G: network steady-state firing rate formean-field firing rate model.

Inh. syn. strength (nA)firing rate, network correlation, and coefficient of variationof the interspike interval (ISI CV) for only a 2D manifold of
the four-dimensional (4D) space of synaptic weights (Jee,Jei, Jie, Jii), these parameter maps serve only as a guide withwhich to start our simulations because the parameter spacebecomes 4D when the dynamic synapses are allowed to vary
the synaptic strengths according to individual dynamic synapse formulas with different AUDF parameters. We always
started our networks with synaptic weight initial conditions
from these 2D maps, so that the networks started in a setting
where the network behavior is asynchronous irregular andthe feedback current magnitudes are on par with or smaller
than the level of the background current injection. Based onthese maps, we used the synaptic weight parameter range
shown in Fig. 5, A, D, and G as the starting point for the restof our conductance- and current-based spiking network
simulations, as well as our sparse, random mean-field simulations of the Eq. 11 system, all of which used dynamic
synapses. Specifically, Jee, Jie 僆 [0, 0.1] ⫻ 10⫺9 A and Jei,
Jii 僆 [0, 1.67] ⫻ 10⫺9 A for current-based spiking networks
and the current-based mean-field model and Jee, Jie 僆 [0,1.67] ⫻ 10⫺9 S and Jei, Jii 僆 [0, 10] ⫻ 10⫺9 S for the
conductance-based spiking networks and the dynamic synapses were set to these values using Eq. 13.
As a verification of our current-based mean-field model(Eqs. 9–11) without dynamic synapses, we used the samesettings as in the networks of current-based spiking neurons(see Fig. 5G) for the mean-field model and compared theaverage firing rate. The mean-field model was in excellent
agreement with the current-based spiking network simulations
across the very wide parameter range sampled, with a mean
absolute error (MAE) of 0.8 Hz and a maximum error of 3.8 Hz
over the parameter range.
Interaction of synaptic dynamics with network activity
To give an idea of how dynamic synapses interact with the
network, we discuss an example in detail. The example in Figs. 6

and 7 shows a burst in a network of current-based spiking neurons
followed by a convergence to a steady-state firing rate (a conductance-based example would work nearly identically, except thatthe effect on the synaptic conductances would be harder tovisually interpret as a result of the effect of the driving force oneach synapse). We also show a sparse, random mean-field simulation using the same parameters. A network was chosen that
normally fires at 20 Hz in the absence of dynamic synapses, Jee ⫽Jie ⫽ 0.05 ⫻ 10⫺9 A and Jei ⫽ Jii ⫽ ⫺0.1 ⫻ 10⫺9 A. The
network was initialized with the dynamic synapse parameter
initial conditions set to a 5-Hz steady-state firing rate and thedesired firing rate for the network was set to 10 Hz, by adjusting
the scaling factors Amn in the Eq. 8 system. The firing rate initialconditions for the network were also set to 5 Hz, as best as
possible. Because a firing rate of 5 Hz gives an overall positive
network gain when the target firing rate is set to 10 Hz, the firing
rate shot up at the beginning of the simulation (Fig. 6, A and B).
The dynamic synapses adapted to this new firing rate and brought
it back down, ultimately reaching nearly the steady-state of 10 Hz
for the excitatory population, whereas the inhibitory population
firing rate remained near 20 Hz. Thus we see that an excitatory
firing rate very near the target firing rate x*e was reached fromquite different initial conditions and then maintained. During
the simulation, the average dynamic synapse weights weremodified and caused the changes shown in the network firing
rates (Fig. 6, C–J). We can assert that the dynamic synapses
caused the change in network firing rate (and not the otherway around) because the time constants of the dynamic
synapses are much longer than any neuronal time constants.
Again there is excellent agreement between the currentbased mean-field model and the current-based spiking network simulation.
A spike raster of a random subset of neurons in thecurrent-based spiking network example (Fig. 7A) shows thatthe network functioned in the asynchronous range, aside
from correlation induced by a common network firing rate.
Figure 7, B–D shows data from a single example neuron in

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

Exc. syn. strength (nA)

G 0.1

E 0.1
Exc. syn. strength (nA)

Exc. syn. strength (nA)

200

0.05 0.1 0.15

4

Network Correlation

Inh. syn. strength (nS)

Firing Rate (Hz)

0.08
0.06
0.04
0.02

0.5

0.5

Inh. syn. strength (nS)

D 0.1

1

Exc. syn. strength (nA)

1.5

C

CV

Exc. syn. strength (nS)

B

Firing Rate (Hz)Exc. syn. strength (nS)

Exc. syn. strength (nS)

A

4087

4088

D. SUSSILLO, T. TOYOIZUMI, AND W. MAASS

Firing Rate (Hz)

A

60

B 60

40

40

20

20

0
µee (nA)

C

100

200

300

400

500

0.08

µie (nA)

300

400

500

0

100

200

300

400

500

H 0.070

100

200

300

400

500

0

100

200

300

400

500

0

100

200
300
Time (msec)

400

500

0.04
0

100

200

300

400

500

0

100

200

300

400

500

0.12

F 0.15

0.1

G

0.1

0.06
0.05
0.04

0.06
0.05
0

µii (nA)

200

100

200

300

400

500

0.16
0.14
0.12
0.1
0.08
0.06

J
0.15
0.1
0.05

0

100

200
300
Time (msec)

400

500

the simulation. The total synaptic current hovered aroundthe background current and the recurrent synapses subtly
influenced this current. When the network achieved thesteady-state firing rate, the total current was slightly below
the normal background current (Fig. 7B), thus creating an
average firing rate of 10 rather than 20 Hz without dynamic
synapses. The individual synaptic currents in Fig. 7C show
differences in magnitude of three- to fourfold as a result of
the short-term plasticity.

Neuron #To investigate the self-tuning properties of dynamic synapses on spiking networks we asked whether dynamic synapses
could help tune a spiking neural network to a steady-statetarget firing rate even if the external input was perturbed
because it occurs, for example, in primary sensory areas of thecortex in response to sensory stimuli. To test this, we used a
conductance-based spiking network (Jee ⫽ Jie ⫽ 0.4 ⫻ 10⫺9 S

Circuit Response

A
200
400
0

B
Total Syn. Inputs(nA)

Self-tuning in networks of spiking neurons

3
2
1
0
−1

0

50

100

150

200

250

300

350

400

450

500

50

100

150

200

250

300

350

400

450

500

Syn. Inputs (nA)

C
0.1
0
−0.1
−0.2
−0.3
0

50

100

150

200

250

300

350

400

450

500

0

50

100

150

200

250
Time (msecs)

300

350

400

450

500

Vm (mV)

D
−40
−60
−80

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

FIG. 7. Detailed analysis of network activity in the network example from Fig. 6. A:
spike raster showing the spiking activity of
500 neurons (out of 5,000): 400 excitatory and
100 inhibitory. B–D: data from a single excitatory neuron. B: total synaptic excitatory andinhibitory currents. Blue line shows the total
excitatory input over time; red line shows thetotal inhibitory input over time; black line
shows the total amount of current into theneuron excluding the noise injection; green
dashed line shows the mean current injection
into the neuron, which corresponds to a 20-Hz
firing rate in the absence of synaptic feedback.
C: individual synaptic currents (each color
represents a separate synapse to the neuron).
Individual synaptic currents are plotted to givean idea of changes in synaptic weight values.
Upward exponential curves represent excitatory currents and downward, inhibitory currents. Peak values of these currents are set by
the dynamic synapses and vary roughly over
twofold the magnitude of the average current.
D: membrane potential of the selected neuron,with action potentials. Thick black line represents the membrane potential of the cell; thinblack line represents the average; blue line
shows the neuron’s firing threshold.

Downloaded from jn.physiology.org on June 12, 2007

100

0.08
0.07
0.06
0.05

0.08

I

0

D 0.0900.06

Eµei (nA)0

FIG. 6. Example showing the firing rates
and synaptic weights over time to demonstrate the effect of dynamic synapses on the
network firing rate. Scaling factors Amn were
set such that the desired steady-state firing
rate is 10 Hz, even though the synaptic
weight values [Je, Ji] ⫽ [0.5, ⫺1] ⫻ 10⫺10 A
gave a steady-state firing rate of 20 Hz without dynamic synapses. Initial conditions of
the dynamic synapses were set to values
corresponding to 5-Hz steady-state firingrate. A: firing rate curve of mean-field model:
xe in blue, xi in red. B: firing-rate curve of
network simulation; the firing rate issmoothed with a low-pass filter (box car filter
with length equivalent to 20 ms). Two firing
rate curves are in good qualitative agreement(the drop in firing rates at 500 ms is only afiltering artifact). C–J: synaptic weights
through time as set by the dynamic synapse
equations. C and D: ␮ee. E and F: ␮ei. G and
H: ␮ie. I and J: ␮ii. C, E, G, and I: mean-fieldmodel. D, F, H, and J: spiking networkmodel. Dashed lines represent the values Jee,Jei, Jie, and Jii. Comparison of the panels on
the left and right shows a good agreement
between the mean-field model and the spiking network.

SELF-TUNING OF NEURAL CIRCUITS

G

E
20

6

8

0

STD External Input (nA)

20

6

8

STD External Input (nA)

0

8

0

20

4

6

8

0

STD External Input (nA)

20

4

6

8

4

0

6

8

0

40

3.5
3
2.5
2
1.5

20

4

6

8

0

STD External Input (nA)

I
40

3.5
3
2.5
2
1.5

20

STD External Input (nA)

F
40

40

3.5
3
2.5
2
1.5

40

3.5
3
2.5
2
1.5

STD External Input (nA)

20

4

6

8

0

STD External Input (nA)

FIG. 8.

Self-tuning with respect to perturbations of the level and variance of external input currents. All networks tested were set to a default steady-state firing
rate of 10 Hz. Scaling factors Amn were set only once; thus each set of simulations (a single panel) where dynamic synapses were used had identical synapse
parameters. Blue denotes the optimal firing rate of 10 Hz (between 9 and 11 Hz), cyan between 8 and 9 and 11 and 12 Hz, green between 7 and 8 and 12 and
13 Hz, and gray denotes networks between 0 and 1 Hz. A–C: networks with conductance-based synapses. A: response (average excitatory firing rate) ofconductance-based networks with static synapses to perturbations of external current mean and noise. B: averaged response of conductance-based networks overdynamic synapse parameters R1–R3. C: response of conductance-based networks with the dynamic synapse parameters based on experimental data. D–F:
networks with current-based synapses. D: response of current-based networks with static synapses to perturbations of external current mean and noise strength.
E: averaged response of current-based networks over dynamic synapse parameters R1–R3. F: response of current-based networks with the dynamic synapse
parameters based on experimental data. G–I: response of current-based mean-field networks. G: response of current-based mean-field model with static synapses
to perturbations of external current mean and noise strength. H: averaged response of mean-field networks over dynamic synapse parameters R1–R3. I: responseof the mean-field networks with the parameters based on experimental data.

and Jei ⫽ Jii ⫽ 8.48 ⫻ 10⫺9 S), a current-based spiking
network (Jee ⫽ Jie ⫽ 0.013 ⫻ 10⫺9 A and Jei ⫽ Jii ⫽ ⫺0.18 ⫻10⫺9 A), and a current-based mean-field model (Jee ⫽ Jie ⫽0.013 ⫻ 10⫺9 A and Jei ⫽ Jii ⫽ ⫺0.18 ⫻ 10⫺9 A), all with asteady-state firing rate of 10 Hz. We additionally performed allof our self-tuning tests with the current-based mean-fieldmodel because it aids us in understanding the self-tuning
properties of networks with dynamic synapses and shows, inparticular, that the self-tuning properties of a circuit withdynamic synapses are a property of the dynamic synapse modeland not something hidden in the many details and parametersof a complicated spiking neural network simulation. Everysimulation was performed using both the parameters based on
experimental data and the R1–R3 dynamic synapse parametersthat were derived from our theoretical considerations (seeMETHODS). The scaling factors Amn of the dynamic synapses
were set only once such that the steady-state dynamic synapse
curves intersected at 10 Hz. Thus each simulation had identical
dynamic synapse parameters. We then perturbed in small increments both the mean and SD of the external input from ⫺50 to⫹50% of their normal values (Iinject ⫽ 2.455 ⫻ 10⫺9 A and noisedrawn from a 0-mean Gaussian, ␴inject ⫽ 6 ⫻ 10⫺9 A SD). The
results, which are the steady-state excitatory firing rates from thelast 1 s of each simulation, are shown in Fig. 8. All the networks
that had a steady-state firing rate between 9 and 11 Hz (adifference of ⱕ1 Hz from the target excitatory firing rate) are
shown in blue, networks with firing rates between 8 and 9 and 11
and 12 Hz (a difference between 1 and 2 Hz) are shown in cyan,

and networks with firing rates between 7 and 8 and 12 and 13 (adifference between 2 and 3 Hz) are shown in green. In Fig. 8, A,
D, and G we show the control case with static synapses for thethree different network models for comparison. There were few
networks that supported a 10-Hz steady-state firing rate. The
maximal firing rates were ⬎40 Hz and many networks were
essentially shut down with a firing rate very close to or equal to 0
Hz. Thus the steady-state firing rates of networks with staticsynapses are quite dependent on mean current injection andinjected noise level. When we ran the same simulations withdynamical synapses turned on (Fig. 8, remaining panels), the large
firing rates were reduced. We examined two different assignmentsof UDF values to the four different types of synaptic connectionsin the network: values that are sufficient for inducing stability
according to the previously (see first section of RESULTS) describedanalytical approach (Fig. 8, B, E, and H) and values taken fromexperimental data (Fig. 8, C, F, and I). Interestingly, for thisperturbation test the dynamic synapse parameters based on experimental data somewhat outperformed the dynamic synapses parameters based on the analytical approach.
We emphasize that when viewing changes in external input
as a perturbation to the network we are examining only thesteady-state response of the network and that there is indeed astrong transient response to external input that could provide
the excitation necessary for cortical computations, while stillallowing for homeostatic regulation on this very short timescale. A transient response to any network perturbation resultsfrom the fact that dynamic synapses do not equilibrate instan-J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

4

6

3.5
3
2.5
2
1.5

H
40

3.5
3
2.5
2
1.5

4

STD External Input (nA)

Mean External Input (nA)

0

40

3.5
3
2.5
2
1.5
4

Mean External Input (nA)8

20

Mean External Input (nA)

Mean External Input (nA)D

6Mean External Input (nA)4

STD External Input (nA)

C
40

3.5
3
2.5
2
1.5

Mean External Input (nA)20

Mean External Input (nA)

B
40

3.5
3
2.5
2
1.5

Mean External Input (nA)

Mean External Input (nA)

A

4089

4090

D. SUSSILLO, T. TOYOIZUMI, AND W. MAASS

Mean External Input (nS)

A

B
3.5
3
2.5
2
1.5

20

4

D
Mean External Input (nS)

40

6

0

8

40
20

4

6

3.5
3
2.5
2
1.5

E

3.5
3
2.5
2
1.5

0

8

C
40
20

4

6

40
20

4

6

Firing Rate (Hz)

10

200

400

600

800

4

6

0

8

40
20

4

6

0

8

STD External Input (nS)

30
20
10
0

0

200

Time (msec)

400

600

800

Time (msec)

taneously and thus the response will follow the timescale of the
dynamic synapses. This is shown in Fig. 9 for the case of
conductance-based spiking networks with dynamic synapse
parameters based on the experimental data shown in Table 1(Fig. 9, A–C) and R1 dynamic synapse parameters (Fig. 9,
D–F) (the R2 and R3 parameters were similar). Figure 9, A–F
shows 50-ms excitatory firing rate “snapshots” for each network at times 1–50, 251–300, and 451–500 ms. Additionally,
Fig. 9G shows a single network’s excitatory and inhibitory
response to a constant input with constant noise varianceexperimentally observed UDF parameters; Fig. 9H shows thesame network with R1 dynamic synapse parameters. In thecase of experimentally observed synapse parameters the network compensated for the changed external input in approximately 600 ms, resulting in an excitatory steady-state firing
rate near 11 Hz and an inhibitory steady-state firing rate of 24
Hz. In the case of the R1 dynamic synapse parameters the
network compensated for the changed external input in approximately 150 ms, resulting in an excitatory steady-state firing
rate and an inhibitory steady-state firing rate of approximately
24 Hz. This 150-ms period is in the same range as the durationof typical transient rate responses to visual stimuli in primaryvisual cortex (see, e.g., Fig. 2 in Lamme and Roelfsema 2000).
The duration of transient rate responses reported for higher
visual areas (see, e.g., Rainer et al. 2004) is somewhat longer
and lies between the durations found in these network simulations for R1 dynamic synapses and dynamic synapses based on
experimental data.
We next examined the effect of long-term changes in synaptic weights on firing rate stability in networks with dynamic
synapses. We used the parameter range of static synaptic
weight values already shown in Fig. 5 and assumed that these
synaptic weights were a perturbation of some unknown but
correct synaptic weights that would lead to a steady-state firing
rate of 10 Hz. Taking each network with its given, possibly
perturbed static synaptic weights, we then used Eq. 13 to set

the synaptic scaling factor for spiking networks with dynamic
synapses assuming a target firing rate of 10 Hz. That is, the Amnvalues were set to Amn ⫽ Jmn/[ R*mn (10 Hz)/U1*
mn (10 Hz)] andwere thus scaled by a single, constant scale factor acrosssimulations using the same dynamic synapse parameters.
The results of the synaptic weight perturbations for conductance-based spiking networks, current-based spiking networks,
and current-based mean-field models are shown in Fig. 10. In
Fig. 10, A, D, and G are the control cases with static synapses,
shown for comparison. With static synapses the steady-statefiring rates for networks with high excitatory synaptic weights
are in excess of 200 Hz. However, as shown in Fig. 10, B, C,
E, F, H, and I, in networks with dynamic synapses, regardlessof conductance-based (A–C) or current-based (D and E) synapses and regardless of which dynamic synapses parameters
are used, the excessive firing rates disappear. Most important,in all but the one case of current-based networks with the
dynamic synapse parameters based on experimental data, thenumber of perturbed networks that nevertheless had firing rates
very close to the target firing rate increased dramatically. The
continuous mean-field predictions of the dynamic synapses
applied to the networks with synaptic weight perturbations areshown in Fig. 10, G–I. The mean-field simulations show
excellent performance across nearly the entire range of perturbations and are qualitatively similar to spiking network simulations using dynamic synapses. Network correlation did notplay a role in the effect the dynamic synapses had on the
networks.
One may notice that some networks that had a firing rate
between 9 and 11 Hz without dynamic synapses no longer had
this firing rate when dynamical synapses were turned on(compare Fig. 10, D and E, bottom left corner). Furthermore,
the mean-field model with dynamic synapses (compare Fig. 10,
G and H) expands this range of target firing rates, whereas thecurrent-based spiking network does not. Typically, the average

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

Firing Rate (Hz)

20

0

0

8

H

30

20

3.5
3
2.5
2
1.5

STD External Input (nS)

G

40

F

3.5
3
2.5
2
1.5

STD External Input (nS)

0

0

8

3.5
3
2.5
2
1.5

FIG. 9. Transient firing rate responses in network models with dynamic synapses. A–C: transient firing rate responses for changes in the mean
and variance of external inputs are shown for conductance-based networks with experimentally observed UDF parameters for dynamic synapses (seeTable 1). A: firing rate of excitatory neurons during
the time interval 1–51 ms (in response to a change
in external inputs at time 0). Black circle indicatesthe specific network that is shown in G. Firing rateof excitatory neurons over periods of 251–300 and
450 –500 ms (B and C, respectively). D–F: samefor conductance-based networks using R1 dynamic
synapse parameters. G: average firing rate of excitatory (blue) and inhibitory (red) neurons as afunction of time (averaged with a 20-ms sliding
window) for the network represented in A as a
black dot. H: same for the network represented in D
as a black dot (the same network with static synapses assumed a steady-state firing rate of 19 Hz
for both excitatory and inhibitory neuron populations; not shown). One sees that in either case thetransient firing rate transmitted substantial information about the new level of external input during asignificant amount of time after the change in theexternal input, suggesting that a network-induced
homeostasis of firing rates is compatible with computational processes that process information about
the change in external inputs.

SELF-TUNING OF NEURAL CIRCUITS

100

0.5
4

6

8

1
0.5

0

2

Inh. syn. strength (nS)

200

0.08
0.06
0.04
0.02

100

0.05 0.1 0.15

2

0

8

0

200

0.08
0.06
0.04
0.02

100

0.05 0.1 0.15

0Inh. syn. strength (nA)

100

0.05 0.1 0.15

6

F 0.1

0

200

0.08
0.06
0.04
0.02

4

Inh. syn. strength (nS)

100

H 0.1
Exc. syn. strength (nA)

Exc. syn. strength (nA)

0200

0.05 0.1 0.15

100

0.5

I

0.1
0.08
0.06
0.04
0.02

0

Inh. syn. strength (nA)

200
100

0.05 0.1 0.15

0Inh. syn. strength (nA)

FIG. 10.

Self-tuning with respect to systematic synaptic strength perturbations. All networks tested were set to a default steady-state firing rate of 10 Hz. Blue
denotes the optimal firing rate of 10 Hz (between 9 and 11 Hz), cyan between 8 and 9 and 11 and 12 Hz, green between 7 and 8 and 12 and 13 Hz, and graydenotes networks between 0 and 1 Hz. Synaptic strength scaling factors Amn for all synaptic perturbations used the same rescaling factor. A–C: response (asmeasured by the steady-state average excitatory firing rate) of the conductance-based spiking networks to systematic synaptic strength perturbations. A: response
for conductance-based spiking networks with static synapses. B: response of conductance-based spiking networks with dynamic synapses using the R1–R3
parameters. C: response for conductance-based spiking networks with dynamic synapses using the parameters based on experimental data. D–F: response ofcurrent-based spiking networks to systematic synaptic strength perturbations. D: response of current-based networks with static synapses. E: response ofcurrent-based networks with dynamic synapses using the R1–R3 parameters. F: response of current-based networks with dynamic synapses using the parameters
based on experimental data. G–I: response of current-based mean-field firing rate models to systematic perturbations in synaptic strength. G: response ofcurrent-based mean-field models with static synapses. H: response of current-based mean-field models with dynamic synapses using the R1–R3 parameters. I:
response of current-based mean-field models with dynamic synapses using the parameters based on experimental data. Mean-field model differs from thecurrent-based network in F; see footnote 5 for discussion.

excitatory firing rate drops by 1 or 2 Hz when dynamic
synapses are turned on in the spiking network. The reason forthis is somewhat subtle and it has to do with the spiking
networks’ synaptic weight values, which are not always inexact agreement with the mean-field model synaptic strength
values5. Although this effect shows up in all three self-tuning
results we show, as is clear from the figures, its effect is small,resulting in differences between the spiking networks and themean-field models of only a few Hertz at most and, nevertheless, the results with dynamic synapses show great improvement over the control case with static synapses.
5

In a spiking network without dynamic synapses there is a distribution of
firing rates, whose mean we have used as the steady-state firing rate of the
network. In these networks even the tails of the neuronal firing rate distribution
have the same synaptic weight values Jmn. However, in the spiking networkwith dynamic synapses turned on, the lower part of the firing rate distribution
will typically have synapses that are stronger than the set value ␮mn(x*) ⫽ Jmn
for the target firing rate and the higher part of the distribution will typicallyhave synapses that are weaker than Jmn. Because the steady-state dynamic
synapse efficacy curves were constructed only to make sure the derivatives had
the right signs there is no reason that the expected value of the steady-state
synaptic weight across all synapses from neuron type n to m should exactly
equal the mean-field average [i.e., that 兰P(xn) ␮*mn (xn) ⫽ Jmn, where P(xn) is
the distribution of firing rates for neuron type n]. Furthermore, changing thetarget firing rate may change the shape of the steady-state firing rate distribution of the network, thus making it effectively impossible to cancel this effect
if one desires to tune the same network to different firing rates. Overall, thiseffect results in no more than a 1- to 2-Hz difference between the mean-field
prediction of the effect of dynamic synapses and the spiking network resultswith dynamic synapses.

Why should simply scaling the Amn for dynamic synapses
according to Eq. 13 have a helpful effect on the target firing
rate of the network if there is an unknown perturbation of thefactor Jmn in that same equation? Assume, for example, that
there are static synaptic strengths (J*e , J*i ) that lead to a target
firing rate of x*. Assume further that there is an unknown
perturbation of these synaptic strengths that can be written as(J*e ⫹ ␧e, J*i ⫺ ␧i). Because we have increased the excitation anddecreased the inhibition we should expect an increase in thetarget firing rate; that is, the perturbed synaptic strengths (J*e ⫹␧e, J*i ⫺ ␧i) lead to x* ⫹ x␧. If we now examine the same
network and synaptic weight perturbation but use dynamic
synapses we see a change in behavior. As already outlined we
set the Amn according to Eq. 13, so that ␮*me (x*) ⫽ J*e ⫹ ␧e and␮*mi (x*) ⫽ J*i ⫺ ␧i because we do not know that there is aperturbation. Now the situation with dynamic synapses is quite
different from with static synapses because at x* the dynamic
synapses equilibrate to values that support a higher firing rate
than x*, specifically x* ⫹ x␧. At a firing rate of x* ⫹ x␧,
however, because the derivatives of the dynamic synapses havethe correct signs, the network supports a firing rate below x* ⫹x␧. Thus we can conclude that the firing rate must be between
x* and x* ⫹ x␧, which is better than the case of static synapses,
where the network firing rate is completely dictated by the
perturbed synaptic weights and thus has a firing rate of x* ⫹ x␧.
The degree to which the final firing rate will be closer to x* or

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

Inh. syn. strength (nA)1

Inh. syn. strength (nA)

100

0.05 0.1 0.15

8

0.08
0.06
0.04
0.02

0

200

0.08
0.06
0.04
0.02

6

E 0.1

Inh. syn. strength (nA)

G 0.1

4

200

1.5

Inh. syn. strength (nS)

Exc. syn. strength (nA)

Exc. syn. strength (nA)

D 0.1

100

C

Exc. syn. strength (nA)

2

200

1.5

Exc. syn. strength (nA)

1

B

Exc. syn. strength (nS)

200

1.5

Exc. syn. strength (nS)

Exc. syn. strength (nS)

A

4091

4092

D. SUSSILLO, T. TOYOIZUMI, AND W. MAASS

0% Exc. Neurons Deleted0% Inh. Neurons Deleted

D

15

40

10

20

5

0

0

50

0

% Inh. Neurons Deleted

60

15

20

5
50
% Inh. Neurons Deleted

0

0

50

0

% Inh. Neurons Deleted

15

40

10

20

5
0

50

0

% Inh. Neurons Deleted

15

40

10

20

5
0

50

10

20

5

0

0

0

50

0

% Inh. Neurons Deleted

60

15

40

10

20

5

0

0

50

0

% Inh. Neurons Deleted

I

60

0

15

40

F

60

0

60

60

15

40

10

20

5

0

0

% Inh. Neurons Deleted

500% Inh. Neurons Deleted

FIG. 11.

Sparse, spiking networks demonstrate self-tuning with respect to systematic neuron inactivation. All networks tested were set to a default steady-statefiring rate of 10 Hz. Scaling factors Amn were set only once; thus each set of simulations (a single panel) where dynamic synapses were used had identical dynamic
synapse parameters. Blue denotes the optimal firing rate of 10 Hz (between 9 and 11 Hz), cyan between 8 and 9 and 11 and 12, green between 7 and 8 and 12
and 13 Hz, and gray denotes networks between 0 and 1 Hz. A–C: response (as measured by the steady-state average excitatory firing rate) of the
conductance-based spiking networks to systematic neuronal inactivation. A: response of conductance-based spiking networks with static synapses. B: averaged
response of conductance-based spiking networks with dynamic synapses using the R1–R3 parameters. C: response of conductance-based spiking networks withthe dynamic synapse parameters based on experimental data. D–F: response (as measured by the steady-state average excitatory firing rate) of the current-based
spiking networks to systemic neuronal inactivation. D: response of current-based spiking networks with static synapses. E: averaged response of current-based
spiking networks with dynamic synapses using the R1–R3 parameters. F: response of current-based spiking networks with dynamic synapses using the parameters
based on experimental data. G–I: response (as measured by the steady-state excitatory firing rate) of the current-based mean-field firing rate models to systemic
neuronal inactivation. G: response of current-based mean-field models with static synapses. H: averaged response of current-based mean-field models withdynamic synapses using the R1–R3 dynamic synapse parameters. I: response of current-based mean-field models with dynamic synapses using the parameters
based on experimental data.

x* ⫹ x␧ depends on the magnitude of these derivatives. We
refer to the last part of the Supplementary material for ananalytical comparison of network stability with static anddynamic synapses.
As a final self-tuning test we also investigated the effects
of systematic changes in the excitability of subpopulationsof neurons in the simulated circuit. Diverse changes in the
excitability of neurons were previously reported to becaused by neuromodulators (especially dopamine) in various cortical areas (Gulledge and Jaffe 1998, 2001; e.g., seeFig. 2B in first reference). Because such changes in the
excitability of neurons cannot be modeled directly in thecontext of LIF neurons, we simply inactivated (i.e., deleted)randomly selected subsets of neurons in the circuit. For
continuity, we chose the same networks as in the input
perturbation self-tuning simulation: conductance-based
spiking networks (Jee ⫽ Jie ⫽ 0.4 ⫻ 10⫺9 S and Jei ⫽ Jii ⫽8.48 ⫻ 10⫺9 S), current-based spiking networks (Jee ⫽ Jie ⫽0.013 ⫻ 10⫺9 A and Jei ⫽ Jii ⫽ ⫺0.18 ⫻ 10⫺9 A), and thecurrent-based mean-field models (Jee ⫽ Jie ⫽ 0.013 ⫻ 10⫺9
A and Jei ⫽ Jii ⫽ ⫺0.18 ⫻ 10⫺9 A), all with a steady-statefiring rate of 10 Hz. We randomly chose excitatory andinhibitory neurons in these networks to inactivate, with thetotal number of inactivated neurons between 0 and 70% of
their respective pool totals, in small increments. The cur-rent-based mean-field model received an analogous perturbation by modifying the Kmn parameter in Eq. 9 according to
the neuronal inactivation percentage in the spiking networks. The control simulations with static synapses areshown in Fig. 11 A, D, and G; results with the R1–R3
dynamic synapse parameters are shown in Fig. 11, B, E, andH; and results with the dynamic synapse parameters basedon experimental data are shown in Fig. 11, C, F, and I. For
the networks tested, the dynamic synapse parameters basedon experimental data did not have a significant effect but theR1–R3 dynamic synapse parameters showed improvement
over the control, with the vast majority of networks within
2 Hz of the target firing rate of 10 Hz. Again, network
correlation effects did not play a role in the results of thesimulations.
Choosing other target firing rates
Implicit in the demonstrated self-tuning properties is theability to set a specific fixed-point firing rate using dynamic
synapses, through use of the scaling in Eq. 13. The network
perturbations can then be viewed as a perturbation awayfrom this fixed point. We thus wondered, in principle, howclose a spiking network could get to a preprogrammed fixed
point, as set by Eq. 13, regardless of what the initial synaptic

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

10

0

0

H

40

0

5

E

60

G
% Exc. Neurons Deleted50% Exc. Neurons Deleted

0

10

20

% Exc. Neurons Deleted

5

15

40

% Exc. Neurons Deleted

10

20

C

60

% Exc. Neurons Deleted

1540% Exc. Neurons Deleted

B% Exc. Neurons Deleted% Exc. Neurons Deleted

A
60

SELF-TUNING OF NEURAL CIRCUITS

(Hz)

A

B

Mean Excitatory Firing Rate

120

0.6

100

0.5

80

0.4

60

0.3

40

0.2

20

0.1

0

0

20

40
60
80
100
Target Firing Rate (Hz)

120

0

4093

Excitatory Network Correlation

0

20

40
60
80
100
Target Firing Rate (Hz)120

strengths of the network were. Briefly, we used five networks from the current-based, spiking network synaptic
parameter range of Fig. 5D. The networks had Je and Ji as
follows (in nA): blue [0.025, ⫺0.05], red [0.025, ⫺0.15],
cyan [0.05, ⫺0.1], magenta [0.075, ⫺0.05], and green[0.075, ⫺0.015], leading to corresponding firing rates of 20,
12, 20, 100, and 20 Hz with static synapses. We then choseas the target firing rates every integral firing rate in the range
of 1–120 Hz and used Eq. 13 to appropriately set the
dynamic synapses in the five networks. We do not mean tosuggest that cortical networks create fixed points at suchhigh firing rates, rather only to test the principle behind Eq.
13. The results, which are the steady-state excitatory firing
rates, are shown in Fig. 12A. All networks showed an ability
to enhance their firing rates toward the desired target firing
rate, although the two networks with the weakest excitatory
synapses (denoted by blue and red) were unable to follow
the target firing rate closely at ⬎20 Hz. Another complication arose in that the two networks with the strongest
excitatory synapses (green and magenta), although able toclosely follow the target firing rate, also showed significant
network correlation at ⬎20 Hz, shown in Fig. 12B. Despite
these caveats, the results clearly demonstrate that the ability
to tune the steady-state firing rate of a network is a generalproperty of the dynamic synapse configuration. The results
lend further evidence to the idea that self-tuning propertiesof spiking networks with dynamic synapses are caused by
setting a firing rate fixed point, to which the networks
attempt to return after a perturbation.
Theoretical analysis
There is little hope that one can give a complete theoreticaljustification for the demonstrated self-tuning properties ofsparsely connected circuits of spiking neurons with dynamic
synapses or for the associated mean-field model defined by
Eqs. 7–9. The reason lies in the nonlinearity of these models
and the fairly large number of variables that are involved.

However, the demonstrated self-tuning properties of thesecomplex systems are clearly related to known self-tuning
properties of substantially simpler systems for which a theoretical analysis is feasible. These analytical results will be
discussed in the following.
An analytical result on the assignment of differential
dynamics to different types of synapses in a
distributed circuit
In the first section of RESULTS we proposed a heuristic forchoosing the signs of the derivatives of the synaptic weight
curves ␮*mn (xn) at x*. We then demonstrated through numericalsimulations of circuits of spiking neurons that this choice wasjustified. This heuristic can also be justified analytically for a
2D mean-field model, similar to that described in Eq. 11(assuming dynamic synapses with very small derivatives andinstantaneously equilibrating dynamic synapses). One canprove that the derivatives dmn ⫽ (d␮mn/dxn)(x*) should havethe following signs for the excitatory population to havehomeostatic properties with respect to changes in the intensityof external input currents: dee ⬍ 0, dii ⬎ 0, die ⬎ 0, dei ⬍ 0 (seeSupplementary material for details). These signs agree with thepreviously proposed heuristic rules for the assignment of depressing and facilitating synapses to the four types of synaptic
connections (see first section of RESULTS; note that the negative
sign is included in ␮ei and ␮ii). For the excitatory population toalso have homeostatic properties with respect to changes in theintensity of external input currents to the inhibitory population
the analytically derived condition is somewhat more complicatedd ce ⬍ 0, d ⬎ 0, d ic ⬎ 0, d ei ⬍ 0

d ee ⬍ 0, d ii ⬍ 0, d ie ⬎ 0, d ei ⬎ 0

if ␮ ee 共x*兲 ⬎1␤eif ␮ ce 共x*兲 ⬍1␤ewhere ␤e is the derivative of the excitatory F–I curve evaluatedJ Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

FIG. 12. Many spiking networks can follow a target firing rate using dynamic synapses. Five selected current-based spiking networks with dynamic synapses
showing the networks’ abilities to follow an excitatory target firing rate using scaling factors Amn as set by the heuristic of Eq. 13. Target firing rate was set toevery integral firing rate between 1 and 120 Hz and the A values for each network were calculated based on this target frequency. A: average target firing rate;
black dashed line shows the line y ⫽ x, which is optimal performance. Mean absolute errors for each network blue, red, cyan, magenta, and green are [20.7, 21.4,
4.5, 2.1, 2.0] Hz, respectively. Networks had Je and Ji as follows (in nA): blue [0.025, ⫺0.05], red [0.025, ⫺0.15], cyan [0.05, ⫺0.1], magenta [0.075, ⫺0.05],
and green [0.075, ⫺0.015] and are also denoted in Fig. 5D. Dynamic synapse parameters used were the R1 parameters. B: excitatory network correlation as afunction of target firing rate. There was negligible network correlation for the blue, red, and cyan networks and significant network correlation for the magenta
and green networks between the firing rates of 20 and 100 Hz.

4094

D. SUSSILLO, T. TOYOIZUMI, AND W. MAASS

at the fixed point x*. Taken together these conditions imply
that regardless of the size of ␮ee(x*), the E 3 E connection
should always be depressing and the E 3 I should always be
facilitating to allow the circuit to return to a steady-state firing
rate in spite of a changed level of external input. Additionally,
if ␮ee(x*) ⬎ ␤⫺1e the I 3 E synapses should be facilitating andthe I 3 I synapses should be depressing; otherwise, the I 3 E
synapses should be depressing and the I 3 I synapses should
be facilitating. The fact that this analytical result does not point
to a general rule for the optimal dynamics of synapses from
inhibitory neurons (from the perspective of firing rate stability)
could be seen as a possible explanation for the diverse distribution of experimentally observed UDF parameters for these
synapses (see Fig. 4, E and F). We additionally give justification for our self-tuning heuristic for the parameter range usedin our current-based networks in the Supplementary materials.

One can prove for the same mean-field model with excitatory and inhibitory neurons that the steady-state firing rate x*
is less affected by small changes in synaptic weights (resulting
for example from LTP or LTD) if the synaptic connections are
dynamic and tuned using an appropriate heuristic. An analytical derivation is provided in the second part of the Supplementary material, with specific reference to the fact that our
self-tuning heuristic is appropriate for the parameter regime
used in the current-based spiking networks herein. An additional example figure (Fig. S1) is also given.
DISCUSSION

We have shown that the inclusion of more realistic models
for synapses, which reflect their experimentally found shortterm plasticity, provides models for cortical neural networks
with interesting stability properties. In particular, it enables
such models to return after a large variety of perturbations to a
low but nonzero spontaneous firing rate, even if the perturbation, such as a change in the level of external input, a change
in synaptic strength, or even neuron deletion, is of a longerlasting nature. This surprising stability property of networkmodels with dynamic synapses, which has apparently not been
reproduced in models with static synapses, may provide a
possible explanation for the surprisingly stable low but nonzero
rate of spontaneous firing reported for cortical neurons in avariety of studies. Several other homeostatic processes, such as
synaptic scaling, LTP and LTD, and genetic regulation ofreceptor numbers are likely to support the stability of the rateof spontaneous firing of cortical neurons on a larger timescale(Turrigiano and Nelson 2004). However, none of these mechanisms would apparently be able to support stability on thetimescale of seconds, say, in response to a changed level ofsensory input or to a changed concentration of neuromodulators in a cortical circuit. On the other hand, we showed in Fig.
9 that the stability property endowed by short-term synaptic
plasticity leaves a sufficiently large time window of one toseveral hundred milliseconds during which the firing rate in thecircuit is affected by the external perturbation. Thus information about such external inputs can be transmitted to other

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

A theoretical analysis of the impact of changes of synaptic
weights on stability properties of networks with static anddynamic synapses

cortical circuits and integrated into cortical computation. In
fact, the dynamic response to an increased external input
shown in Fig. 9H is on a timescale similar to that of typicalresponses to visual inputs in primary visual cortex (see, e.g.,
Fig. 2 in Lamme and Roelfsema 2000). Even on a larger
timescale, however, the information about the external (orinternal) perturbation remains accessible to the neural system,
although the firing rate of excitatory neurons is returned
through synaptic dynamics to a given target level. The samefiring rate of excitatory neurons is then accompanied by adifferent firing rate of inhibitory neurons (see Fig. 9, G and H)
and this rate may therefore contain substantial information
about the nature and level of the external perturbation for amuch longer time period. This observation also implies that thesame firing rate of excitatory neurons can be produced by avirtually infinite combination of concentration levels of neuromodulators and activity levels of inhibitory neurons, analogously to the large sets of different circuit parameters that werereported to support the generation of a target firing pattern inthe pyloric circuit in the crustacean stomatogastric ganglion(Marder and Goaillard 2006).
The target value of 10 Hz for the spontaneous firing rate
excitatory neurons in cortical circuits on which we have focused in our network simulations was chosen rather arbitrarily,
although it lies within a range of spontaneous firing rates
reported in a number of experimental studies (see INTRODUCTION). We show in Fig. 12 that a large range of other target
values can be achieved by choosing a suitable global scaling
parameter for synaptic weights, which might for example beunder genetic control. However, to demonstrate stability atreally low rates (say, 1 Hz), one needs to add neuron models to
the network that fire spontaneously even in the absence ofsynaptic input (or during a release from inhibition); otherwise,
the network cannot recover from a temporary pause in firing
and will remain indefinitely silent.
We have shown in this article that the stability of the firing
rate of excitatory neurons endowed by short-term synaptic
dynamics is rooted in general principles of control theory,
recently proposed for much simpler dynamical systems. In fact,
activity-dependent rescaling of feedback strength was previously postulated as a powerful mechanism for self-tuning incontrol theory for purely mathematical reasons (Moreau andSontag 2003). We demonstrated that this theoretical analysis ofstability properties provided by dynamic synapses can beextended to the case of a distributed system that consists ofseveral interacting dynamical modules, such as an interacting
population of excitatory and inhibitory neurons in a corticalcircuit, each with thousands of neurons. In particular, we haveexhibited an analytical method for deriving rules regarding thetype of synaptic dynamics that should be assigned to different
types of synapses for the purpose of network stability. This
theoretical approach may help to provide an explanation for thediversity of the empirically found short-term plasticity ofconnections between modules with different dynamic roles(such as excitatory and inhibitory neurons). It turns out that thistheoretical analysis provides clearer predictions for the dynamics of synapses from excitatory neurons than for GABAergic
synapses. This might provide an explanation for the experimentally found diversity in the dynamics of GABAergic synapses.

SELF-TUNING OF NEURAL CIRCUITS

Finally, we emphasize that the investigation of the relationship between stability properties of networks of spiking neurons and the dynamics of specific types of synapses that wehave initiated in this article should be seen only as a first step.
In a second step one also needs to investigate through computer
simulations more detailed network models consisting of different types of excitatory and especially inhibitory neurons
with specific firing properties and more specific connectivity
patterns. A quite challenging open problem is whether one canalso extend the theoretical analysis of stability properties induced by dynamic synapses to more complex networks.
ACKNOWLEDGMENTS

The authors thank Drs. Stephan Haeusler and Eduardo Sontag for fruitful
discussions that greatly aided this paper and Dr. Rodney Douglas for very
helpful critical comments.
GRANTSREFERENCES

Abbott LF, Nelson SB. Temporal dynamics of biological synapses. In: The
Handbook of Brain Theory and Neural Networks, edited by Arbib MA.
Cambridge, MA: The MIT Press, 2003, p. 1156 –1158.
Abbott LF, Regehr WG. Synaptic computation. Nature 431: 796 – 803, 2004.
Abbott LF, Varela JA, Sen K, Nelson SB. Synaptic depression and cortical
gain control. Science 275: 220 –224, 1997.
Abeles M. Corticonics: Neural Circuits of the Cerebral Cortex. Cambridge,
UK: Cambridge Univ. Press, 1991.
Baddeley R, Abbott LF, Booth MC, Sengpiel F, Freeman T, Wakeman
EA, Rolls ET. Responses of neurons in primary and inferior temporal visual
cortices to natural scenes. Proc Biol Sci 264: 1775–1783, 1997.
Brunel N. Dynamics of sparsely connected networks of excitatory and inhibitory spiking neurons. J Comput Neurosci 8: 183–208, 2000.
Buzsáki G. Rhythms of the Brain. New York: Oxford Univ. Press, 2006.
Destexhe A, Contreras D. Neuronal computations with stochastic network
states. Science 314: 85–90, 2006.
Destexhe A, Rudolph M, Paré D. The high-conductance state of neocorticalneurons in vivo. Nat Rev 4: 739 –751, 2003.
Galarreta M, Hestrin S. Frequency-dependent synaptic depression and thebalance of excitation and inhibition in the neocortex. Nat Neurosci 1:
587–594, 1998.
Gulledge AT, Jaffe DB. Dopamine decreases the excitability of layer V
pyramidal cells in the rat prefrontal cortex. J Neurosci 18: 9139 –9151, 1998.
Gulledge AT, Jaffe DB. Multiple effects of dopamine on layer V pyramidal
cell excitability in rat prefrontal cortex. J Neurophysiol 86: 586 –595, 2001.
Gupta A, Wang Y, Markram H. Organizing principles for a diversity ofGABAergic interneurons and synapses in the neocortex. Science 287:
273–278, 2000.
Haeusler S, Maass W. A statistical analysis of information-processing properties of lamina-specific cortical microcircuit models. Cereb Cortex 17:
149 –162, 2007.

Lamme VA, Roelfsema PR. The distinct modes of vision offered by feedforward and recurrent processing. Trends Neurosci 23: 571–579, 2000.
Legenstein R, Maass W. Edge of chaos and prediction of computational
performance for neural microcircuit models. Neural Networks In press.
Lisman JE. Bursts as a unit of neural information: making unreliable synapses
reliable. Trends Neurosci 20: 38 – 43, 1997.
Maass W, Markram H. Synapses as dynamic memory buffers. Neural
Networks 15: 155–161, 2002.
Maass W, Natschlager T, Markram H. Fading memory and kernel propertiesof generic cortical microcircuit models. J Physiol (Paris) 98: 315–330,
2004.
Marder E, Goaillard JM. Variability, compensation and homeostasis inneuron and network function. Nat Rev 7: 563–574, 2006.
Markram H, Wang Y, Tsodyks M. Differential signaling via the same axonof neocortical pyramidal neurons. Proc Natl Acad Sci USA 95: 5323–5328,
1998.
Moreau L, Sontag E. Balancing at the border of instability. Phys Rev 68:
020901, 2003a.
Moreau L, Sontag E. Feedback tuning of bifurcations. Syst Control Lett 50:
229 –239, 2003b.
Rainer G, Lee H, Logothetis NK. The effect of learning on the function ofmonkey extrastriate visual cortex. PLoS Biol 2: e44, 2004.
Shriki O, Hansel D, Sompolinsky H. Rate models for conductance-based
cortical neuronal networks. Neural Comput 15: 1809 –1841, 2003.
Steriade M. The Intact and Sliced Brain. Cambridge, MA: The MIT Press,
2001.
Steriade M, Timofeev I, Grenier F. Natural waking and sleep states: a view
from inside neocortical neurons. J Neurophysiol 85: 1969 –1985, 2001.
Thomson AM. Activity-dependent properties of synaptic transmission at twoclasses of connections made by rat neocortical pyramidal axons in vitro.
J Physiol 502: 131–147, 1997.
Thomson AM. Presynaptic frequency- and pattern-dependent filtering. J Comput Neurosci 15: 159 –202, 2003.
Tsodyks M, Pawelzik K, Markram H. Neural networks with dynamic
synapses. Neural Comput 10: 821– 835, 1998.
Turrigiano GG, Nelson SB. Homeostatic plasticity in the developing nervous
system. Nat Rev 5: 97–107, 2004.
van Vreeswijk C, Sompolinsky H. Chaos in neuronal networks with balanced
excitatory and inhibitory activity. Science 274: 1724 –1726, 1996.
van Vreeswijk C, Sompolinsky H. Chaotic balanced state in a model ofcortical circuits. Neural Comput 10: 1321–1371, 1998.
Varela JA, Sen K, Gibson J, Fost J, Abbott LF, Nelson SB. A quantitativedescription of short-term plasticity at excitatory synapses in layer 2/3 of rat
primary visual cortex. J Neurosci 17: 7926 –7940, 1997.
Vogels TP, Abbott LF. Signal propagation and logic gating in networks ofintegrate-and-fire neurons. J Neurosci 25: 10786 –10795, 2005.
Vogels TP, Rajan K, Abbott LF. Neural network dynamics. Annu Rev
Neurosci 28: 357–376, 2005.
Wang Y, Markram H, Goodman PH, Berger TK, Ma J, Goldman-Rakic
PS. Heterogeneity in the pyramidal network of the medial prefrontal cortex.
Nat Neurosci 9: 534 –542, 2006.
Yoshimura Y, Callaway EM. Fine-scale specificity of cortical networks
depends on inhibitory cell type and connectivity. Nat Neurosci 8: 1552–
1559, 2005.
Yoshimura Y, Dantzker JL, Callaway EM. Excitatory cortical neurons form
fine-scale functional networks. Nature 433: 868 – 873, 2005.

J Neurophysiol • VOL 97 • JUNE 2007 • www.jn.org

Downloaded from jn.physiology.org on June 12, 2007

D. Sussillo was supported in part by a Fulbright Research Grant from the
U.S. State Department. W. Maass was partially supported by the Austrian
Science Fund (Projects P17229 and S9102-N13) and by the European Union
Project Fast Analog Computing with Emergent Transient States (FP6-015879).

4095