LETTER
Communicated by Bard Ermentrout
Firing Rate of the Noisy Quadratic Integrate-and-Fire Neuron
Nicolas Brunel
brunel@biomedicale.univ -paris5.fr
CNRS, NPSM, Universit´e Pa ris Ren´e Des cartes, 75270 Paris Cedex 06, France
Peter E. Latham
pel@ucla.edu
Department of Neurobiology, University of California at Los Angeles,
Los Ang eles, CA 90095, U.S.A.
We calculate the ring rate of the quadratic integrate-and-re neuron in
response to a colored noise input current. Such an input current is a good
approximation to the noise due to the random bombardment of spikes,
with the correlation time of the noise corresponding to the decay time
of the synapses. The key parameter that determines the ring rate is the
ratio of the correlation time of the colored noise,
¿
s
, to the neuronal time
constant,
¿
m
. We calculate the ring rate exactly in two limits: wh en the
ratio,
¿
s
=¿
m
, goes to zero (white noise) and when it goes to innity. The
correction to the short correlation time limit is O
.¿
s
=¿
m
/
, which is qualita-
tively different from that of the leaky integrate-and-re neuron, where the
correction is O
.
p
¿
s
=¿
m
/
. The difference is due to the different boundary
conditions of the probability density function of the membrane potential
of the neuron at ring threshold. The correction to the long correlation
time limit is O
.¿
m
=¿
s
/
. By combining the short and long correlation time
limits, we derive an expression that provides a good approximation to the
ring rate over the whole range of
¿
s
=¿
m
in the suprathreshold regime—
that is, in a regime in which the average current is sufcient to make the
cell re. In the subthreshold regime, the expression breaks down some-
what when
¿
s
becomes large compared to
¿
m
.
1 Introduction
A major challenge in computational neuroscience is to understand the be-
havior of large recurrent networks of spiking neurons. A rst step in this
endeavor is to be able to compute the location and stability of a network’s
equilibria, given connectivity and single neuron properties. For networks
operating in the asynchronous regime, the equilibria can be determined by
solving a set of algebraic equations in which the ring rate of each neuron
in the network is a function of the ring rates of all the other neurons. To
solve these equations, which can be done using mean-eld techniques, it is
Neural Computation
15, 2281–2306
(2003)
c
° 2003 Massachusetts Institute of Technology
2282 N. Brunel and P. Latham
necessary to be able to compute the mapping from presynaptic ring rates
to postsynaptic rates. For leaky integrate-and-re neurons, the mapping is
known (Ricciardi,1977; Tuckwell, 1988), and equilibria have been co mputed
for networks of irregularly ring leaky integr ate-and-re neurons (Amit &
Brunel, 1997a, 1997b; Brunel, 2000). Here we compute this mapping for a
more realistic reduced neuron model, the quadratic integrate-and-re neu-
ron (Ermentrout & Kopell, 1986).
Cortical neurons in vivo re in an irregular fashion (Burns & Webb, 1976;
Softky & Koch, 1993), and intracellular recordings reveal large uctuations
in the membrane potential (Destexhe & Par´e, 1999; Anderson, Lampl, Gille-
spie, & Ferster, 2000). These experimental observations suggest that the
input to a neuron can be divided into two components, a mean current and
uctuations around that mean, both of which depend on the ring rates of
the presynaptic neurons. It is often reasonable to approximate the uctu-
ations as colored noise with a correlation time proportional to the synap-
tic time constants (Amit & Tsodyks, 1991; Brunel & Sergi, 1998; Destexhe,
Rudolph, Fellous, & Sejnowski, 2001), and this is the approach we take here.
In addition, we assume that there is a unique synaptic d ecay time,
¿
s
, so that
the correlation time of the noise is equal to
¿
s
. With these assumptions, the
input of the neuron is full y described by three parameters: the mean current,
¹
; the variance of the uctuations,
¾
2
; and the correlation time constant,
¿
s
.
The ring rate of the leaky integrate-and-re neuron versus
¹
and
¾
was
calculated by Brunel and Sergi (1998; see also Fourcaud & Brunel, 2002),
in the limit
¿
s
¿
¿
m
where
¿
m
is the membrane time constant. Although
the leaky integrate-and-re neuron is a popular model, it suffers from two
drawbacks: it cannot be derived from conductance-based models using a
rigorous reduction procedure, and its
f
-I curve in the absence of noise has
a pathological behavior close to threshold. For these reasons, it is useful
to consider another simplied model neuron: the quadratic integrate-and-
re neuron. This neuron, which is related by a change of variables to the
µ
-neuron (Ermentrout, 1996; Gutkin & Ermentrout, 1998) represents the
normal form of a neuron with a type I bifurcation leading to spike generation
(Ermentrout & Kopell, 1986; Ermentrout, 1996). The quadratic integrate-
and-re neuron is thus expected to describe the dynamics of any type I
neuron close to bifurcation, where ring rates are low. In particular, its ring
rate in the absence of noise scales as
p
I
¡
I
thr
where
I
thr
is the threshold
current. This is the behavior exhibited by all type I neurons near threshold.
Here we compute the ring rate of the quadratic integrate-and-re neu-
ron in two limits: small and large
¿
s
=¿
m
. Interpolation between these limits
provides a good approximation to the ring rate o ver the whole range of
¿
s
=¿
m
in the suprathreshold regime and a reasonable zeroth-order approx-
imation in the subthreshold regime. The ring rate must be computed nu-
merically, but for applications such as mean-eld analysis, where speed is
important, two lookup tables are sufcient to characterize the ring rate for
all values of
¹
,
¾
,
¿
m
, and
¿
s
.
Firing Rate of the Noisy Quadratic Integrate-and-Fire Neuron 2283
2 Neuron Model and Analytical Methods
The quadratic integrate-and-re neuron with synaptic noise obeys the dif-
ferential equations
¿
m
dv
dt
D
f .v/
C
w
(2.1a)
¿
s
dw
dt
D ¡
w
C
¿
1
=
2
m
¾ ´.t/
(2.1b)
where
¿
m
is the neuronal time constant (hereafter called the membrane time
constant),
v.t/
is a “voltage” variable (hereafter called the membrane poten-
tial),
f .v/
D
v
2
C
¹
(2.2)
with
¹
the mean input current owing into the cell,
w
represents the stochas-
tic component in the synaptic cur rents,
¿
s
is the synaptic time constant,
´
is
white noise with h
´.t/´.t
0
/
i D
±.t
¡
t
0
/
, and
¾
is the overall amplitude of the
noise.
In the noiseless version of the model, when the membrane potential
becomes positive, the quadratic term on the right-hand side of equation 2.1a
ensures that the membrane potential diverges to innity in nite time. The
time at which the membrane potential reaches innity denes the time
a spike is emitted. Another way of dening the spike time would be to
set a threshold for ring and then formally send this threshold to innity.
Immediately after ring, the membrane potential is reset to ¡1. Again,
the quadratic term in the right-hand side of equation 2.1a ensures that the
membrane potential comes back to a nite value in nite time. To determine
the ringrate from equation 2.1, we reformulate the problem using a Fokker-
Planck equation (Risken, 1984),
¿
m
@
t
P.v; w; t/
C
@
v
[
. f.v/
C
w/P
] C
¿
m
¿
s
@
w
[¡
wP
] D
¾
2
2
¿
2
m
¿
2
s
@
2
w
P;
where
P.v; w; t/
is the joint probability density function (pdf) of
.v; w/
at
time
t
. Once
P.v; w; t/
is known, the ring rate is given by the total proba-
bility ux at
v
D C1 (see equation 2.5 below). For other studies of the pdf
approach to neural modeli ng, mostly in the contex t of integrate-and-re
neurons, see Treves (1993), Abbott and van Vreeswijk (1993), Brunel and
Hakim (1999), Brunel (2000), Knight, Omurtag, and Sirovich (2000), and
Nykamp and Tranchina (2000).
The Fokker-Planck equation can be rewritten as a continuity equation,
@
t
P.v; w; t/
C
@
v
J
v
.v; w; t/
C
@
w
J
w
.v; w; t/
D 0
;
(2.3)
2284 N. Brunel and P. Latham
where
J
v
and
J
w
are probability uxes,
J
v
.v; w; t/
D
1
¿
m
. f .v/
C
w/P.v; w; t/
(2.4a)
J
w
.v; w; t/
D
1
¿
s
³
¡
¿
m
¿
s
¾
2
2
@
w
P.v; w; t/
¡
wP.v; w; t/
´
:
(2.4b)
The ring rate,
º.t/
, is the total probability ux at
v
D C1; it is given by
º. t/
D
Z
C1
¡1
dw
lim
v
!1
J
v
.v; w; t/:
(2.5)
We are interested in computi ng the ring rate in the limit
t
! 1. In
this limit, the probability distribution becomes time independent. We thus
set
@
t
P
to zero in equation 2.3, resulting in a continuity equation involving
only
v
and
w
. That equation cannot be solved exactly, but it can be solved
as an expansion in the ratio of the synaptic to the membrane time constant
(or its inverse). To facilitate this expansion, we introduced a new variable,
z
´
kw=¾
, where
k
´
.¿
s
=¿
m
/
1
=
2
. The advantage of using
z
rather than
w
is that its variance remains O
.
1
/
in both the large and small
k
limits, a
property not shared by
w
. With this change of variable and with
@
t
P
set to
zero, equation 2.3 may be written as
LP.v; z/
D
k¾ z@
v
P.v; z/
C
k
2
@
v
[
f .v/P.v; z/
] (2.6)
where
L
¢ ´
1
2
@
2
z
¢ C
@
z
.z
¢
/:
In the next two sections, we solve this equation in the short and long corre-
lation time limits.
3 Short Correlation Time Limit
To solve equation 2.6 in the shor t correlation time limit (
k
¿ 1), we expand
both the probability distribution and the r ing rate in powers of
k
,
P.v; z/
D
1
X
i
D0
k
i
P
i
.v; z/
(3.1a)
º
D
1
X
i
D0
k
i
º
iS
;
(3.1b)
where the subscript
S
on
º
iS
is to remind us that we are working in the short
correlation time limit. (Note that
P
i
should also have a subscript
S
; we do
Firing Rate of the Noisy Quadratic Integrate-and-Fire Neuron 2285
not include it for clarity. It should be clear fro m the context whether we are
referring to the short or long correlation time limit; the latter is discussed
in section 4.) Since
P.v; z/
is a probability distribution, it must integrate to
one, independent of
k
. Consequently, we must have
Z
dv dz P
0
.v; z/
D 1 (3.2a)
Z
dv dz P
i
.v; z/
D 0
;
for
i
¸ 1
:
(3.2b)
To derive equations for the
P
i
, we insert equation 3.1a into equation 2.6
and match powers of
k
. This results in a coupled set of equations:
LP
0
.v; z/
D 0 (3.3a)
LP
1
.v; z/
D
¾ z@
v
P
0
.v; z/
(3.3b)
LP
i
.v; z/
D
¾ z@
v
P
i
¡1
.v; z/
C
@
v
[
f .v/P
i
¡2
.v; z/
]
;
for
i
¸ 2
:
(3.3c)
Using equations 2.4a, 2.5, and 3.1b, the ring rate at order
i
is
º
iS
D
1
¿
m
Z
C1
¡1
dz
lim
v
!1
[
f .v/P
i
.v; z/
C
¾ zP
i
C1
.v; z/
]
D
1
¿
m
Z
C1
¡1
dz
lim
v
!1
f .v/P
i
.v; z/;
(3.4)
where the last equality follows because
P.v; z/
must vanish as
v
! 1;
otherwise, it will not be nor malizable.
3.1 Probability Distribution in the Short Correlation Time Limit. To com-
pute the probability distribution as an expansion in powers of
k
, we simply
solve equations 3.3a through 3.3c order by order. The zeroth-order solu-
tion,
P
0
, corresponds to white noise; the next orders give corrections asso-
ciated with the nite synaptic time constant. It turns out, as we will see,
that we have to go to the fourth order to determine the lo west nonvan-
ishing correction to the ring rate. In this section, however, we sketch the
solutions to equations 3.3a through 3.3c only through second order; the
solutions through fourth order are derived in appendix A. We begin with
equation 3.3a.
Examining equation 3.3a, we see that the zeroth-order contributi on to the
probability distribution,
P
0
, can be written as a linear combination of two
functions: one even in
z
and the other o dd. Those two functions, denoted
Á
0
and
Á
1
, are given by
Á
0
.z/
D
¼
¡1
=
2
exp
.
¡
z
2
/
Á
1
.z/
D exp
.
¡
z
2
/
Z
z
0
exp
.u
2
/ du:
