In a population of neurons, each neuron may be in a different internal state.
In this section we derive partial differential equations that describe how the
distribution of internal states evolves as a function of time. We start in
Section 6.2.1 with a population of integrate-and-fire
neurons. Since the state of an integrate-and-fire neuron is characterized by
its membrane potential, we describe the dynamics of the population as the
evolution of membrane potential densities. In Section 6.2.2 we
turn to neurons described by the Spike Response Model. SRM neurons are
characterized by their state of refractoriness so that we have to introduce
refractory densities. We will see that the solution of the dynamics of
refractory densities leads to a macroscopic population activity equation for
escape noise models. In Section 6.2.3 we compare the two
approaches (i.e. membrane potential densities and refractory densities) and
show that they are equivalent in the low-noise limit. The formulation of the
dynamics of a population of integrate-and-fire neurons on the level of
membrane potential densities has been developed by Abbott and van Vreeswijk (1993),
Brunel and Hakim (1999), Fusi and Mattia (1999), Nykamp and Tranchina (2000), and Omurtag et al. (2000).
The closely related formulation in terms of refractory densities has been
studied by Wilson and Cowan (1972), Gerstner and van Hemmen (1992), Bauer and Pawelzik (1993), and
Gerstner and van Hemmen (1994). Generalized density equations have been discussed by
6.2.1 Integrate-and-Fire Neurons with Stochastic Spike Arrival
We study a homogeneous population of integrate-and-fire neurons. The internal
state of a neuron i is determined by its membrane potential which changes
ui = - ui + RIi(t) .
Here R is the input resistance,
= RC the membrane time constant,
and I(t) the total input (external driving current and synaptic input).
ui = the membrane potential is reset to
ui = ur < .
In a population of N integrate-and-fire neurons, we may ask how many of the
neurons have at time t a given membrane potential. For
fraction of neurons i with membrane potential
u0 < ui(t)u0 + u is
= p(u, t) du
where p(u, t) is the membrane potential density;
cf. Chapter 5.5. The aim of this section is to describe the
evolution of the density p(u, t) as a function of time. As we will see, the
equation that describes the dynamics of p(u, t) is nearly identical to that
of a single integrate-and-fire neuron with diffusive noise;
cf. Eqs. (5.88) and (5.89).
There are three subtle differences though. First, while p(u, t) was
introduced in Chapter 5.5 as probability density for
the membrane potential of a single neuron, it is now interpreted as the
density of membrane potentials in a large population of neurons.
Second, the normalization is different. In Chapter 5.5 the
p(u, t) du 1 was
interpreted as the probability that the neuron under consideration has not yet
fired since its last spike at . The value of the integral decreases
therefore over time. On the other hand, if a neuron in the population fires,
it remains part of the population. Apart from a reset of the membrane
potential nothing changes. Thus the integral over the density remains constant
over time, i.e.,
p(u, t) du = 1 .
The normalization to unity expresses the fact that all neurons have a membrane
potential below or equal to threshold.
Third, the fraction of neurons that `flow' across threshold per unit of time
is the (expected value of) the population activity A(t). If we denote the
flux across threshold as
J(, t), we have
A(t) = J(, t) .
Due to the reset, the neurons that `disappear' across threshold, `reenter' at
the reset potential ur. Hence, the membrane potential density at u = ur
increases at a rate proportional to A(t). More specifically, we have a
A(t) (u - ur) at the reset potential ur that
balances the loss that is due
to the movement across the threshold. The
value of A(t) is given by Eq. (6.13).
We assume that all neurons in the population receive the same driving current
Iext. In addition each neuron receives stochastic background
input. We allow for different types of synapse. An input spike at a synapse
of type k causes a jump of the membrane potential by an amount wk. The
effective spike arrival rate (summed over all synapses of type k) is denoted
as . While the mean spike arrival rates (t) are identical for
all neurons, we assume that the actual input spike trains at different neurons
and different synapses are independent6.1. With these
assumptions, the dynamics for
is in analogy to Eq. (5.88)
p(u, t) - - u + RIext(t) p(u, t)
+ (t) p(u - wk, t) - p(u, t) + A(t) (u - ur) .
The first two terms on the right-hand side describe the
continuous drift, the third term the jumps caused by stochastic spike arrival,
and the last term describes the reset.
Because of the firing condition, we have
p(u, t) = 0 for
u > .
In order to calculate the population activity A(t), we need to determine the
flux across threshold. To keep the argument slightly more general, we will
consider the flux J(u0, t) across an arbitrary reference potential u0,
J(u0, t) = Jdrift(u0, t) + Jjump(u0, t) ,
Jdrift accounts for the continuous drift of the membrane
potential during the time when no input spike arrives.
Jjump is due
to excitatory and inhibitory spike arrival.
Figure 6.2:A. All trajectory that are less than wk
below u0 cross u0 upon spike arrival.
B. The drift
Jdrift(u0, t) depends on
density of trajectories and on the slope
with which the trajectories cross the boundary
Jjump, let us consider excitatory input wk > 0 first.
All neurons that have a membrane potential ui with
u0 - wk < uiu0
will jump across the reference potential u0 upon spike arrival at synapse
k; cf. Fig. 6.2A. Since the rate of spike arrival at
synapse k is , the total flux caused by input spikes at all
Jjump(u0, t) = p(u, t) du .
Jdrift(u0, t) through the reference potential u0 is
given by the density p(u0, t) at the potential u0 times the momentary
du/dt; cf. Fig. 6.2B. With
du/dt = [- u + RIext(t)]/ we have
Jdrift(u0, t) = [- u0 + RIext(t)] p(u0, t) .
The total flux at the threshold
u0 = yields the population activity
A(t) = [- + RIext(t)] p(, t) + p(u, t) du .
Since the probability density vanishes for
u > , the sum over the
can be restricted to all excitatory synapses.
Eqs. (6.14) and (6.18) define the
dynamics in a population of integrate-and-fire neurons with stochastic
188.8.131.52 Continuity equation
In this subsection we aim at an interpretation of
Eqs. (6.14) - (6.18). Let us consider the
portion of neurons with a membrane potential between u0 and u1;
cf. Fig. 6.3. The fraction of neurons with u0 < u < u1
increases if neurons enter from below through the boundary u0 or from above
through the boundary u1. A positive flux J(u, t) > 0 is defined as a flux
towards increasing values of u. Since trajectories cannot simply end, we
have the conservation law
p(u', t) du' = J(u0, t) - J(u1, t)
Taking the derivative with respect to the upper boundary u1 and changing
the name of the variable from u1 to u yields the continuity equation,
p(u, t) = - J(u, t) foruurandu ,
which expresses the conservation of the number of trajectories. At u = ur
u = special care has to be taken because of the reset.
u > the flux vanishes
because neurons that pass the threshold are reset.
neurons that have fired start a new trajectory at ur, we have a `source of
new trajectories' at u = ur, i.e., new trajectories appear in the interval
[ur - , ur + ] that have not entered the interval through one
of the borders. Adding a term
A(t) (u - ur) on the right-hand side of
(6.20) accounts for this source of trajectories.
If we insert the explicit
form of the flux that we derived in Eqs. (6.16)
and (6.17) into the continuity equation (6.20),
we arrive once again at Eq. (6.14). For a numerical
implementation of Eq. (6.20) we refer the reader to the
literature (Nykamp and Tranchina, 2000; Omurtag et al., 2000).
The number of trajectories in the interval [u0, u1]
changes if one of the trajectories crosses the
boundary u0 or u1.
For a large number of neurons this fact is described
by the continuity equation; cf. Eq. (6.20).
Schematic figure where only three trajectories are shown.
184.108.40.206 Diffusion approximation
In the limit of small jump amplitudes wk, the density dynamics
(6.14) can be approximated by a diffusion equation. To show
this we expand the right-hand side of Eq. (6.14) into a
Taylor series up to second order in wk. The result is the Fokker-Planck
- - u + RIext(t) + (t) wkp(u, t)
+ (t) wk2p(u, t)
+ A(t) (u - ur) + () .
The term with the second derivative describes a `diffusion' in terms of the
membrane potential. The firing threshold acts as an absorbing boundary so
that the density at threshold vanishes,
p(, t) = 0 .
In order to calculate the flux through the threshold we expand
Eq. (6.18) in wk about
u = and obtain
A(t) = - ,
where we have defined
(t) = (t) wk2 .
Eqs. (6.21) - (6.23) together
with the normalization (6.12) define the dynamics of a
homogeneous population of integrate-and-fire units with `diffusive' noise.
For a more detailed discussion of the diffusion limit
see Chapter 5.5,
in particular Eqs. (5.89) and (5.91).
We assume that the total input
h0 = RIext + wk is constant. In the stationary state, the temporal derivative on
the left-hand-side of Eq. (6.21) vanishes. The terms on
the right-hand side can be transformed so that the stationary Fokker-Planck
0 = - J(u) + A0 (u - ur) ,
where A0 is the population activity (or mean firing rate) in the stationary
J(u) = p(u) - p(u)
is the total flux; cf. Eq. (6.20). The meaning of
Eq. (6.25) is that the flux is constant except at u = ur
where it jumps by an amount A0.
Similarly, the boundary condition
p(, t) = 0
implies a second discontinuity of the flux at
u = .
With the results from Chapter 5.5 in mind,
we expect that
the stationary solution approaches a Gaussian
u - . In fact, we can check
easily that for any constant c1
with some constant c2. We have written the above expression as a product
of two terms. The first factor on the right-hand side is a standard
Gaussian while the second factor guarantees that
p(u) 0 for
u. If we insert Eq. (6.28) in
(6.25) we can check that it is indeed a solution. The constant
c2 is proportional to the flux,
c2 = 2 J(u) forur < u .
The solution defined by Eqs. (6.27) and
(6.28) must be continuous at u = ur. Hence
c1 = exp dx .
Finally, the constant c2 is determined by the normalization condition
(6.12). We use Eqs. (6.27),
(6.28), and (6.30) in
(6.12) and find
= f (x, u) dx du + f (x, u) dx du = f (x, u) du dx ,
f (x, u) = exp - exp .
Figure 6.4B shows the stationary density p(u) for
different amplitudes of the noise.
The activity A0 is identical to the flux J(u) between ur and
and therefore proportional to the constant c2; cf. Eq. (6.29). If we express the integral over u in Eq. (6.31) in
terms of the error function,
erf(x) = exp(- u2) du, we obtain
A0-1 = expx2 1 + erf(x) dx ,
which is identical to expression (5.104) obtained in
Figure 6.4:A. Membrane potential trajectories of 5 neurons (R = 1 and
= 10 ms) driven by a constant background current I0 = 0.8 and
stochastic background input with
= = 0.8 kHz and
w± = ±0.05. These parameters correspond to h0 = 0.8 and
= 0.2 in the
diffusive noise model. B. Stationary membrane potential
distribution in the diffusion limit for
= 0.2 (solid line),
= 0.1 (short-dashed line), and
= 0.5 (long-dashed line).
= 1). C. Mean activity of a population of
integrate-and-fire neurons with diffusive noise as a function of h0 for
four different noise levels, viz. (from top to bottom)
= 1.0, = 0.5, = 0.2 (solid line),
= 0.1, = 0.0.
6.2.2 Spike Response Model Neurons with Escape Noise
In this section we develop a density formalism for spike response neurons,
similar to the membrane potential density approach for integrate-and-fire
neurons that we have discussed in the preceding section. The main difference
is that we replace the membrane potential density p(u, t) by a refractory
density q(r, t), to be introduced below.
We study a homogeneous population of SRM neurons with escape noise. The
membrane potential of the neurons,
u(t) = (t - ) + hPSP(t|) ,
depends on their refractory state,
r = t - 0 ,
i.e., on the time that has passed since the last spike. If we know r and the
total input current in the past, we can calculate the membrane potential,
u(t) = (r) + hPSP(t| t - r) .
Given the importance of the refractory variable r, we may wonder how many
neurons in the population have a refractory state between r0 and
r0 + r.
For a large population (
fraction of neurons with a momentary value of r in the interval
[r0, r0 + r] is given by
= q(r, t) dr ,
where q(r, t) is the refractory density.
The aim of this section is to describe the dynamics
of a population of SRM neurons by
the evolution of q(r, t).
We start from the continuity equation,
q(r, t) = - Jrefr(r, t) ,
where we have introduce the flux
Jrefr along the axis of the
refractory variable r.
As long as the neuron does not fire, the variable
r = t - increases at a speed of
dr/dt = 1. The flux is the
density q times the velocity, hence
Jrefr(r, t) = q(r, t) = q(r, t) .
The continuity equation (6.38) expresses the fact that, as
long as a neuron does not fire, its trajectories r(t) can neither start nor
end. On the other hand, if a neuron fires, the trajectory stops at the
current value of r and `reappears' at r = 0. In the escape rate formalism,
the instantaneous firing rate of a neuron with refractory variable r is
given by the hazard
(t| t - r) = f[(r) + hPSP(t| t - r)] .
If we multiply the hazard (6.40) with the density q(r, t), we get the
loss per unit of time,
Jloss = - (t| t - r) q(r, t) .
The total number of trajectories that disappear at time t due to firing
is equal to the population activity, i.e.,
A(t) = (t| t - r) q(r, t)dr .
The loss (6.41) has to be added as a `sink' term on the right-hand
side of the continuity equation, while the activity A(t) acts as a source at
r = 0. The full dynamics is
This partial differential equation is the analog of the Fokker-Planck equation
(6.21) for the membrane potential density of
integrate-and-fire neurons. The relation between the two equations will be
discussed in Section 6.2.3.
Equation (6.43) can be rewritten in form of an integral equation
for the population activity. The mathematical details of the integration will
be presented below. The final result is
A(t) = PI(t|) A() d ,
PI(t|) = (t|) exp - (t'|) dt'
is the interval distribution of neurons with escape noise;
cf. Eq. (5.9). Thus, neurons that have fired their last spike at
time contribute with weight
PI(t|) to the activity at time
t. Integral equations of the form (6.44) are the starting
point for a formal theory of population activity;
cf. Section 6.3. For a numerical implementation of population
dynamics, it is more convenient to work directly on the level of the density
equations (6.43). A simple discretization scheme for numerical
implementations is discussed below.
All neurons that have fired together at time form a group that moves
along the r-axis at constant speed. To solve Eq. (6.43) we turn
to a frame of reference that moves along with the group.
We replace the
variable r by
t - r and define a new density
Q(, t) = q(t - , t) ,
t. The total derivative of Q with respect to t is
Q(, t) = q(r, t) + q(r, t)
dr/dt = 1. We insert Eq. (6.43) on the right-hand side
of (6.47) and obtain
Q(, t) = - (t|) Q(, t) .
The partial differential equation (6.43) has thus been transformed
into an ordinary differential equation that is solved by
Q(, t) = Q(, t0) exp - (t'|) dt' ,
Q(, t0) is the initial condition, which is still to be fixed.
From the definition of the refractory density q(r, t) we conclude
that q(0, t) is the proportion of neurons at time t that have just fired,
q(0, t) = A(t) or, in terms of the new refractory density,
Q(t, t) = A(t). We can thus fix the initial condition in Eq. (6.49) at
t0 = and find
For a numerical implementation of Eq. (6.43)
we discretize the refractory variable r in steps of t
nk(t) = q(kt + s, t) ds.
The normalization is
nk(t) = 1 for all t.
The probability that a neuron with
r = kt fires in a time step t is
If refractory effects are negligible for
then we can truncate the summation at
All neurons with
kkmax have the same firing probability
If we introduce the normalization
nk = 1 -
we arrive at the update rule
n0(t + t) = Pfree + Pfree - PF(k) nk(t)
nk(t + t) = 1 - PF(k - 1) nk-1 .
Note that n0 is the number of neurons that
fire in the interval
[t, t + t]. Hence we set
A(t + t) = n0(t + t)/t.
The above algorithm allows for a rapid integration of the
population density equation (6.43).
We simulate a population of integrate-and-fire neurons with escape rate given
by Eq. (5.116). At t = 100ms a time-dependent input is switched on.
The population activity A(t) calculated numerically by iterating
Eqs. (6.56) and (6.57) responds to
the input in a non-linear fashion; cf. Fig. 6.5. The
population activity of the model with escape noise can be compared to that of
a model with diffusive noise. The results are strikingly similar. We note
that even the rapid components of the input current are, at least partially,
visible in the population activity.
6.2.3 Relation between the Approaches
Density methods have proven to be a useful tool
for the simulation and analysis of the behavior
of large populations of neurons.
The two approaches that we have discussed in this section,
viz., membrane potential densities and refractory densities,
are closely related.
For noise-free neurons driven by a constant super-threshold
stimulus, the two mathematical formulations
are, in fact, equivalent and related
by a simple change of variables.
But even for noisy neurons with sub-threshold stimulation
the two approaches are comparable;
cf. Fig. 6.5.
Both methods are linear in the densities
and amenable to efficient numerical implementations.
The formal mathematical relation between the two
approaches is shown at the end of this section.
Figure 6.5:A. Activity of a population of integrate-and-fire neurons
with diffusive noise (solid line) or escape noise
(dashed line). At t = 150ms a time-dependent current
is switched on.
B. Time course of the input current.
The current contains 5 randomly chosen frequencies components
between 5 and 500Hz.
= 1, R = 1,
= 10 ms
and ur = 0; diffusive noise model with
escape rate given by Eq. (5.116)
the diffusive noise model has been simulated
using membrane potential densities while
the escape noise model has been simulated using
What are the main differences between the two approaches?
Both are single-variable density methods.
In the refractory density method
the relevant variable is the refractory variable
r. The refractory method is therefore
not compatible with diffusive noise
since diffusive noise generates a
distribution of membrane potentials
so that we would need a description by two-dimensional densities.
Diffusive noise can, to a certain degree,
be replaced by a suitably chosen
escape function f in a noisy threshold model.
But while escape rate models can mimic
the consequence of stochastic spike arrival
for the population activity,
escape rates cannot serve as a model
for membrane potential
fluctuations that are commonly seen
in intracellular recordings.
On the other hand, refractory density methods combined with escape noise are
well suited to describe neuronal refractoriness. It is, for example,
straightforward to implement a spiking neuron model with two or more time
scales, a slow one for the spike after-potential and a fast one for the
postsynaptic potential -- as it is typical for Hodgkin-Huxley type
neurons; cf. Chapters 2.2 and 4.3. A
finite rise time of the postsynaptic potential can be included in the analysis
without any additional effort. While it is difficult to include reversal
potentials and adaptation in the analysis of refractory density models, both
effects can be incorporated phenomenologically in numerical implementations by
a few additional macroscopic variables. An advantage of refractory densities
combined with escape noise is that the density equations can be integrated
formally. While the integral representation is not useful for numerical
solutions, it is useful for an interpretation of the dynamics in terms of
interspike interval distributions.
Membrane potential density methods
work best for integrate-and-fire neurons.
Since the membrane potential of an integrate-and-fire
neuron is described by a single differential
the derivation of the density equation
A numerical integration of the density equations
causes no major problems if the threshold and reset
mechanism as well as the `drift' and `jump' terms
are carefully implemented (Nykamp and Tranchina, 2000).
Reversal potential can be included
at no extra cost -- both in simulations
and in the analysis (Abbott and van Vreeswijk, 1993; Nykamp and Tranchina, 2000; Omurtag et al., 2000).
arrival can be described explicitly
and both the amplitude and
the frequency spectrum of the
fluctuations of the membrane potential
can be predicted.
While the formal approach
can be extended to detailed neuron models,
the implementation of high-dimensional
density equations requires efficient
numerical implementations (Knight, 2000).
Most implementations, so far, have been
restricted to integrate-and-fire
neurons with or without reversal potentials.
220.127.116.11 From membrane potential densities to
phase densities (*)
For constant input current I0
it is possible to transform
membrane potential densities
into phase densities.
A description by phase variables is interesting in itself;
at the same time this section will introduce the methods
that we will use below
for the transformation
from membrane potential densities to refractory densities.
We consider a population of noise-free integrate-and-fire neurons
stimulated by a constant super-threshold current.
To keep the arguments general, we
consider a non-linear integrate-and-fire neuron
(cf. Chapter 4.1)
u = F(u) = - + I0 .
The voltage dependence of and r can account for voltage dependent
conductances and reversal potentials. For
(u) = and R(u) = R
we have the standard equation of an integrate-and-fire neuron.
For I0 sufficiently large, the neuron will reach the firing threshold
= 1. After firing, the membrane potential is reset to ur = 0.
Integration of the differential equation (6.58) yields the
u(t) = (t - ) ,
d/dt = F[(t - )]. Here denotes the last firing
time of the neuron; cf. Fig. 6.6. For constant
super-threshold input, the process of integration and reset repeats with a
period T, i.e.,
u(t + T) = u(t). We may therefore introduce a phase
variable that increases from zero to T and is after each firing reset
to zero, i.e.,
= tmodT. This allows us to rewrite the
membrane potential trajectory as
u(t) = ().
Figure 6.6:A. With constant input, a noise-free (non-linear) integrate-and-fire
neuron fires regularly with period T.
We define a phase variable by
= t - and write the trajectory as
u = ().
B. The membrane potential density p(u, t)
is related to the phase density q(, t)
d/d. The number of neurons
with refractory variable between and
(shaded region) is equal to the number of neurons
membrane potential between
u0 = () and
u1 = ().
We can introduce a phase densityq(, t) so that
q(, t) d is the fraction of neurons with
phase variable in the interval
[,]. The phase density q is
related to the membrane potential density p by
q(, t) d = p(u, t) du (see Fig. 6.6B) so that
q(, t) = p(u, t) = p[(), t] F[()]
where we have used
d/dt = 1.
The normalization is therefore
q(, t) d = p((), t) F[(r)] d = p(u, t) du = 1
We now want to derive the continuity equation for the phase variable
q(, t). We start from the continuity equation (6.20)
of the membrane potential densities which reduces in the absence of noise to
p(u, t) = - p(u, t) F(u) , for 0 < u < 1 .
The term in square brackets is the drift current with
F(u) = du/dt.
We use the product rule to evaluate the derivative on the right-hand side and
multiply by F(u). The result is
The left-hand side is identical to
cf. Eq. (6.60). Taking the partial derivative of
Eq. (6.60) with respect to yields the right-hand side of
(6.63). Thus Eq. (6.63) can be rewritten as
q(, t) = - q(, t) for 0 < < T ;
cf. (Abbott and van Vreeswijk, 1993). The phase variable plays a role similar to the
refractory variable r. In fact, Eq. (6.64) is identical to
the noise-free refractory density equation (6.43). We emphasize,
however, that phase variables are restricted to periodic behavior and require
therefore constant super-threshold input - in contrast to refractory
In this paragraph we want to show the formal relation between the dynamics of
p(u, t) and the evolution of the refractory densities q(r, t). We focus on
a population of standard integrate-and-fire neurons with escape noise. The
equation of the membrane potential
u(t|) = exp - + exp - I(t') dt'
can be used to define a transformation from voltage to refractory variables:
ur = t - . It turns out that the final equations are
even simpler if we take instead of r as our new variable. We
therefore consider the transformation
Before we start, we calculate the derivatives of Eq. (6.65).
The derivative with respect to t yields
u/t = [- u + RI(t)]/ as expected for integrate-and-fire neurons. The
derivative with respect to is
The densities in the variable are denoted as
Q(, t). From
Q(, t) d = p(u, t) du we have
Q(, t) = p[u(t|), t] F(t,) .
We now want to show that the differential equation for the density
Q(, t) that we derived in (6.48),
Q(, t) = - (t|) Q(, t) , for < t ,
is equivalent to the partial differential equation for the membrane potential
densities. If we insert Eq. (6.67) into Eq. (6.68) we
F + F + p = - pF .
For the linear integrate-and-fire neuron we have
F/t = - F/. Furthermore for
RI(t) > we have F 0. Thus we
can divide (6.69) by F and rewrite Eq. (6.69) in the form
= - p(u, t) - f (u - ) p(u, t) , forur < u < ,
where we have used the definition of the hazard via the escape function
(t|) = f[u(t|) - ] and the definition of the reset
ur = . If we compare Eq. (6.70) with the
Fokker-Planck equation (6.21), we see that the main
difference is the treatment of the noise. For noise-free integrate-and-fire
(t|) = 0 for
u) the equation
(6.21) for the membrane potential densities is
therefore equivalent to the density equation
Q(, t)/t = 0; cf. Eq. (6.68).