ARTICLE

Communicated by Nicolas Brunel

Mean-Field Approximations for Coupled Populationsof Generalized Linear Model Spiking Neurons
with Markov RefractorinessTaro Toyoizumi
taro@neurotheory.columbia.edu
Department of Neuroscience, Columbia University, New York, NY 10032, U.S.A.

Kamiar Rahnama Rad
kamiar@stat.columbia.edu
Department of Statistics, Columbia University, New York, NY 10027, U.S.A.

Liam Paninski
liam@stat.columbia.edu
Department of Neuroscience and Department of Statistics, Columbia University,
New York, NY 10027, U.S.A.

There has recently been a great deal of interest in inferring network connectivity from the spike trains in populations of neurons. One class ofuseful models that can be fit easily to spiking data is based on generalized linear point process models from statistics. Once the parameters
for these models are fit, the analyst is left with a nonlinear spiking network model with delays, which in general may be very difficult to understand analytically. Here we develop mean-field methods for approximating the stimulus-driven firing rates (in both the time-varying andsteady-state cases), auto- and cross-correlations, and stimulus-dependent
filtering properties of these networks. These approximations are valid
when the contributions of individual network coupling terms are smalland, hence, the total input to a neuron is approximately gaussian. These
approximations lead to deterministic ordinary differential equations thatare much easier to solve and analyze than direct Monte Carlo simulationof the network activity. These approximations also provide an analytical
way to evaluate the linear input-output filter of neurons and how thefilters are modulated by network interactions and some stimulus feature.
Finally, in the case of strong refractory effects, the mean-field approximations in the generalized linear model become inaccurate; therefore,
we introduce a model that captures strong refractoriness, retains all of
the easy fitting properties of the standard generalized linear model, andleads to much more accurate approximations of mean firing rates andcross-correlations that retain fine temporal behaviors.

Neural Computation 21, 1203–1243 (2009)
C 2008 Massachusetts Institute of Technology

1204

T. Toyoizumi, K. Rad, and L. Paninski

1 Introduction
One of the fundamental problems in the statistical analysis of neural data is
to infer the connectivity and stimulus dependence of a network of neurons
given an observed sequence of extracellular spike times in the network(Brown, Kass, & Mitra, 2004). Estimating the parameters of network models
from data is in general a computationally difficult problem because of thehigh dimensionality of the parameter space and the possible complexity ofthe model’s objective function given an observed data set of network spiketimes.
One class of models that has proven quite useful is known in the statistics
literature as the generalized linear model (GLM; Brillinger, 1988; Chornoboy,
Schramm, & Karr, 1988; see also Martignon et al., 2000; Kulkarni & Paninski,
2007; Nykamp, 2007), variants of which have been applied successfullyin the hippocampus (Harris, Csicsvari, Hirase, Dragoi, & Buzsaki, 2003;
Okatan, Wilson, & Brown, 2005), motor cortex (Paninski, Fellows, Shoham,
Hatsopoulos, & Donoghue, 2004; Truccolo, Eden, Fellows, Donoghue, &Brown, 2005), retina (Pillow et al., 2008), and cultured cortical slice (Rigat,
de Gunst, & van Pelt, 2006). In its simplest form, this model incorporates
both stimulus-dependence terms and direct coupling terms between eachobserved neuron in the network; fitting the model parameters leads to an
inhomogeneous (stimulus-driven), nonlinear coupled spiking model withdelay terms. Statistically speaking, the model is attractive due to its explicitly probabilistic nature and because fitting the model parameters is surprisingly simple: under certain simple conditions, the log-likelihood function
with respect to the model parameters is concave, and the maximum likelihood estimation of the parameters can be easily performed using standard
ascent methods (Paninski, 2004; Paninski, Pillow, & Lewi, 2007).
Once the model parameters are obtained, we are left with an obvious
question: What do we do next? One of the key applications of such a networkmodel is to better understand the input-output properties of the network.
For example, we would like to be able to predict the mean firing response of
the network given a novel input and to carve out the impact of the network
coupling terms on this stimulus-response relationship (e.g., what is local
inhibition’s impact on the stimulus filtering properties of the network?). We
would also like to know how the correlation properties of spike trains in
the network might depend on the stimulus, and in general how correlations
might encode stimulus information (Riehle, Grün, Diesmann, & Aertsen,
1997; Hatsopoulos, Ojakangas, Paninski, & Donoghue, 1998; Oram, Hatsopoulos, Richmond, & Donoghue, 2001; Nirenberg, Carcieri, Jacobs, &Latham, 2002; Schnitzer & Meister, 2003; Kohn & Smith, 2005; Pillow et al.,2008).
We can in general study these questions with direct Monte Carlo simulations. However, simulation of a large-scale probabilistic spiking network
is computationally expensive, since we often need to draw many samples

Mean-Field Approximations for GLM Neurons

1205(i.e., run many simulated trials) in order to compute the quantities of interest to the desired precision. More importantly direct numerical simulation
often provides limited analytical insight into the mechanisms underlying
the observed phenomena.
The goal of this study is to investigate how much of the behaviorsof these GLM networks can be understood using standard analyticalmean-field approximations (Renart, Brunel, & Wang, 2003; Wilson &Cowan, 1972; Amit & Tsodyks, 1991; Ginzburg & Sompolinsky, 1994; Hertz,
Krogh, & Palmer, 1991; Kappen & Spanjers, 2000; Gerstner & Kistler, 2002;
Meyer & van Vreeswijk, 2002). In particular, we develop approximations
for the mean firing rates of the network given novel stimuli, as well as
the auto- and cross-correlation and input-output filtering properties ofthese networks. These approximations are valid when the contributionsof individual network coupling terms are small and lead to deterministic
ordinary differential equations that are much easier to solve and analyze
than direct Monte Carlo simulation of the network activity.
However, in the case of strong refractory effects, these mean-field approximations become inaccurate, since the spike history terms in the generalized linear model must be large to induce strong refractoriness, and thispushes our approximations beyond their region of accuracy. Therefore we
introduce a modified model: a generalized linear model with Markovian
refractoriness. This model has several advantages in this setting: it capturesstrong refractoriness, retains all of the easy-fitting properties of the standard
generalized linear model, and leads to much more accurate mean-field approximations.
2 The Generalized Linear Point Process Model
We begin by precisely defining the generalized linear point process model(see Figure 1). We consider N recurrently connected spiking neurons, modeled as a multivariate point process (Snyder & Miller, 1991). The firing rate
(conditional intensity function) of neuron i at time t is a function of the total
input, ui (t), the neuron receives:
λi (t) = f (ui (t)),(2.1)

where f (.) is a smooth, nonnegative, and monotonically increasing function.
For numerical simulations, we use exponential nonlinarity, f (u) = e u , but
the analysis in this article is not limited to this exponential nonlinearity. The
total input to the neuron is expressed as the sum of external input Ii andrecurrent input Hi :ui (t) = Ii (t) + Hi (t).

(2.2)

1206

T. Toyoizumi, K. Rad, and L. Paninski

Figure 1: A schematic illustration of the generalized linear point process model.

The recurrent input is modeled as a sum of identical responses caused byeach presynaptic spike and written as
Hi (t) =N jwi (t − t j,n ),(2.3)

j=1 t j,n <tjwhere wi (t) is a coupling term from neuron j (upper index) to i (lower
index), and t j,n is the nth spike time of neuron j. Without loss of generality,j
we may express each coupling term wi (t) as a weighted sum of exponentialfunctions of different time constants. As we will see below, this sum-ofexponentials representation is extremely useful: we will take full advantageof the Markovian nature of the decaying exponential function. To keep
notation under control, we will restrict our attention to the case of a singlejexponential function with time constant τi ,jjj

wi (t) = J i exp(−t/τi )(t),(2.4)

where  is the Heaviside step function, which takes one for positive arguj
ments and zero otherwise; generalizations to the case where wi (t) are given
by a sum of exponentials with different time constants will be left implicit
and will be straightforward in all cases. If we write the output spike trainof neuron i as
Si (t) =δ(t − ti,n ),(2.5)
nMean-Field Approximations for GLM Neurons

1207

the recurrent input of equation 2.3 is rewritten as the convolution of thesynaptic filter and the presynaptic spike train:

Hi (t) =N jds wi (t − s)S j (s)

j=1

=

N


j

j

J i h i (t)

(2.6)

j=1

with
j

h i (t) =

 t
−∞

j

e −(t−s)/τi S j (s) ds.

(2.7)

This model is a generalization of the inhomogeneous Poisson process, since
past output spikes modulate the firing intensity function. If we set all thejcoupling weights J i to zero, we recover the the inhomogeneous Poissonprocess model. The model is also closely related to the spike-response modelof Gerstner and Kistler (2002) and the inhomogeneous Markov interval
process of Kass and Ventura (2001), in which the firing rate depends on
only the last observed spike. Note that in the GLM, the firing intensityof a neuron may depend not only on the last spike of the neuron but on
the past several spikes of this neuron (see, e.g., Paninski, 2004; Paninski
et al., 2007, for further discussion). As mentioned in section 1, this modelis attractive because it is fairly easy to fit to data and seems to capture thefiring statistics of neural populations in a variety of experimental settings.
However, in this article, we will not discuss the estimation problem but
rather limit our attention to the problem of determining the network firing
rate, auto- and cross-correlation functions, and input-output properties,
given fixed, known model parameters.
3 Master Equation
The main object we would like to compute here is P(h, t): this is the disj
Ttribution of the state variable h = (h 11 , . . . , h NN ) at time t. These h i are thesolutions of the multivariate coupled Markov process,jj

dh i (t)h (t)= − i j + S j (t),
dt
τi

(3.1)

and therefore solving for P(h, t) will allow us to compute a number ofimportant quantities. For example, to compute the mean firing intensity

1208

T. Toyoizumi, K. Rad, and L. Paninski

νi (t) = E[Si (t)] of neuron i given some input {Ii (t)}, we may computeνi (t) = P(h, t)λi (t)d h=P(h, t) f  Ii (t) +J i h i  d h,jjjwhere E[.] is the average over the spike history {Si (t)|i = 1, . . . , N, 0 ≤ t <T} in an entire duration T given external input {Ii (t)}. Hence, we can readthe mean firing rate off directly from P(h, t).
We begin by writing down the master equation—the Kolmogorov forward equation that governs the time evolution of this Markov process(Oksendal, 2002): j
∂ P(h, t)  1 ∂ h i P(h, t)=
λ j (h − e j , t)P(h − e j , t)
+
j
j
∂t
τi
∂h i
i, j
j

λ j (h, t)P(h, t),
(3.2)
−
j
j

j

with λi (h, t) = f (Ii (t) + j J i h i ) as in equation 2.1 and the term e j denoting the vector with its {i  j  } slot δ j j  . Note that there are no diffusion
terms here; the random nature of the stochastic differential equation entersinstead through the jump terms.
While this jump-advection partial differential equation (PDE) in principle completely describes the behavior of this system, unfortunately thisequation remains difficult to handle analytically (due to the nonlinearity inthe λ terms) and numerically (due to the large dimensionality of the statevariable h). We will pursue a direct analysis of this PDE elsewhere; here, wewill attack this system using a different approach, based on a self-consistent
expansion of the first and second moments of Si (t).
In the one-dimensional case (i.e., a single neuron, with one h(.)), however, we may easily solve the PDE numerically to obtain an exact solution
for the mean firing rate. Figure 2 compares this exact solution versus themean-field approximations derived in the following section. The meanfield approximation is reasonably good here but not perfect, because of thefluctuation in the recurrent input (J = −1, τ = 10 ms).
4 Mean-Field Approximation of GLM Statistics
Now we turn to the main topic of this article. As emphasized above,
we would like to calculate the mean-firing intensities of neurons, νi (t) =E[Si (t)], or the spike-cross-correlation function between neurons, φi j (t, t  ) =E[Si (t)S j (t  )]. The calculation of these statistics is difficult, however, because

Mean-Field Approximations for GLM Neurons

1209

Figure 2: Direct solution of the master equation in the one-dimensional case.
(A) Top: A single neuron is simulated with a pulse input I (t) shown. Bottom:
Comparison of the firing rates computed via direct Monte Carlo sampling (simulation); the inhomogeneous Poisson case, with no history dependence, that is,
w = 0 (noH); the master equation of equation 3.2 (exact); and the mean-field approximation methods in the limit of Nc → ∞ introduced in section 4. The meanfield approximation gives good approximations of the output firing intensity. (B)Top: The probability distribution P(h, t); grayscale axis indicates the height ofthe probability density. The self-history term h(t) was the convolution of the output spike train and a single exponential filter with time constant τ = 10 ms andweight J = −1; thus, the neuron was inhibited for a brief period following each
spike. Note the discontinuous behavior visible at stimulus onset, where probability mass jumps significantly toward h = 1, then decays back down toward
zero, according to equation 3.1. Bottom: The mean of h calculated from the probability distribution P(h, t). This E[h] also increases following the input pulse,which has the effect of slowing the firing rate. Since the weight J is negative, thehistory term is inhibitory here. (A color version of this figure is available online at
http://www.mitpressjournals.org/doi/suppl/10.1162/neco.2008.04-08-757.)of the nonlinearity f and the recurrent input term H in the model. A typicalway of approximating these recurrent effects is to replace the interaction
term H by its mean field (Renart et al., 2003; Wilson & Cowan, 1972; Amit
& Tsodyks, 1991; Ginzburg & Sompolinsky, 1994; Hertz et al., 1991; Kappen& Spanjers, 2000; Gerstner & Kistler, 2002; Meyer & van Vreeswijk, 2002).
Let us assume that each neuron receives input from some subset, Nc , ofthe total population of neurons N. If the synaptic strength scales J ∼ 1/Nc
with Nc , and assuming that the population is in the “asynchronous” state,then it is known that the cross-covariance function between two neurons φi j (t, t  ) = Cov[Si (t), S j (t  )] scales with 1/Nc for large Nc (Ginzburg& Sompolinsky, 1994). This means that the contribution of each network

1210

T. Toyoizumi, K. Rad, and L. Paninski

interaction term is negligible and all the neuron activities are independent
in the Nc → ∞ limit. Hence, calculating the mean-firing intensities andcross-correlations is easy in this limit.
We first derive the mean-firing intensity in the limit of Nc → ∞ and thenevaluate the finite-size effect of order 1/Nc . Note that this mean-field solution is a good approximation not only when Nc is very large; equivalently,
this is a good approximation when the synaptic coupling J is small for fixed
Nc , where the contribution of an individual synapse is small.
Now our first key approximation is to appeal to the central limit theorem
applied to the sum of equation 2.6, and to assume that the recurrent input
H is roughly gaussian. We can express the mean firing intensities of theneurons in our network as a function of two sufficient state variables: themean recurrent input E[H] and the variance of the recurrent input Var[H],νi (t) = E[λi (t)]= E[ f (Ii (t) + Hi (t))]≈ F (Ii (t) + µi (t), σi (t)),(4.1)

where, in the third line, the distribution of the recurrent input, Hi (t), is
assumed to be gaussian with mean µi (t) = E[Hi (t)] and variance σi2 (t) =Var[Hi (t)], and the expectation was replaced by a gaussian integral:f (u) − (u−µ)2 2
F (µ, σ ) =e 2σ du.
(4.2)
√
2πσ 2
Note for an exponential nonlinearity f (u) = e u , the gaussian integral isgiven by F (µ, σ ) = exp(µ + σ 2 /2).
Next we calculate the cross-correlation function,
φi j (t, t  ) = E[Si (t)S j (t  )]= δi j δ(t − t  )E[Si (t)] + (t − t  )E[λi (t)S j (t  )]+ (t  − t)E[Si (t)λ j (t  )]= δi j δ(t − t  )νi (t) + (t − t  )νi j (t | t  )ν j (t  )+ (t  − t)ν ji (t  | t)νi (t),(4.3)

where in the second line, we decomposed the correlation function into
the three terms: the first term shows the simultaneous correlation with
its amplitude proportional to the firing intensity due to the point-process
nature of the spike train, the second term shows the correlation when t > t (the spike variable Si (t) is averaged given past spike history to yield λi (t)),
and the third term is the correlation for t  > t. In the third line, we introduced
conditional firing intensity,
νi j (t | t  ) = E[Si (t) | spike of j at t  ] =E[Si (t)S j (t  )],E[S j (t  )](4.4)

Mean-Field Approximations for GLM Neurons

1211

which describes the probability of neuron i firing at t given the spike of j att  . We also introduce an abbreviation in the following and write E[.|S j (t  )]instead of E[.|spike of j at t  ]. Similar to the approximation in equation 4.1,
we assume a gaussian distribution of H given a spike of neuron j at time t 
to findνi j (t | t  ) = E[λi (t) | S j (t  )]= E[ f (Ii (t) + Hi (t)) | S j (t  )]≈ F (Ii (t) + µi j (t | t  ), σi j (t | t  ))(4.5)

for t > t  with conditional mean, µi j (t | t  ) = E[Hi (t) | S j (t  )], and conditional variance, σi2j (t | t  ) = Var[Hi (t) | S j (t  )]. Note that equation 4.5 is valid
only for t > t  ; if t < t  instead, we evaluate ν ji (t  | t) and use the Bayes
theorem to calculate νi j (t | t  ) = ν ji (t  | t)νi (t)/ν j (t  ). The unconditional andconditional mean and variance of H in equations 4.1 and 4.5 are evaluatedasµi (t) =J ik αik (t),kµi j (t | t  ) =J ik αikj (t | t  ),k

σi2 (t) =J ik J il Aikl (t),k,lσi2j (t | t  ) =J ik J il Aiklj (t | t  ),

(4.6)

k,l

where we defined
αik (t) = E h ik (t) = t−∞ds e −(t−s)/τi νk (s),αikj (t | t  ) = E h ik (t) | S j (t  ) =

k

 t
−∞

ds e −(t−s)/τi νk j (s | t  ),kAikl (t) = Cov[h ik (t), h li (t)]
=

 t

 t

−∞

ds


−∞



l



l

ds  e −(t−s)/τi −(t−s )/τi φkl (s, s  ),kAiklj (t | t  ) = Cov h ik (t), h li (t) | S j (t  )
=

 t

 t

−∞

ds

−∞

ds  e −(t−s)/τi −(t−s )/τi Cov[Sk (s), Sl (s  ) | S j (t  )],
k

(4.7)

1212

T. Toyoizumi, K. Rad, and L. Paninski

with φkl (s, s  ) = φkl (s, s  ) − νk (s)νl (s  ). The evolution of αik , αikj , and Aikl inequation 4.7 follows simple ordinary differential equations (ODEs) and iseasy to simulate,
α k (t)d kαi (t) = − i k + νk (t),dtτiαi j (t | t  )d k+ νk j (t | t  ),αi j (t | t  ) = −
dt
τik
k

 t+| |
d kl
1
1
k
kl
Ai (t) = −
A
+
(t)
+
lim
e −(t−s)/τi φkl (s, t) ds
i
→0
dt
τik
τil
−∞

 t−| |le −(t−s )/τi φkl (t, s  ) ds +−∞ t+| |Aikl (t)k+lime −(t−s)/τi νkl (s | t)νl (t) ds
kl
→0
τi
−∞

 t−| |le −(t−s )/τi νlk (s  | t)νk (t) ds 
+

=−

−∞

=−

Aikl (t)l+ αilk (t + | t)νl (t) + αik(t − | t)νk (t),
τikl

(4.8)

where 1/τikl = 1/τik + 1/τil , νi j (t | t  ) = νi j (t | t  ) − νi (t), φkl (s, s  ) = νkl(s | s  )νl (s  ) = νlk (s  | s)νk (s), αikj (t | t  ) = αikj (t | t  ) − αik (t), and αilk (t ± |t) = lim →0 αilk (t± | || t). Note that the differentiation of Aikl (t) for k = l isa little tricky because φkl (s, s  ) has a delta peak at s = s  . While the definition of
Aikl (t) in equation 4.7 is clear, we have to make sure in the evaluationd Aikl (t)of dt not to double-count the integration of the delta peak. Hence, inthe third equation of equation 4.8, we included the delta peak along the
integration over s (including the peak at s = t) but excluded this peak fromthe integration along s  (excluding the peak at s  = t).
4.1 Solution in the Mean-Field Limit. In this section we evaluate mean
firing intensities, νi , and cross-correlation functions, φi j , in the mean-field
limit of Nc → ∞ (which is also valid in the noncoupled limit J → 0 for finite
Nc ), that is, synaptic strengths are scaled J ∼ 1/Nc so that the individual
contribution of a spike is negligible for large Nc ; and we also assume, as
discussed above, that the network is in the asynchronous state. Under these
conditions, the cross-covariances between neurons are known to scale as
φi j ∼ 1/Nc for i = j and φi j ∼ 1 for i = j (Ginzburg & Sompolinsky,1994). Then we can easily see σi2 (t) ∼ 1/Nc because Aikl ∼ 1 for k = l and
Aikl ∼ 1/Nc for k = l from equation 4.7. Hence, in the limit of Nc → ∞, weMean-Field Approximations for GLM Neurons

1213

obtain from equations 4.1 and 4.8,
νi (t) → f (Ii (t) + µi (t)) ,J ik αik (t),µi (t) =k

α k (t)d kαi (t) = − i k + νk (t).
dt
τi

(4.9)

This self-consistent update rule corresponds to a well-known mean-field
dynamics of the mean-firing intensity (Ginzburg & Sompolinsky, 1994).
Similarly, the scaling of cross-covariance φi j ∼ 1/Nc for i = j implies
αikj (t | t  ) ∼ 1/Nc for k = j and αikj (t | t  ) ∼ 1 for k = j from equation 4.7.
Moreover the conditional covariance Cov[Sk (s), Sl (s  ) | S j (t  )] scales with
1/Nc for k = l = j and thus σi2j (t | t  ) ∼ 1/Nc . Hence, in the limit of Nc → ∞,we find νi j (t | t  ) → νi (t) for t > t  and, from equation 4.3,
φi j (t, t  ) → δi j δ(t − t  )νi (t).

(4.10)

These equations are very easy to simulate, and the approximation of themean-firing intensity is good even for a small Nc if J is small (see Figure 2A).
On the other hand, the mean-field approximation is not good for large J(see Figures 3 and 4B). This is because, for fixed Nc , the fluctuation of therecurrent input increases with J , making the mean-field assumption (which
neglects these fluctuations) invalid. Figure 4 shows the cross-correlation
functions between two neurons. Because the self-interaction term J ∼ 1/Nc
vanishes for the large Nc limit, the autocorrelation function in this limit cancapture the Poisson peak only at t = t  ; it does not capture any nontrivial
correlation caused by the interaction term J for finite Nc (see Figure 4).
4.2 Estimating the Finite Size Effect. In the previous section, we calculated the mean-firing intensity, νi (t), and the cross-covariance function,
φi j (t, t  ), in the Nc → ∞ limit. Ideally, once provided the true σi2j (t | t  )term, which is of order 1/Nc2 different from σi2 (t), we can iteratively
evaluate equations 4.1 to 4.8, starting from the mean-field limit solutionof equations 4.9 and 4.10 to find the precise approximation of νi (t) andφi j (t, t  ) as long as the gaussian H assumption is valid. However, we need
the third-order correlations to evaluate σi2j (t | t  ), and the equations are notclosed. Here, we give up evaluating O(1/Nc2 ) terms and just evaluate thefinite size corrections of νi (t) and φi j (t, t  ) up to O(1/Nc ) terms. Up to thisorder, it can be argued that σi2j (t | t  ) ≈ σi2 (t), and under this approximation,
we obtain a self-consistent expression for νi (t) and φi j (t, t  ). Equations 4.1

1214

T. Toyoizumi, K. Rad, and L. Paninski

Figure 3: Comparison of the direct Monte Carlo simulation of a single neuron (solid), the firing intensity without recurrent input H (dotted), and themean-field approximation (dashed). (A) Low step input (baseline I = 2, peakI = 4) is applied. The self-inhibition term is set to J = −5 and τ = 2 ms. (B)High step input (baseline I = 2, peak I = 8) is applied. The self-inhibition
term is set to J = −5 and τ = 2 ms. Because of the large nongaussian fluctuations of input caused by the strong self-inhibition, the mean-field approximation loses accuracy here. (A color version of this figure is available online at
http://www.mitpressjournals.org/doi/suppl/10.1162/neco.2008.04-08-757.)

and 4.5 yield, up to the first order of 1/Nc ,νi (t) = F (Ii (t) + µi (t), σi (t)),νi j (t | t  ) = F (Ii (t) + µi j (t | t  ), σi (t)),(4.11)

where σi2j (t | t  ) ≈ σi2 (t) to the first order of 1/Nc . Equation 4.11 can be evaluated using equations 4.6 to 4.8, which we repeat here for convenience:µi (t) =J ik αik (t),kµi j (t | t ) =J ik αikj (t | t  ),k(σi (t)) =2

J ik J il Aikl (t),k,lα k (t)d kαi (t) = − i k + νk (t),dtτiαi j (t | t  )d k+ νk j (t | t  ),αi j (t | t  ) = −
dt
τik
k

Akl (t)d kllAi (t) = − i kl + αilk (t + | t)νl (t) + αik(t − | t)νk (t).
dt
τi(4.12)Mean-Field Approximations for GLM Neurons

1215

Figure 4: Comparison of cross-correlation functions of two neurons calculatedby direct simulation (solid), calculated from the mean-field approximation (seeequations 4.9 and 4.10; dashed), and the mean-field approximation with finite size correction (see equations 4.11 and 4.12; dot-dashed). (A) Weak selfinteraction. Parameters are set to J 11 = −1, J 22 = −1, J 21 = 1, J 12 = −0.5, τ11 =10 ms, τ22 = 10 ms, τ21 = 10 ms, τ12 = 20 ms, and I1 = 2, I2 = 4. The result withfinite size correction is obtained after five iterations of the self-consistent equations. The calculation in the mean-field limit captures only the steady-state firing
intensity, not any temporal cross-correlation function. Cross-correlation functions are well approximated with the finite-size correction terms. (B) Strong selfinteraction: parameters are set to J 11 = −5, J 22 = −5, J 21 = 1, J 12 = −0.5, τ11 =2 ms, τ22 = 2 ms, τ21 = 10 ms, τ12 = 20 ms, and I1 = 2, I2 = 4. Because of thelarge nongaussian fluctuation of the input caused by the strong self-inhibition
term, the mean-field approximation does not work very well. In this case thefinite size solution converges rather slowly. Hence, the result is after 10 iterations of the self-consistent equations. (A color version of this figure is available
online at http://www.mitpressjournals.org/doi/suppl/10.1162/neco.2008.0408-757.)

In principle, the above equations are valid even with time-dependent
input, and several iterations of the above set of equations yield timedependent cross-correlation functions. However, evaluation of the above
equations is computationally expensive with time-dependent input. If one
neuron is connected to Nc surrounding neurons and time is discretized into
T time bins, it takes O(N2 · Nc · T 2 ) operations and memory to evaluate
αikj (t | t  ). If the input is constant in time, αikj (t | t  ) is a function of t − t  only,
and just O(N2 · Nc · T) operations are required to solve the equations. Therefore, we focus on calculating the finite size effect for constant input in this1216

T. Toyoizumi, K. Rad, and L. Paninski

article.1 Generally several iterations of the above set of equations are
required for the convergence of mean firing intensities and crosscorrelations. However, by initially setting those variables to the solution
in the limit of Nc → ∞, we found good approximations after a couple ofiterations in many cases.
The mean firing rate and cross-correlation functions calculated fromequations 4.11 and 4.12 are precise for small J and large Nc because therecurrent input is close to gaussian in this case (see Figure 4A). On the other
hand, for large J and small Nc , the gaussian input approximation is poor,and therefore the approximation becomes worse (see Figure 4B). One key
result here is that both the autocorrelations and the cross-correlations inthese networks depend on the inputs I (t) in general, as can be seen fromequations 4.6 to 4.8. Changing just the baseline (mean) of the input I (t) cansignificantly change the correlation structure in the network as seen in a
number of physiological preparations (e.g., Hatsopoulos et al., 1998; Kohn
& Smith, 2005), even if the coupling terms J and w are unchanged.
4.3 Linear Input-Output Filter. The ODE derived in equation 4.9 isgeneral and may be applied given arbitrary time-dependent input {Ii (t)}.
In this section, we exploit these equations to explore the linear response
properties of these spiking networks. A standard physiological technique
is white noise analysis (Marmarelis & Marmarelis, 1978; Rieke, Warland,de Ruyter van Steveninck, & Bialek, 1997), in which we inject white noise
input in order to study a network’s input-output filtering properties. When
this analysis
is applied to our GLM system, the external input is given byIi (t) = K i (t − s)ξi (s) ds, where K i is neuron i’s linear stimulus filter and(0)ξi (t) is the white noise stimulus with mean ξ (t) = ξi and small variance
(0)
(0)
2


(ξi (t) − ξi )(ξ j (t ) − ξ j ) = σξ δi j δ(t − t ). We keep the input variance small
here in order to treat the input as almost constant and develop a consistentlinear response theory. Note that . describes the average over the fluctuating stimulus ξi (t) here. We use δ X = X − X for any function X of input
stimulus and neglect O(δ I 2 ) terms in the following calculations.
Let us summarize the calculations in the mean-field (Nc → ∞) limit. The
mean-firing intensity of a neuron is given by the following self-consistentequations:νi (t) = f (Ii (t) + µi (t)),
wik (t − s)νk (s) ds.
(4.13)µi (t) =k
1 Note that the computational complexity to calculate the finite size effect is not N2 but
N2 · Nc because each synaptic time constant, τik , is distinct. Hence, we needed to evaluateeach αikj = E[h ik | S j ] term separately. However, if all the synaptic constants onto eachneuron are constant, that is, τik = τi , we can directly evaluate µi j = k J ik αikj without
evaluating αikj separately, and the computational complexity to evaluate the finite size
effect reduces to O(N2 · T).

Mean-Field Approximations for GLM Neurons

1217

Figure 5: Comparison of (A) a simulated firing rate of a single neuron and(B) the self-consistent solutions of equation 4.13 for different input levels,
I (0) , and self-inhibition strengths, J . (C) The percentage errors of the meanfield approximation are shown with respect to the simulated firing intensity. We set τ = 10 ms. (A color version of this figure is available online at
http://www.mitpressjournals.org/doi/suppl/10.1162/neco.2008.04-08-757.)(0)

In particular, for a constant stimulus, ξi (t) = ξi , we obtain
 (0)
(0)
(0) 
νi = f Ii + µi ,

(0)
(0)
µi =
τik J ik νk(4.14)kfrom equation 2.4. Note that X(0) describes responses at constant stimulus
(0)
(0)
(0)
K (t) dt. We can solve the zeroth-order terms
ξi , for example, Ii = ξi(0)in equation 4.14 self-consistently by finding the roots of the equations νi =(0)
k k (0)f (Ii + k τi J i νk ). In Figure 5, we compare this self-consistent mean-field
solution with numerical simulation of a single neuron as a function of constant input, I (0) , and self-inhibition strength J . We see that the approximation is good for small values of self-interaction J , but that the error increases
with | J | and I (0) , as expected given the analysis in the previous section.
To perform the linear response analysis, we require the first-order terms.
Differentiation of equation 4.12 with respect to stimulus perturbation, δξ ,yields (0)(0) δνi (t) = f  Ii + µi (δ Ii (t) + δµi (t))

 (0)
 (0)
(0) 
(0) 
= f  Ii + µi
K i (t − s)δξi (s) ds + f  Ii + µi
×

wik (t − s)δνk (s) ds.

(4.15)k

Now we can rewrite equation 4.15 by introducing vector notation asδν(t) =

(K(t − s)δξ (s) + w(t − s)δν(s)) ds,(4.16)

1218

T. Toyoizumi, K. Rad, and L. Paninski

Figure 6: Comparison of the analytically derived and simulated input-output
filters. (A) A simulated normalized linear input-output filter for I (0) = 2 (thindotted) and I (0) = 4 (thick-dotted) is compared with analytical results at I (0) = 2(thin-solid) and I (0) = 4 (thick-solid) of equation 4.18. Because the strengthof the self-inhibition component changes with the baseline input throughf  (I (0) + µ(0) ), the input filter is sharpened for high baseline input. Other parameters are set to J = 1, τ = 10 ms, K (t) = e −t/τ I (t) with τ I = 20 ms and σξ = 1 Hz.
All filter amplitudes are normalized by their maximum amplitudes. (B) Left panels show the comparison of analytically derived input-output filters G 11 (t) andG 12 (t) of equation 4.20 (solid lines) and simulated input-output filters (dottedlines) for I (0) = 2 (thin lines) and I (0) = 4 (thick lines). Right panels show neuron
1’s spatiotemporal input-output filter G̃ 1 (t, z) for I (0) = 2 (top) and I (0) = 4 (bottom). See the text for parameter values. (A color version of this figure is available
online at http://www.mitpressjournals.org/doi/suppl/10.1162/neco.2008.0408-757.)

(0)

(0)

(0)

where we defined [K(s)]i j = δi j f  (Ii + µi )K i (s) and [w(s)]i j = f  (Ii +(0)jµi )wi (s). Therefore, the Fourier transformation of equation 4.16 yieldsδ ν̂(ω) = K̂(ω)δ ξ̂ (ω) + ŵ(ω)δ ν̂(ω)= Ĝ(ω)δ ξ̂ (ω),(4.17)

e iωt [I − ŵ(ω)]−1 K̂(ω). Note that I is anwith a gain function G(t) = dω
2πidentity matrix here (not the input vector). We should note that a mathematically identical linear response filter was derived for a related model inGinzburg and Sompolinsky (1994). One difference is that in their model,
spinlike binary variables flip depending on input level with a fixed time
constant. This corresponds, in our model, to the case that all the synaptic
jtime constants τi are identical. Hence, the linear response filter in equation 4.17 is a natural generalization of their result.

Mean-Field Approximations for GLM Neurons

1219

We show in Figure 6A a simulated normalized receptive field of a single
neuron and the analytical result of equation 4.17. In this single-neuron case,where K is set to be a single exponential function with time constant τ I , forsimplicity, the Fourier transformations of the synaptic filter and the input
filter are ŵ(ω) = τ/(1 + iωτ ) and K̂ (ω) = τ I /(1 + iωτ I ), respectively. In thiscase, we obtain, after a Fourier inverse transformation by a Cauchy integral,

G(t) =

J f 2 τ− −(1/τ −J f  )tf−t/τ Ie+e,
1 + J f  τ−
1 + J f  τ−

(4.18)

where τ− = (1/τ I − 1/τ )−1 and f  = f  (I (0) + µ(0) ). We see that the linear
response may be described as a combination of the stimulus-filter term K
along with an additional term due to the recurrent filter w. Of course, inthe limit J → 0, we recover G → f  K . Note also that the input-output filter
changes with the baseline input through f  . In the case of the exponential
nonlinearity, f , and J < 0, the larger the input baseline, the sharper the
linear filter. This is because the effective strength of of the spike history
effect is modulated by the input-output nonlinearity and the baseline input
level (see, e.g., Bryant & Segundo, 1976, for some similar effects in an invitro physiological preparation).
Next, we calculate how interactions between neurons can change theinput-output filter. For concreteness, we consider two neurons that receive a spatiotemporal white noise stimulus ξ (z, ∞t) with unit variance.
Input to the neuron i is described by ξi (t) = −∞ dzL i (z)ξ (t, z), where2
L i (z) = √12π e −(z−zi ) /2 are spatial filters with means z1 = 1 and z2 = −1.
For simplicity, we assume that two neurons receive the same baseline of(0)(0)
input, I1 = I2 = I (0) with an identical temporal stimulus filter, K (t) =−t/τ Ie(t), and have symmetric synaptic interactions J 12 = J 21 = J = −1,
J 11 = J 22 = 0, and τik = τ I = τ = 10 ms. The derivative of the nonlinear(0)(0)
function is described by f i = f (Ii + µi ) for i = 1, 2. Proceeding as above,
we can calculate the spatiotemporal input-output filter of neuron i asG̃ i (t, z) = j G i j (t)L j (z). Similar to the previous single-neuron case, thefi
Fourier transformation of the temporal stimulus
filter is [ K̂(ω)]i j = δi j 1/τ +iω,
f i J i jwhile the synaptic filter is [ŵ(ω)]i j = 1/τ +iω , which yield the Fourier transformation of the temporal input-output filter,Ĝ(ω) = [I − ŵ(ω)]−1 K̂(ω)=

1/(1/τ + iω)
1 − ( f 1 f 2 J 2 )/(1/τ + iω)2

f 1f 1 f 2 J /(1/τ + iω)×f 1 f 2 J /(1/τ + iω)
f 2

.

(4.19)

1220

T. Toyoizumi, K. Rad, and L. Paninski

After a Fourier inverse transformation using Cauchy integration, we findG(t) = e

−t/τ

(t)


f 1 cosh( f 1 f 2 J t)  f 1 f 2 sinh( f 1 f 2 J t)
f 1 f 2 sinh( f 1 f 2 J t).
f 2 cosh( f 1 f 2 J t)
(4.20)

Figure 6B compares the analytically derived input-output filter of equation 4.20 with numerical simulation for different input baseline I (0) = 2 andI (0) = 4. Figure 6B, also shows the spatiotemporal input-output filter ofneuron 1, G̃ 1 (t, z) = G 11 (t)L 1 (z) + G 12 (t)L 2 (z), for the two input baselines.
Under large input baseline conditions, we see stronger interactions between
the two neurons and, hence, significant differences in the effective inputoutput filtering properties of the neuron in the two conditions. Thus, asemphasized above, the filtering properties of the network can change in a
stimulus-dependent way, without any changes in the system parameters
K , J , or w. (See, e.g., Pillow et al., 2008, for some related results concerning
the effects of interneuronal coupling terms on the filtering properties of
populations of primate retinal ganglion neurons.)
5 The Generalized Linear Point Process Model with Markov
RefractorinessIn section 4, we saw that we could derive simple approximate expressions
for many important quantities (firing rate, correlations) in the basic generalized linear point process model as long as the individual synaptic couplingterms | J | are not so large that the accuracy of the mean-field expansion iscompromised (see Figures 3 and 4).
In particular, the small | J | condition is likely acceptable for the multineuronal coupling terms, since in many cases, the inferred coupling
parameters—in motor cortex (Paninski et al., 2004; Truccolo et al., 2005)
and retina (Pillow et al., 2008), for example—have been empirically foundto be small (though larger network effects are found in hippocampus; seeHarris et al., 2003; Okatan et al., 2005). Similarly, we might expect the “slow”
self-inhibition terms in the GLM (the terms responsible for adaptation in
the firing rate, for example) to be relatively small as well (Pillow, Paninski,
Uzzell, Simoncelli, & Chichilnisky, 2005; Truccolo et al., 2005).
However, it is clear that the generalized linear model with small weights| J | will not be able to account for the strong refractoriness that is a fundamental feature of neural point-process dynamics at fine timescales: in the GLmodel, strong, brief inhibition is produced by large, negatively weighted,
and sharply decaying history effects wii (t), and these large wii (t) terms spoil

Mean-Field Approximations for GLM Neurons

1221

our mean-field approximation. Thus, it is clear that we need to adapt themethods introduced above in order to handle strong refractory effects.
One possible solution would be to take a discrete-time approach: insteadof modeling neural responses as point processes in continuous time, we
could simply model each cell’s spike count within a discrete-time bin as a
Bernoulli random variable whose rate νi (t) is an increasing function of ui (t),very much as in the continuous-time model described above (Chornoboy
et al., 1988; Okatan et al., 2005). This discrete-time Bernoulli model sharesmuch of the ease of fitting (including concavity of the log likelihood; Escola& Paninski, 2007) as the continuous-time model. In fact, the discrete-timemodel converges in a natural sense to the continuous-time model as the timebinwidth dt becomes small. The advantage of this discrete-time formalism isthat the firing rate νi (t) can never exceed 1/dt (since the Bernoulli probabilities νi (t) dt cannot exceed one), and therefore, in a crude sense, the discretetime model has a refractory effect of length dt. A mean-field analysis canbe developed for this discrete-time model, which exactly mirrors that developed above in the continuous-time model. We simply need to exchange
our ordinary differential equations for discrete-time difference equations.
Of course, this discrete-time approach is unsatisfactory in at least one
respect, since we would like, if possible, to model the firing behavior ofneurons down to a millisecond timescale (Berry & Meister, 1998; Keat,
Reinagel, Reid, & Meister, 2001; Pillow et al., 2005; Paninski et al., 2007), andthe discrete-time approach, by construction, ignores these fine timescales.
Thus, we introduce a model that allows us to incorporate strong refractory
effects directly. We set the rates,
λi (t) = f (ui (t))I (xi (t) = M),(5.1)

where f (.) and ui (t) are defined as in section 4 and I (.) is the indicatorfunction that takes the value 1 if the argument is true. We have introduced
an auxiliary refractory variable xi (t), which takes a discrete state from 1 to
M. We assume that this variable xi (t) is itself Markovian, with transition
rates


−1/τr


 1/τr
Wi (t) =  0



0

0

0

−1/τr

0

0

1/τr
..
.

−1/τr

0

Si (t)
...

0
0

..

.





;



(5.2)

−Si (t)that is, when xi (t) is in state M, it will transition to state 1 with each spike,and then transitions occur from state m to state m + 1 with rate 1/τr until xi (t) has reached the spiking state M once again. Refractoriness in this1222

T. Toyoizumi, K. Rad, and L. Paninski

model is enforced because the neuron is silent whenever xi (t) is in one ofthe postspike quiescent states xi (t) < M. It is easy to see that this is a kind of
inhomogeneous renewal model (the delay τd required for the state variable
xi to move from state 1 to first reach the active state M is a sum of (M − 1)independent and identically (i.i.d.) distributed exponential random variables of mean τr , and so τd is an i.i.d. gamma variable with parameters(M − 1, 1/τr )) and may be considered a special case of the “inhomogeneous Markov interval” model introduced by Kass and Ventura (2001). By
adjusting the number of states M and the rate 1/τr , we may adjust the refractory properties of the model: for example, M = 1 implies that we havean inhomogeneous Poisson model, while M = 2 provides an exponentiallydecaying relative refractory effect, and if we let M be large, with τr scaling
like τr ∼ τ/M, we obtain an absolute refractory effect of length τ .
Finally, it is easy to show that the additional Markov term in the definitionof the rate λi (t) does not have a negative impact on the estimation of the
GL model parameters J given spiking data; as discussed in Paninski (2004),
maximizing the likelihood in this model (or constructing an EM algorithm,
as in Escola & Paninski, 2007) requires that we maximize a nonnegatively
weighted version of the standard point-process log likelihood, and thisweighted log likelihood retains all of the concavity properties of the original
GL model. Thus, this new model requires just a single concave optimizationand is as easy to fit as the standard GL point process model. The advantage
is that a strong relative refractory effect is intrinsic to the model and does
not need to be enforced by a large, brief, negative wii (t) term. Instead, we
may use small adjustments to the J ii terms to fine-tune the model to matchthe short-time details of the observed interspike interval density.2 As wewill see below, this leads to much more accurate mean-field approximationsof the firing rates and correlation functions in these networks.
Let us define the probability pi (t) = (P(xi (t) = 1), . . . , P(xi (t) = M))Tthat the neuron i is in state xi (t) = 1, . . . , M. Note that the boldface characters denote vectors. The dynamics of pi is described byd pi (t)= Wi (t) pi (t),
dt

(5.3)

2 In this article, we mostly use M = 3 (except for Figure 7B, where the effect of abrupt
refractoriness is studied), because this is the smallest value of M that requires the following
mean-field formulation (we do not need a vector formulation for an M = 2 case because
there is only one degree of freedom that corresponds to the probability of being in theactive state). It is worth noting that Escola and Paninski (2007) discuss methods for
additionally estimating an optimal transition matrix Wi (t) via an EM algorithm. This
provides another method for adjusting the short-time details of the model’s responses. In
addition, we may extend many of the mean-field methods developed below to the caseof more general rate matrices Wi (t). However, for simplicity, we will not pursue these
extensions here.

Mean-Field Approximations for GLM Neurons

1223

where Wi (t) has the same elements as Wi (t) of equation 5.2 except for[Wi (t)]1M = −[Wi (t)] MM = E[Si (t) | xi (t) = M]
E[Si (t)I (xi (t) = M)]
E[I (xi (t) = M)]
νi (t)=,pi (t)=(5.4)

which describes the probability of neuron i firing at time t given xi (t) = M.
Note that we used νi (t) = E[Si (t)] = E[Si (t)I (xi (t) = M)] and wrote the lastcomponent of pi (t) as [ pi (t)] M = pi (t).
We also define [ pi j (t | t  )]m = E[I (xi (t) = m) | S j (t  )] as the probability ofneuron i being in state m given a spike of neuron j at time t  . In particular,
the last component is written as pi j (t | t  ) = [ pi j (t | t  )] M . The evolution ofpi j (t | t  ) is described by a similar equation,d pi j (t | t  )= Wi j (t | t  ) pi j (t | t  ),
dt

(5.5)

where the transition matrix Wi j (t | t  ) has the same elements as Wi (t) of
equation 5.2 except for[Wi j (t | t  )]1M = −[Wi j (t | t  )] MM = E[Si (t) | xi (t) = M, S j (t  )]E[Si (t)I (xi (t) = M) | S j (t  )]E[I (xi (t) = M) | S j (t  )]
νi j (t | t  ).
=pi j (t | t  )
=

(5.6)

Before turning to “mean-field” approximation, it is worth noting that wemay solve for the firing rates and correlations exactly in the special caseof no synaptic couplings, J = 0. (The analogous case in the standard GL
model is the inhomogeneous Poisson case, which is of course trivial.) We
can proceed by exploiting either the renewal nature of the model (Gerstner& Kistler, 2002), which leads to convolution or infinite-sum formulas forthe firing rate, or the Markovian structure, which leads to somewhat moreintuitive ordinary differential equations. We pursue the second approachhere.
In this case, the mean firing intensity is given by
νi (t) = E[ f (Ii (t))I (xi (t) = M)] = f (Ii (t)) pi (t).

(5.7)

1224

T. Toyoizumi, K. Rad, and L. PaninskiThis in turn gives, together with equation 5.4, [Wi (t)]1M = −[Wi (t)] MM =f (Ii (t)). Hence, we can find the dynamics of mean firing intensity by solving
equation 5.3. In the special case of constant input Ii , it is easy to solve thefixed point of pi (t); we find for the last component that pi = [1 + (M −1)τr f (Ii )]−1 .
The autocorrelation function in this case is similarly easy to derive. (Notethat the cross-covariance functions are zero in this J = 0 case, since thecross-coupling terms are set to zero.) In case of no couplings, J = 0, theconditional firing intensity is given by
νi j (t | t  ) = E[ f (Ii (t))I (xi (t) = M) | S j (t  )] = f (Ii (t)) pi j (t | t  )(5.8)

for t > t  . We find those corner elements of Wi j (t | t  ) as
[Wi j (t | t  )]1M = −[Wi j (t | t  )] MM = f (Ii (t))(5.9)

for t > t  . This means that Wi j (t | t  ) = Wi (t) for t > t  . Hence, pi j (t | t  ) follows the same differential equation as pi (t) in this J = 0 case but startingfrom the initial conditionpi j (t | t − ) = lim pi j (t | t− |→0|)

= δi j (1, 0, . . . , 0)T + (1 − δi j ) pi (t)(5.10)

because the state xi is reset to state 1 just after a given spike of neuron i.
In the following sections, we apply mean-field approximations for thefiring rates and correlations in the nonzero J case.
5.1 Mean-Field Approximation of the GL Model with Markov Refractoriness. Let us approximate the mean firing intensity and the crosscorrelations of GL model with Markov refractoriness assuming that thecontribution of individual synaptic coupling is small. Similar to the calculation without the Markov refractoriness, we assume a gaussian distributionof recurrent input Hi (t) given xi (t) = M to findνi (t) = E[ f (ui (t))I (xi (t) = M)]= E[ f (ui (t)) | xi (t) = M)] pi (t)≈ F (Ii (t) + µi (t), σi (t)) pi (t),(5.11)

with conditional mean µi (t) = E[Hi (t) | xi (t) = M] and conditional variance σi2 (t) = Var[Hi (t) | xi (t) = M]. Also, assuming a gaussian distributionof recurrent input Hi (t) given xi (t) = M and a spike of neuron j at time t  ,Mean-Field Approximations for GLM Neurons

1225

we findνi j (t | t  ) = E[ f (ui (t))I (xi (t) = M) | S j (t  )]= E[ f (ui (t)) | xi (t) = M, S j (t  )]E[I (xi (t) = M) | S j (t  )]≈ F (Ii (t) + µi j (t | t  ), σi j (t | t  )) pi j (t | t  )(5.12)

for t > t  with conditional mean µi j (t | t  ) = E[Hi (t) | xi (t) = M, S j (t  )] andvariance σi2j (t | t  ) = Var[Hi (t) | xi (t) = M, S j (t  )]. The dynamics of pi (t) andpi j (t | t  ) are described by equations 5.3 and 5.5, respectively.
Next we calculate the conditional means and variances of the recurrent
input. We findµi (t) =



J ik

k

µi j (t | t  ) =



J ik

k

σi2 (t) =βik (t),pi (t)
βikj (t | t  )J ik J ilk,lσi2j (t | t  ) =J ik J ilk,l,pi j (t | t  )Bikl (t),pi (t)

Biklj (t | t  )pi j (t | t  )(5.13)

with auxiliary variables defined as kβ i (t) m = E[h ik (t) | xi (t) = m] · [ pi (t)]m , k
β i j (t | t  ) m = E[h ik (t) | xi (t) = m, S j (t  )] · [ pi j (t | t  )]m , kl
B i (t) m = Cov[h ik (t), h li (t) | xi (t) = m] · [ pi (t)]m , kl
B i j (t | t  ) m = Cov[h ik (t), h li (t) | xi (t) = m, S j (t  )] · [ pi j (t | t  )]m .

(5.14)

Note that we use the final component of the above quantities without anindex (e.g., βik (t) = β ik (t) M ). After some calculation (see the appendix), we
obtainβ ik (t) =β ikj (t | t  ) =

 t
−∞

 t

−∞

ds e −(t−s)/τi pik (t | s)νk (s),kds e −(t−s)/τi pik j (t | s, t  )νk j (s | t  ),
k

1226

T. Toyoizumi, K. Rad, and L. PaninskiB ikl (t) m =

 t

 t

−∞

−

ds

−∞

 −(t−s)/τik −(t−s  )/τil

ds e



[ pikl (t | s, s  )]m φkl (s, s  )[ pik (t | s)]m [ pil (t | s  )]m
νk (s)νl (s  ) ,[ pi (t)]m(5.15)

where [ pikl (t | s, s  )]m = E[I (xi (t) = m) | Sk (s), Sl (s  )] and, particularly, thelast component is written as pikl (t | s, s  ) = [ pikl (t | s, s  )] M .
5.2 Solution in the Mean-Field Limit. Similar to the case without
Markov refractoriness, we first derive a self-consistent equation for themean firing intensity and cross-covariances in the mean-field limit Nc → ∞,where J ∼ 1/Nc and assuming an asynchronous state so that φi j ∼ 1/Nc .
Because all the cross-covariance functions vanish and σi2 (t) → 0 in the limit,
we find, from equations 5.11 to 5.13, for t > t  ,νi (t) → f (Ii (t) + µi (t)) pi (t),νi j (t | t  ) → f (Ii (t) + µi (t)) pi j (t | t  ),µi (t) =


k

J ik

βik (t),pi (t)(5.16)

while the cross-correlation function is given by equation 4.3. For the corner terms of the transition matrices, we find [Wi (t)]1M → f (Ii (t) + µi (t))
from equation 5.4 and [Wi j (t | t  )]1M → f (Ii (t) + µi (t)) for t > t  from equation 5.6. Hence, the evolution of pi (t) and pi j (t | t  ) is described bydpi (t)= Wi (t) pi (t),dt
dpi j (t | t  )= Wi (t) pi j (t | t  )dt(5.17)

for t > t  , where the initial condition of pi j (t | t  ) is given by equation 5.10.
Finally, from equation 5.15, the evolution of β i (t) is written as an ODE,d k
1
β (t) = − k + Wi (t) β ik (t) + pik (t | t − )νk (t),dt i
τi

(5.18)

where we can apply equation 5.10 to evaluate pik (t | t − ) in this limit. Figure 7
plots the mean firing intensity calculated from the above equations for various input, and the mean-field approximation provides good approximationof the time dependent firing intensities. Figure 8 compares the above approximation (Nc → ∞ limit) with the mean-field approximation, including
the finite size effect that we discuss in the next section. In the mean-field
limit, all the cross-correlations vanish and the autocorrelation functions

Mean-Field Approximations for GLM Neurons

1227

Figure 7: Comparison of the mean firing intensity obtained from direct simulation of a single neuron (solid), without recurrent input H and no Markov
refractory effect (dotted), and the mean-field equation in the limit of Nc → ∞(dashed). In all cases, the Markov refractory effect with τr = 2 ms/(M − 1)
was included. (A) Weak step input (baseline I = 2, peak I = 4) and refractoriness M = 3; J = −1, τ = 10 ms. Although the model includes strong
Markov refractory effect, the mean-field approximation is not abolished. (B)Strong step input (baseline I = 2, peak I = 8) and refractoriness M = 10;
J = −1, τ = 10 ms. The mean-field approximation also approximates very
well the transient oscillation caused by the the strong refractory effect andsteep current. (C) The mean-field approximation works for colored noise inI (t)= −I (t) + ξ (t) with E[ξ (t)] = 4, Cov[ξ (t), ξ (t  )] = 0.1δ(t − t  ), andput: τ I d dt
τ I = 20 ms; M = 3, J = −1, τ = 10 ms. The initial error is due to the preset initial values that deviate from the true ones. (D) Sinusoidal input input:
I (t) = 2 + sin(2π · 10 Hz · t); M = 3, J = −1, τ = 10 ms. The mean-field results
give good approximations. (A color version of this figure is available online at
http://www.mitpressjournals.org/doi/suppl/10.1162/neco.2008.04-08-757.)

are approximated by including the Markov refractory effect but not the
self-interaction terms J ii . We can also calculate the linar input-output filter
discussed in section 4.3 with this Markov refractoriness (see the appendix).
5.3 Estimating the Finite Size Effect. Similar to the case without
Markov refractoriness, we evaluate the mean firing intensity and crosscovariance functions up to 1/Nc terms. We find, as in the case without
Markov refractoriness, φi j (t, t  ) ∼ 1/Nc for i = j,  pik (t | s) ∼ 1/Nc fori = k, and the third-order covariances, such as  pikl (t | s, s  ), are of order1228

T. Toyoizumi, K. Rad, and L. Paninski

Figure 8: Comparison of spike cross-correlation functions of two neurons obtained by direct numerical simulation (solid), mean-field approximation in
the limit of Nc → ∞ (dashed), and with finite size correction terms (dotdashed). While the approximation in the limit of Nc → ∞ captures only thestrong Markov refractory effect, the weaker refractory effect caused by the
self-interaction terms is well captured with the finite size correction. The finite size correction provides a good approximation despite the strong refractoriness because of an explicit evaluation of the Markov refractory effect. Mean-field equations are iteratively solved over five iterations. (A) Two
neurons are directly coupled: parameters are set to J 11 = −1, J 22 = −1, J 12 =−0.5, J 21 = 1, τ11 = 10 ms, τ22 = 10 ms, τ12 = 20 ms, τ21 = 10 ms, M = 3, I1 =2, I2 = 4. (B) Two neurons are not directly connected but receive common input from a third neuron: parameters are set to J 13 = J 23 = 2, τ13 = τ23 = 10 ms,
M = 3, I1 = I2 = I3 = 2. (A color version of this figure is available online at
http://www.mitpressjournals.org/doi/suppl/10.1162/neco.2008.04-08-757.)

1/Nc2 under the asynchronous state if the contributions of individual synaptic weights are small, J ∼ 1/Nc . This implies σi2j (t | t  ) − σi2 (t) ∼ 1/Nc2 .
The dynamics of pi (t) and pi j (t | t  ) are described by equations 5.3 and
5.5, respectively. We now want to calculate µi (t), µi j (t | t  ), and σi2 (t) to thefirst order of 1/Nc . First, direct differentiation of β ik (t) in equation 5.15 yields
an ODE update equation (see the appendix):d k1
β (t) = − k + Wi (t) β ik (t) + pik (t | t − )νk (t)
dt i
τi

νi (t)−m βik (t) − αiik (t − | t) pi (t),pi (t)(5.19)

with m = (1, 0, . . . , 0, −1)T . The origin of the second term on the right-handside of equation 5.19 is due to the correlation of spike variable Sk (s) and theMean-Field Approximations for GLM Neurons

1229

refractory variable xi (t). Next, we need to evaluate β ikj (t | t  ) to the zerothorder of 1/Nc for k = i or k = j because J ∼ 1/Nc makes the contributionof them ∼ 1/Nc , and we need to evaluate β ikj (t | t  ) to the first order of 1/Nc
for k = i and k = j case. Except for the case i = j = k, we find (see theappendix)
β ikj (t | t  ) ≈ β ik (t) + αikj (t | t  ) pi j (t | t  ) − αik (t) pi (t)(5.20)

to the order described above. The precise evaluation of this term for i = j =k is much harder but, to a good approximation, equation 5.20 holds (see theappendix). Note that the approximation of the i = j = k term affects the order 1/Nc term of the autocovariance function, φii (t, t  ), but affects only theO(1/Nc2 ) terms of the mean firing intensities and cross-covariance functions.
Finally, dropping higher-order terms that result as O(1/Nc2 ) on the quantities in equation 5.13, we find (see the appendix)
B ikl (t) ≈ pi (t)Aikl (t).

(5.21)

Altogether, using equation 4.12, we obtain the following update equations to evaluate the finite size effect up to order 1/Nc :νi (t) = F (Ii (t) + µi (t), σi (t)) pi (t),νi j (t | t  ) = F (Ii (t) + µi j (t | t  ), σi (t)) pi j (t | t  ),
βik (t),pi (t)k 
β k (t) − αik (t) pi (t)J ik αikj (t | t  ) + i,µi j (t | t  ) =pi j (t | t  )kJ ik J il Aikl (t),σi2 (t) =µi (t) =

J ik

k,l

α k (t)d kαi (t) = − i k + νk (t),dtτi
αi j (t | t )d kα (t | t  ) = −+ νk j (t | t  ),dt i jτikkAkl (t)d kll(t − | t)νk (t),Ai (t) = − i kl + αilk (t + | t)νl (t) + αik
dt
τi
d k
1
β i (t) = − k + Wi (t) β ik (t) + pik (t | t − )νk (t)dtτiνi (t),−m βik (t) − αiik (t − | t) pi (t)pi (t)(5.22)

1230

T. Toyoizumi, K. Rad, and L. Paninski

and the dynamics of pi (t) and pi j (t | t  ) follows equations 5.3 and 5.5. As
is the case without the Markov refractoriness, calculating this finite size
effect is computationally expensive for time-dependent input, although inprinciple, this equation should work under that case as well. We plot themean firing intensities of two neurons and cross-correlations in Figure 8.
The finite size effect well captures the cross-correlation functions. It also
successfully captures a peak in the cross-correlation function if two neurons
are not directly connected but receive a common input from a third neuron(see Figure 8). Again, if the values of the above equations are initially set tothe solution in the limit of Nc → ∞ given by equations 5.16 to 5.18, only acouple of iterations of the above equations are sufficient for the evaluationof cross-correlation functions.
6 Discussion
We have introduced mean-field methods for analyzing the dynamics of acoupled population of neurons whose activity may be well approximated as
the output of a generalized linear point process model. Our approximations
for the mean time-varying firing rate and correlations in the population areexact in the mean-field limit (Nc → ∞ under J ∼ 1/Nc scaling), though wehave found numerically that the finite size correction to the mean-field equations is useful even for physiologically relevant connectivity strengths. The
approximations may be computed by solving a set of coupled ordinary differential equations. This approach is much more computationally tractable
than direct Monte Carlo sampling and may lead to greater analytical insight
into the behavior of the network. In addition, we have introduced a newmodel, the generalized linear point process model with Markovian refractoriness, that captures strong refractoriness, retains all of the easy-fitting
properties of the standard generalized linear model, and whose firing rate
dynamics are much more amenable to mean-field analysis than the standard modells. Note that most of our illustrations involved the simulationof only a couple of mutually connected neurons; this small-Nc case can beconsidered a kind of worst-case analysis, since with more asynchronous
neurons (larger Nc ), the central limit theorem becomes progressively applicable, the distribution of recurrent input approaches the assumed gaussian
distribution, and our approximations become more accurate.
This kind of mean-field analysis—or more generally, approaches for
reducing the complexity of the dynamics of networks of noisy, nonlinear
neurons—has a long history in computational neuroscience, as reviewed(Hertz et al., 1991; Gerstner & Kistler, 2002; Renart et al., 2003). While thepoint-process models we have discussed here are somewhat distinct fromthe noisy integrate-and-fire-type models that have been analyzed most extensively in this literature (e.g., Mattia & Del Giudice, 2002; Shriki, Hansel,& Sompolinsky, 2003; Moreno-Bote & Parga, 2004, 2006; Fourcaud-Trocme,
Hansel, van Vreeswijk, & Brunel, 2003; Chizhov & Graham, 2007; Doiron,Mean-Field Approximations for GLM Neurons

1231

Lindner, Longtin, Maler, & Bastian, 2004; Lindner, Doiron, & Longtin, 2005),
it is more important to stress the differences in the motivation of this workand previous work. Our main goal here was to provide analytical tools thatan experimentalist who has fit a point process model to his or her data (asin, e.g., Paninski, 2004; Truccolo et al., 2005; Okatan et al., 2005; Pillow et al.,2008) may use to understand the behavior of the network or single-cellmodels that have just been inferred. In particular, we were interested in
predicting, for example, the mean firing rate of a specific GLM network to
a novel arbitrary dynamic input stimulus. This contrasts with the literaturecited above, which has for the most part focused on the mean firing
rates, first-order response, and fixed-point stability properties behavior ofidealized, infinitely large populations (Mattia & Del Giudice, 2002, is an
exception here) with homogeneous membrane and connectivity properties.
In addition, most analytical studies of the input-output property of modelneurons with a dynamical input start from a Fokker-Planck formalism andrely on a direct numerical simulation of the partial differential equation(Nykamp & Tranchina, 2000; Knight, Omurtag, & Sirovich, 2000; Chizhov& Graham, 2007; this direct approach is feasible only in the case that the system dynamics may be reduced to a state space of dimension at most two or
so, and therefore does not apply in the networks studied here), or on linear
response theory that assumes a small time-dependent component relative
to the baseline component (Shriki et al., 2003; Fourcaud-Trocme et al., 2003;
Doiron et al., 2004; Lindner et al., 2005), or on some quasi-stationary assumption (Knight et al., 2000; Mattia & Del Giudice, 2002; Fourcaud-Trocme& Brunel, 2005; Moreno-Bote & Parga, 2004, 2006) restricting the applicable
input stimulus to slowly changing one or two simple step stimuli with
sufficiently long gaps between steps that the network may reach equilibrium. Thus, it is difficult to apply the insights developed in these previous
analyses directly to obtain a quantitative prediction of the dynamic firing
rates of a specific nonhomogeneous, nonsymmetric network.
We derived, in this article, the finite size correction to the mean-field
equation that captures cross-correlation functions between neurons. This
finite size effect originatesfrom the J i j ∼ 1/Nc scaling of synaptic strengths√
that guarantees O(1/ Nc ) fluctuation of inputs under asynchronous states.
Hence, the finite size effect described in this article is the contribution ofsmall fluctuation in the input around the mean. Another kind of finite size
effect has also been discussed in the literature. In a sparsely connected
network, correlations between two neurons disappear if Nc is sufficiently
smaller than the total number of neurons, N (Derrida, Gardner, & Zippelius,
1987). The finite Nc effect has been evaluated using stochastic Fokker-Planck
equations (Brunel & Hakim, 1999; Mattia & Del Giudice, 2002). One shouldnote, however, that the evaluation of this finite size effect for an experimentally estimated network structure is generally not straightforward because
the evaluation of a Fokker-Planck equation with many state-space variables is computationally hard. Finally, J i j ∼ 1/Nc is not the unique way to1232

T. Toyoizumi, K. Rad, and L. Paninski√
scale synapses. Under the balanced input assumption, J i j ∼ 1/ Nc yields
order one fluctuation of input even in the limit of Nc → ∞ (Sompolinksy,
Crisanti, & Sommers, 1988; van Vreeswijk & Sompolinsky, 1996, 1998). However, nontrivial cross-correlation functions cannot be captured under the
standard asynchronous assumption of mean-field analysis in the Nc →√∞limit. The calculation of the finite size effect of a network with J i j ∼ 1/ Nc
is not the scope of this article.
Two works that bear a stronger mathematical resemblance to ours are
Ginzburg and Sompolinsky (1994) and Meyer and van Vreeswijk (2002).
These authors applied mean-field approximations to coupled model neurons to study the mean firing rates and autocorrelation and cross-correlation
functions in the network. However, the two-state neuron model discussed
in Ginzburg and Sompolinsky (1994) (where the neurons flip between “active” and “inactive” firing states according to a Markov process with somefinite rate constant) is somewhat distinct from the point-process models
we have treated here (where the neuron remains in the “active” state fora negligible time—i.e., the neuron spikes instantaneously), and we havenot been able to translate these earlier results to the problems of interest in
this article. Renewal point spiking neuorns are analyzed in Meyer and van
Vreeswijk (2002), and their results are closer to ours. However, according totheir refractory model, one has to consider infinitely many refractory statesin the continuous time limit, whereas our spiking neuron model has M
refractory state; M can be as small as 2 while keeping the order 1 strong refractory effect. Accordingly, their cross-correlations should be evaluated byintegral equations that are computationally (and conceptually) somewhat
more involved, while the methods developed here require us to evaluate
just a couple of differential equations.
Finally, the effect of common input on the cross-correlation function waspreviously studied in a related model (Nykamp, 2007), using an expansionof the output f (.) nonlinearity arround J = 0. One major difference between
the analysis we have presented here and this previous work is that a simple
expansion around J = 0 does not lead to a good approximation of thefiring rate, even in the mean-field limit of Nc → ∞, because the recurrent
input changes the baseline firing rate (cf. Figure 5); thus, it is much moreaccurate to expand around the zeroth-order firing rate given by the rootsof equation 4.13, as discussed in section 4.3. (On the other hand, our meanfield approach does require a gaussian approximation that is known to beinaccurate in the case of large J terms.) It would be interesting to explore
whether the methods developed here could help lead to more accurateinference of the common-input effects discussed in Nykamp (2007).
There are many possible applications of this mean-field method for problems that require fast evaluation of mean firing rates and cross-correlation
functions. One example is to evaluate the information coded by spikingof a recurrently connected network about the input stimulus. This kind of

Mean-Field Approximations for GLM Neurons

1233

information calculation usually requires averaging over recent spike history (Toyoizumi, Aihara, & Amari, 2006), but the gaussian approximationof the input described here could greatly ease the computational complexity to evaluate the information of stimulus coded by the network. We hopeto pursue this direction in future work.
Appendix: Mathematical Details
In this appendix, we collect the details of the expansion analyses summarized in the main text.
k
A.1 Calculation of β ki , β kij , and B kl
i . In this appendix, we use Di s =e −(t−s)/τi ds to simplify some expressions. Direct evaluation of equation 5.14
yieldsk k
β i (t) m = E h ik (t) | xi (t) = m · [ pi (t)]m= E h ik (t)I (xi (t) = m)
 t
Dik s E[Sk (s)I (xi (t) = m)]
=
−∞

=

 t

−∞

Dik s [ pik (t | s)]m νk (s), kβ i j (t | t  ) m = E[h ik (t) | xi (t) = m, S j (t  )] · [ pi j (t | t  )]m(A.1)= E[h ik (t)I (xi (t) = m) | S j (t  )], t
Dik s E[Sk (s)I (xi (t) = m) | S j (t  )],
=
−∞

=

 t

−∞

Dik s [ pik j (t | s, t  )]m νk j (s | t  ),where [ pik (t | s)]m = E[I (xi (t) = m) | Sk (s)](xi (t) = m) | Sk (s), S j (t  )]. Similarly, we findand[ pik j (t | s, t  )]m = E[I
B ikl (t) m = Cov h ik (t), h li (t) | xi (t) = m · [ pi (t)]m = E h ik (t)h li (t) | xi (t) = m − E h ik (t) | xi (t) = m

×E h li (t) | xi (t) = m · [ pi (t)]m[β k (t)]m [β li (t  )]m= E h ik (t)h li (t)I (xi (t) = m) − i[ pi (t)]m1234

T. Toyoizumi, K. Rad, and L. Paninski

=

 t
−∞

Dik s

 t
−∞


Dil s  E[Sk (s)Sl (s  )I (xi (t) = m)][ pik (t | s)]m [ pil (t | s  )]m
νk (s)νl (s  )−[ pi (t)]m
 t
 t

Dik s
Dil s  [ pikl (t | s, s  )]m φkl (s, s  )
=
−∞

−∞


[ pik (t | s)]m [ pil (t | s  )]m
νk (s)νl (s  ) .
−[ pi (t)]m(A.2)

A.2 Approximation of β ki (t), β kij (t | t ) and B kl
i (t). In order to evaluate
the finite size effect—order 1/Nc terms of the mean firing intensity andcross-covariance functions—we need to evaluate µi (t) pi (t) = k J ik β ik (t),µi j (t | t  ) pi j (t | t  ) = k J ik β ikj (t | t  ), and σi2 (t) pi (t) = k,l J ik J il B ikl (t) to thefirst order of 1/Nc . Assuming that the synaptic strengths scale as J ∼ 1/Nc ,we need to evaluate
β ik (t) = t−∞

Dik s pik (t | s)νk (s)to the first order of 1/Nc ,β ikj (t | t  ) = t−∞

Dik s pik j (t | s, t  )νk j (s | t  )to the zeroth order of 1/Nc if k = i or k = j because J ∼ 1/Nc but to the firstorder of 1/Nc otherwise (in the following, we divide these into five casesand consider each case separately), and

B ikl (t) m =

 t
−∞

−

Dik s

 t
−∞


Dil s  [ pikl (t | s, s  )]m φkl (s, s  )[ pik (t | s)]m [ pil (t | s  )]m
νk (s)νl (s  )[ pi (t)]mto the first order of 1/Nc if i, j, k are all different, and to the zeroth order of
1/Nc if k = l, k = i, or l = i.
First, direct differentiation of β ik (t) yields
β k (t)d kβ i (t) = − i k + pik (t | t − )νk (t) +
dt
τi

 t
−∞

Dik s Wik (t | s) pik (t | s)νk (s)1= − k + Wi (t) β ik (t) + pik (t | t − )νk (t)τiMean-Field Approximations for GLM Neurons

+

 t
−∞

1235

Dik s Wik (t | s) pik (t | s)νk (s)1= − k + Wi (t) β ik (t) + pik (t | t − )νk (t)τi t
νi (t)pik (t | s) νk (s)Dik s νik (t | s) −+mpi (t)−∞
1= − k + Wi (t) β ik (t) + pik (t | t − )νk (t)τi tνk (s)+m
Dik s νki (s | t) −pik (t | s) νi (t)pi (t)−∞
1= − k + Wi (t) β ik (t) + pik (t | t − )νk (t)τi+ m αiik (t − | t) −βik (t)νi (t),pi (t)(A.3)

where m = (1, 0, . . . , 0, −1)T , and we used in the above calculation, fromequations 5.4 and 5.6,
Wik (t | s) pik (t | s) = [Wik (t | s) − Wi (t)] pik (t | s)=mνik (t | s)νi (t)−pik (t | s)pi (t)pik (t | s).

(A.4)

Next, to evaluate β ikj (t | t  ), we use  pikl (t | s, s  ) = [ pikl (t | s, s  ) − pi (t) − pik (t | s) −  pil (t | s  )] ∼ 1/Nc2 if i, k, l are all different;  pik (t | s) = pik (t |s) − pi (t) ∼ 1/Nc for i = k. As we discussed at the beginning of this section,in order to evaluate µi j (t | t  ) to the first order of 1/Nc , we have to consider
the number of combinations of indexes. For example, there are Nc terms of
β iii and Nc2 terms of β ii j . Hence, we consider the following five possibilities
for combining indexes:

r Case 1 (i = j, i = k); evaluated up to the first order of 1/Nc :β ikj (t | t  ) ≈
=

 t

−∞

 t

−∞

Dik s [ pi (t) +  pik (t | s) +  pi j (t | t  )]νk j (s | t  )Dik s [ pi (t)νk j (s | t  ) +  pik (t | s)(νk (s) + νk j (s | t  ))+ pi j (t | t  )(νk (s) + νk j (s | t  ))] t
Dik s [ pi (t)νk j (s | t  ) +  pik (t | s)νk (s)
≈
−∞

1236

T. Toyoizumi, K. Rad, and L. Paninski+ pi j (t | t  )νk (s)] t=
Dik s [ pi (t)νk j (s | t  ) + ( pik (t | s) − pi (t))νk (s)−∞+ pi j (t | t  )νk (s)]= pi (t)αikj (t | t  ) + (β ik (t) − pi (t)αik (t)) +  pi j (t | t  )αik (s)= β ik (t) + pi (t)αikj (t | t  ) +  pi j (t | t  )αik (s)≈ β ik (t) + pi j (t | t  )αikj (t | t  ) − pi (t)αik (s).

(A.5)

Note that in the second line, we used pik ∼ 1/Nc and νi j ∼ 1/Nc ,for example, and neglected higher-order terms such as pik νk j .
r Case 2 (i = j, k = i); evaluated up to the first order of 1/Nc : tβ iik (t | t  ) ≈Dik s [ pii (t | t  ) +  pik (t | s)]νki (s | t  )
−∞

=

 t

−∞

Dik s [ pii (t | t  )νki (s | t  ) + ( pik (t | s) − pi (t))νk (s)]= pii (t | t  )αiik (t | t  ) + (β ik (t) − pi (t)αik (t))= β ik (t) + pii (t | t  )αiik (t | t  ) − pi (t)αik (t).

(A.6)

r Case 3 (i = k, i = j); evaluated up to the zeroth order of 1/Nc :β ii j (t | t  ) ≈

 t

−∞

Dii s pii (t | s)νi (s)= β ii (t)≈ β ii (t) + pi j (t | t  )αii j (t | t  ) − pi (t)αii (s).

(A.7)

r Case 4 ( j = k, i = j); evaluated up to the zeroth order of 1/Nc :β i j (t | t  ) ≈
j

≈

 t

−∞

 t

−∞

Di s [ pi (t) +  pi j (t | s) +  pi j (t | t  )]ν j j (s | t  )jDi s pi (t)ν j j (s | t  )j= pi (t)αi j (t | t  )j≈ β i (t) + pi (t)αi j (t | t  ) − pi (t)αi (t).
j

j

j

(A.8)

r Case 5 (i = j = k); evaluated up to the first order of 1/Nc :β iii (t | t  ) = t−∞

Dii s piii (t | s, t  )νii (s | t  ).

(A.9)

Mean-Field Approximations for GLM Neurons

1237

This is a difficult situation to analyze precisely. However, unless theneuron i receives extremely strong input, the conditional firing intensity
νii (s | t  ) is almost zero if s and t  are as close as the refractory time constantof the neuron (M − 1)τr , so the contribution of this interval on the integral
is small. If s and t  are separated more than the refractory time constant, theinteraction between the spike at s and t  becomes small, that is, piii (t | s, t  ) ≈pi (t) +  pii (t | s) +  pii (t | t  ), and νii (s | t  ) = νii (s | t  ) − νi (s) becomes
small. Hence, to a good approximation,
β iii (t | t  ) ≈
=

 t
−∞

 t

−∞

Dii s [ pii (t | t  ) +  pii (t | s)]νii (s | t  )Dii s [ pii (t | t  )νii (s | t  )+ pii (t | s)νi (s)+  pii (t | s)νii (s | t  )]
 t
Dii s [ pii (t | t  )νii (s | t  ) +  pii (t | s)νi (s)]≈−∞= β ii (t) + pii (t | t  )αiii (t | t  ) − pi (t)αii (t).

(A.10)The validity of this intuitive approximation should be examined by numerical simulations. Note that β iii (t | t  ) contributes to the order 1/Nc term of the
autocovariance function φii (t, t  ) but contributes only to order 1/Nc2 termsof the mean firing intensity νi (t) and cross-covariance functions φi j (t, t  )with i = j.
Finally, we can rewrite B ikl as

B ikl (t) m =

 t
−∞

Dik s

 t
−∞


Dil s  [ pikl (t | s, s  )]m φkl (s, s  )[ pik (t | s)]m [ pil (t | s  )]m
νk (s)νl (s )−[ pi (t)]m

 t
 t
=
Dik s
Dil s  [ pikl (t | s, s  )]m φkl (s, s  )−∞−∞+ [ pikl (t | s, s  )]m −
=

 t
−∞

Dik s

 t
−∞

Dil s [ pik (t | s)]m [ pil (t | s  )]m[ pi (t)]mνk (s)νl (s  )[ pikl (t | s, s  )]m φkl (s, s  )+ [ pikl (t|s, s  )]m −[ pik (t|s)]m [ pil (t|s  )]m
νk (s)νl (s  ) ,[ pi (t)]m(A.11)

1238

T. Toyoizumi, K. Rad, and L. Paninski

where we used  pik (t | s) = pik (t | s) − pi (t) and  pikl (t | s, s  ) = pikl (t |s, s  ) − pi (t) −  pik (t | s) −  pil (t | s  ). We need to evaluate B ikl up to thefirst order of 1/Nc if i, k, l are all different and up to the zeroth order of 1/Nc
if k = l, k = i, or l = i. Up to these orders, B ikl can be approximated as
B ikl (t) ≈

 t
−∞

Dik s

 t
−∞

Dil s  pi (t)φkl (s, s  )= pi (t)Aikl (t).

(A.12)

Note that  pik ∼ 1/Nc for i = k,  pikl (t | s, s  ) ∼ 1/Nc2 if i, k, l are different,and  piik (t | s, s  ) ∼ 1/Nc for i = k.
A.3 Calculation of the Linear Filter in the Generalized Linear Modelwith Markov Refractoriness. We assume that the input is given byIi (t) =K i (t − s)ξi (s) ds(A.13)(0)with stimulus ξi = ξi + δξi (t), and calculate the linearized input-output(0)
filter for small stimulus fluctuation δξi (t) about a constiant baseline ξi .
In the limit of Nc → ∞, the mean-field equation is given by
νi (t) = f (Ii (t) + µi (t)) pi (t),µi (t) =



J ik

k

βik (t),pi (t)dpi (t) = Wi (t) pi (t),dtd k1
β (t) = − k + Wi (t) β ik (t) + pi (t)νk (t).
dt i
τi(A.14)(0)First, for a constant input Ii , we obtain(0)(0) (0)νi = f i pi ,
 β k (0)(0)
J ik i (0) ,
µi =
pi
k
(0)

(0)

0 = Wi pi ,

−1
1
(0)
(0)
(0) (0)pi νk ,
β ik = k − Wiτi(A.15)Mean-Field Approximations for GLM Neurons
(0)

(0)

1239

(0)

where f i = f (Ii + µi ) and the second to last equation is easily solvable(0)with respect to pi ; we find (0)pi m =








(0)

τr f i(0)

1 + (M − 1)τr f i

1




(0)
1 + (M − 1)τr f i(m = 1, . . . , M − 1)
(A.16)(m = M).

Next, the linear response to a small perturbation δ Ii (t) =s)δξi (s) ds is given by

K i (t −δνi (t) = f i (δ Ii (t) + δµi (t)) pi + f i δpi (t),


(0)
k

βik δpi (t)k δβi (t)Ji−
δµi (t) =
,
(0)
(0)
pi
( pi )2
k
(0)

(0)

(0)

d
(0)
(0)
δ pi (t) = δWi (t) pi + Wi δ pi (t),
dt
1
d
(0)
(0)
δβ k (t) = − k + Wi
δβ ik (t) + δWi (t)β ik
dt i
τi
(0)

(0)

+δ pi (t)νk + pi δνk (t),(A.17)

where f i (0) = f  (Ii0 + µi ) and δWi (t) = meTM [δνi (t)/ pi − νi δpi (t)/( pi )2 ].
Remember that m = (1, 0, . . . , 0, −1)T and e M = (0, . . . , 0, 1)T . Fourier transformation of the above equations gives
(0)

(0)

(0)

(0)

δ ν̂i (ω) = f i (δ Îi (ω) + δ µ̂i (ω)) pi + f i δ p̂i (ω),


(0)
k

βik δ p̂i (ω)
k δ β̂i (ω)δ µ̂i (ω) =
,
Ji
−
(0)
(0)
pi
( pi )2
k
(0)

(0)

(0)

δ p̂i (ω) = L̂ i (ω)mδ ν̂i (ω),
kkδ β̂ i (ω) = q̂ik (ω)δ ν̂i (ω) + R̂i (ω)δ ν̂k (ω),
(0)

(0)

(0)

(A.18)(0)where L̂ i (ω) = [iω + meTM νi / pi − Wi ]−1 , q̂ik (ω) = [iω + 1/τik − Wi ]−1k(0)(0) (0)(0)(0)(0)[mβik / pi − meTM L̂ i (ω)mβik νi /( pi )2 + L̂ i (ω)mνk ], and R̂i (ω) = [iω +
(0)

(0)

1/τik − Wi ]−1 pi . Note that i is the imaginary unit.
The above equation is still very complicated. So it is worth thinking
about the special case J ik = 0. In this case, the linear response is described

1240

T. Toyoizumi, K. Rad, and L. Paninskiby(0)δ ν̂i (ω) =f i (0) pi(0)
1 − f i eTM L̂ i (ω)m

δ Îi (ω).

(A.19)

The calculation of L̂ i (ω) is straightforward thanks to the bidiagonal nature
(0)

(0)

(0)

of [iω + meTM νi / pi − Wi ]. After direct matrix inversion, we find


1



 ρ

 ρ2

L̂ i (ω) = τr ρ 
 ..
 .

 M−2
 ρ
ρ M−1 η

0

0

1

0

ρ

1

···

..

.

···

ρ

1

···

ρ η

ρη

2


0
.. 

.



,




0
η

(A.20)

with ρ = 1/(1 + iωτr ) and η = 1 + 1/(iωτr ). Hence, the linear response issimplified as(0)δ ν̂i (ω) =f i (0) pi(0)

1 − f i τr ρη(ρ M−1 − 1)

δ Îi (ω).

(A.21)

In particular, if M = 2, we find that the Fourier inverse transformation of the(0)
f i (0) pi
has an analytical form,
gain function Ĝ(ω) = δ ν̂i (ω)/δ Îi (ω) =
(0)
1+
f
τ
/(1+iωτ
r
r)
i
(0)
(0)
that is, using pi = 1/(1 + τr f i ), we findG(t) =

f i (0)(0)
1 + τr f i


(0)
(0)
δ(t) − f i e −( fi +1/τr )t (t) .

(A.22)

Note that in the limit of τr → 0, the trivial gain function is G(t) → f i (0) δ(t).
Hence, we can see that the suppressive kernel is added to this instantaneous
gain function due to the Markov refractoriness. This linear response for
J ik = 0 case is a special case of a more general linear response for a renewal
neuron (Gerstner & Kistler, 2002). Note that in equation A.18, we discussed
a more general case and considered the interaction between the refractory
effect and spike interaction effect of a network of recurrently connected
neurons.
Acknowledgments
We thank E. Shea-Brown for helpful conversations. T.T. is supported by
the Robert Leet and Clara Guthrie Patterson Trust Postdoctoral Fellowship,

Mean-Field Approximations for GLM Neurons

1241

Bank of America, Trustee. L.P. is supported by an NSF CAREER award, an
Alfred P. Sloan Research Fellowship, the McKnight Scholar award, and NEI
grant EY018003.
References
Amit, D. J., & Tsodyks, M. V. (1991). Quantitative study of attractor neural networks
retrieving at low spike rates. I: Substrate—spikes, rates, and neuronal gain.
Network, 2, 259–273.
Berry, M., & Meister, M. (1998). Refractoriness and neural precision. J. Neurosci., 18,
2200–2211.
Brillinger, D. (1988). Maximum likelihood analysis of spike trains of interacting
nerve cells. Biological Cyberkinetics, 59, 189–200.
Brown, E., Kass, R., & Mitra, P. (2004). Multiple neural spike train data analysis:
State-of-the-art and future challenges. Nature Neuroscience, 7, 456–461.
Brunel, N., & Hakim, V. (1999). Fast global oscillations in networks of integrateand-fire neurons with low firing rates. Neural Comput., 11, 1621–1671.
Bryant, H., & Segundo, J. (1976). Spike initiation by transmembrane current: A
white-noise analysis. Journal of Physiology, 260, 279–314.
Chizhov, A. V., & Graham, L. J. (2007). Population model of hippocampal pyramidal
neurons, linking a refractory density approach to conductance-based neurons.
Phys. Rev. E Stat. Nonlin. Soft Matter Phys., 75, 011924.
Chornoboy, E., Schramm, L., & Karr, A. (1988). Maximum likelihood identificationof neural point process systems. Biological Cybernetics, 59, 265–275.
Derrida, B., Gardner, E., & Zippelius, A. (1987). An exactly solvable asymmetricneural network model. Europhys. Lett., 4, 167–173.
Doiron, B., Lindner, B., Longtin, A., Maler, L., & Bastian, J. (2004). Oscillatory
activity in electrosensory neurons increases with the spatial correlation of thestochastic input stimulus. Phys. Rev. Lett., 93(4), 048101.
Escola, S., & Paninski, L. (2007). Hidden Markov models applied toward the inferenceof neural states and the improved estimation of linear receptive fields. Manuscript
submitted for publication.
Fourcaud-Trocme, N., & Brunel, N. (2005). Dynamics of the instantaneous firing
rate in response to changes in input statistics. J. Comput. Neurosci., 18, 311–321.
Fourcaud-Trocme, N., Hansel, D., van Vreeswijk, C., & Brunel, N. (2003). How spike
generation mechanisms determine the neuronal response to fluctuating inputs.
J. Neurosci., 23, 11628–11640.
Gerstner, W., & Kistler, W. (2002). Spiking neuron models: Single neurons, populations,plasticity. Cambridge: Cambridge University Press.
Ginzburg, I., & Sompolinsky, H. (1994). Theory of correlations in stochastic neuralnetworks. Phys. Rev. E, 50(4), 3171–3191.
Harris, K., Csicsvari, J., Hirase, H., Dragoi, G., & Buzsaki, G. (2003). Organizationof cell assemblies in the hippocampus. Nature, 424, 552–556.
Hatsopoulos, N., Ojakangas, C., Paninski, L., & Donoghue, J. (1998). Information
about movement direction obtained by synchronous activity of motor corticalneurons. PNAS, 95, 15706–15711.

1242

T. Toyoizumi, K. Rad, and L. Paninski

Hertz, J., Krogh, A., & Palmer, R. G. (1991). Introduction to the theory of neuralcomputation. Redwood City, CA: Addison-Wesley.
Kappen, H. J., & Spanjers, J. J. (2000). Mean field theory for asymmetric neuralnetworks. Phys. Rev. E, 61(5 Pt B), 5658–5663.
Kass, R., & Ventura, V. (2001). A spike-train probability model. Neural Comput., 13,
1713–1720.
Keat, J., Reinagel, P., Reid, R., & Meister, M. (2001). Predicting every spike: A modelfor the responses of visual neurons. Neuron, 30, 803–817.
Knight, B., Omurtag, A., & Sirovich, L. (2000). The approach of a neuron population
firing rate to a new equilibrium: An exact theoretical result. Neural Comput., 12,1045–1055.
Kohn, A., & Smith, M. (2005). Stimulus dependence of neuronal correlation inprimary visual cortex of the macaque. Journal of Neuroscience, 25, 3661–3673.
Kulkarni, J., & Paninski, L. (2007). Common-input models for multiple neuralspike-train data. Network: Computation in Neural Systems, 18, 375–407.
Lindner, B., Doiron, B., & Longtin, A. (2005). Theory of oscillatory firing inducedby spatially correlated noise and delayed inhibitory feedback. Phys. Rev. E Stat.
Nonlin. Soft Matter Phys., 72(6 Pt 1), 061919.
Marmarelis, P., & Marmarelis, V. (1978). Analysis of physiological systems: The
white-noise approach. New York: Plenum Press.
Martignon, L., Deco, G., Laskey, K., Diamond, M., Freiwald, W., & Vaadia, E. (2000).
Neural coding: Higher-order temporal patterns in the neurostatistics of cell
assemblies. Neural Comput., 12, 2621–2653.
Mattia, M., & Del Giudice, P. (2002). Population dynamics of interacting spiking
neurons. Phys. Rev. E Stat. Nonlin. Soft Matter Phys., 66, 051917.
Meyer, C., & van Vreeswijk, C. (2002). Temporal correlations in stochastic networksof spiking neurons. Neural Comput., 14, 369–404.
Moreno-Bote, R., & Parga, N. (2004). Role of synaptic filtering on the firing responseof simple model neurons. Phys. Rev. Lett., 92, 028102.
Moreno-Bote, R., & Parga, N. (2006). Auto- and cross-correlograms for the spike
response of leaky integrate-and-fire neurons with slow synapses. Phys. Rev. Lett.,
96, 028101.
Nirenberg, S., Carcieri, S., Jacobs, A., & Latham, P. (2002). Retinal ganglion cells act
largely as independent encoders. Nature, 411, 698–701.
Nykamp, D. (2007). A mathematical framework for inferring connectivity inprobabilistic neuronal networks. Math. Biosci., 205, 204–251.
Nykamp, D., & Tranchina, D. (2000). A population density approach that facilitateslarge-scale modeling of neural networks: Analysis and an application toorientation tuning. Journal of Computational Neuroscience, 8, 19–50.
Okatan, M., Wilson, M., & Brown, E. (2005). Analyzing functional connectivity
using a network likelihood model of ensemble neural spiking activity. NeuralComput., 17, 1927–1961.
Oksendal, B. (2002). Stochastic differential equations. Berlin: Springer.
Oram, M., Hatsopoulos, N., Richmond, B., & Donoghue, J. (2001). Excess synchrony
in motor cortical neurons provides redundant direction information with thatfrom coarse temporal measures. Journal of Neurophysiology, 86, 1700–1716.

Mean-Field Approximations for GLM Neurons

1243

Paninski, L. (2004). Maximum likelihood estimation of cascade point-process neural
encoding models. Network: Computation in Neural Systems, 15, 243–262.
Paninski, L., Fellows, M., Shoham, S., Hatsopoulos, N., & Donoghue, J. (2004).
Superlinear population encoding of dynamic hand trajectory in primary motor
cortex. J. Neurosci., 24, 8551–8561.
Paninski, L., Pillow, J., & Lewi, J. (2007). Statistical models for neural encoding,decoding, and optimal stimulus design. In P. Cisek, T. Drew, & J. Kalaska (Eds.),
Computational neuroscience: Progress in brain research. Dordrecht: Elsevier.
Pillow, J., Paninski, L., Uzzell, V., Simoncelli, E., & Chichilnisky, E. (2005). Prediction
and decoding of retinal ganglion cell responses with a probabilistic spiking
model. Journal of Neuroscience, 25, 11003–11013.
Pillow, J. W., Shlens, J., Paninski, L., Sher, A., Litke, A. M., & Chichilnisky, et al.
(2008). Spatio-temporal correlations and visual signalling in a complete neuronal
population. Nature, 454, 995–999.
Renart, A., Brunel, N., & Wang, X.-J. (2003). Computational neuroscience: A
comprehensive approach. Boca Raton, FL: CRC Press.
Riehle, A., Grün, S., Diesmann, M., & Aertsen, A. (1997). Spike synchronization andrate modulation differentially involved in motor cortical function. Science, 278,
1950–1953.
Rieke, F., Warland, D., de Ruyter van Steveninck, R., & Bialek, W. (1997). Spikes:
Exploring the neural code. Cambridge, MA: MIT Press.
Rigat, F., de Gunst, M., & van Pelt, J. (2006). Bayesian modelling and analysis ofspatio-temporal neuronal networks. Bayesian Analysis, 1, 733–764.
Schnitzer, M., & Meister, M. (2003). Multineuronal firing patterns in the signal fromeye to brain. Neuron, 37, 499–511.
Shriki, O., Hansel, D., & Sompolinsky, H. (2003). Rate models for conductance-based
cortical neuronal networks. Neural Comput., 15, 1809–1841.
Snyder, D., & Miller, M. (1991). Random point processes in time and space. Berlin:
Springer-Verlag.
Sompolinksy, H., Crisanti, A., & Sommers, H. (1988). Chaos in random neuralnetworks. Physical Review Letters, 61, 259–262.
Toyoizumi, T., Aihara, K., & Amari, S. (2006). Fisher information for spike-based
population decoding. Phys. Rev. Lett., 97, 098102.
Truccolo, W., Eden, U., Fellows, M., Donoghue, J., & Brown, E. (2005). A point
process framework for relating neural spiking activity to spiking history, neuralensemble and extrinsic covariate effects. Journal of Neurophysiology, 93, 1074–1089.
van Vreeswijk, C., & Sompolinsky, H. (1996). Chaos in neuronal networks withbalanced excitatory and inhibitory activity. Science, 274, 1724–1726.
van Vreeswijk, C., & Sompolinsky, H. (1998). Chaotic balanced state in a model ofcortical circuits. Neural Comput., 10, 1321–1371.
Wilson, H. R., & Cowan, J. D. (1972). Excitatory and inhibitory interactions inlocalized populations of model neurons. Biophys. J., 12(1), 1–24.

Received April 16, 2008; accepted September 18, 2008.