Notes on Generalized Linear Models of Neurons
Jonathon Shlens∗
Google
Research
Mountain View, CA 94043
(Dated: April 9, 2014, version 1.51)
Experimental neuroscience increasingly requires tractable models for analyzing and predicting the behavior of neurons and networks. The generalized linear model (GLM) is an increasingly popular statistical framework for analyzing neural data that is flexible, exhibits rich dynamic behavior and is

arXiv:1404.1999v1 [] 8 Apr 2014

computationally tractable (Paninski, 2004; Pillow et al., 2008; Truccolo et al., 2005). What follows is a
brief summary of the primary equations governing the application of GLM’s to spike trains with a few
sentences linking this work to the larger statistical literature. Latter sections include extensions of a basic GLM to model spatio-temporal receptive fields as well as network activity in an arbitrary numbers
of neurons.

The generalized linear model (GLM) is a powerful framework for modeling statistical relationships in complex data
sets. One application of the GLM is to relate the activity of a neuron (or network) to a sensory stimulus and surrounding network dynamics. A GLM of a neuron requires a probabilistic model of spiking activity paired with a
function relating extrinsic factors (e.g. preceding stimulus, refractory dynamics, network activity) to a neuron’s internal drive. In what follows, the former is provided by a Poisson likelihood while the latter is provided by a (cleverly
constructed) conditional intensity. Finally we discuss extensions of the GLM for modeling spatio-temporal receptive
fields and correlated network activity.1

I. POISSON LIKELIHOOD

A common starting point for modeling a neuron is the idealization of a Poisson process (Dayan and Abbott, 2001).
If y is the observed number of spikes, the probability of such an event is given by a Poisson distribution
P (y|λ) =

(λ∆)y
exp(−λ∆)
y!

(1)

where λ∆ is the expected number of spikes in a small unit of time ∆ and λ is the intensity of the Poisson process.
A spike train Y is defined as a vector of spike counts {yt } binned at a brief time resolution ∆ and indexed by time
t. The likelihood that the spike train arises from a time-varying (inhomogenous) Poisson process with a conditional

∗ Electronic
1

address: jonathon.shlens@gmail.com
This manuscript is based on notes by Jonathan Pillow and Liam Paninski and conversations with Eero Simoncelli and EJ Chichilnisky. The
manuscript assumes that the reader is familiar with sensory physiology and has a working knowledge of probability, linear algebra and
calculus. Bold letters k, h, s are column vectors while capitalized bold letters X are matrices.

2
(time-varying) intensity λt is the product of independent observations,
P (Y |θ) =

Y (λt ∆)yt
yt !

t

exp(−λt ∆)

where θ = {λt } refers to the collection of intensities over the entire spike train. The log-likelihood of observing the
entire spike train is
log P (Y |θ) =

X

yt log λt +

X

t

yt log ∆ −

X

t

log yt ! − ∆

t

X

λt

t

In a statistical setting our goal is to infer the hidden variables θ from the observed spike train Y . Many criterion exist
for judging the quality of values for θ (Kay, 1993; Moon and Stirling, 1999). Selecting a particular set of values for θ
that maximize the likelihood (or log-likelihood) is one such criterion that is intuitive and and relatively easy to calculate. The goal of this manuscript is to calculate the maximum likelihood estimate of hidden variables by optimizing
(or maximizing) the log-likelihood.
In our setting, the spike train has been observed and is therefore fixed. The log-likelihood is solely a function of the
(unknown) conditional intensities θ. To emphasize our application of likelihood, we relabel L(θ) ≡ log P (Y |θ) and
group the terms independent of θ as some arbitrary normalization constant c,
L(θ) =

X

yt log λt − ∆

t

X

λt + c.

t

In practice, models of neurons are at fine temporal resolutions, thus yt is effectively binary, if not largely sparse. The
binary representation of yt simplifies the log-likelihood further
L(θ) =

X
t=spike

log λt − ∆

X

λt .

(2)

t

In the final form of the log-likelihood, we ignore the normalization constant c because the goal of the optimization is
to determine the hidden variables θ (not the absolute value of the likelihood).

II. GENERALIZED LINEAR MODEL

A GLM can attribute the variability of a spike train to a rich number of factors including the stimulus, spike history
and network dynamics. The GLM extends the maximum likelihood procedure to more interesting variables than
λt by positing that the conditional intensity is directly related to biophysically interesting parameters. For example,
pretend there is some parameter k, a neuron’s receptive field. If we posit that the receptive field is linearly related to
some known (and invertible) function of λt , then in principle it is just as easy to infer k as it is to infer λt . By positing
a linear relationship, the stimulus filter is quite simple to estimate, but surprisingly provides a rich structure for
modeling the neural response.

3
+

A

B
nonlinearity

stim filter (k)

+

stochastic
spiking

neuron 1

ex

coupling
filters (h ij )

post-spike filter (h)

+
neuron 2

FIG. 1 Diagram of a generalized linear model of a single neuron (a) and a network of neurons (b).

More precisely, the goal of a single neuron GLM is to predict the current number of spikes yt using the recent spiking
history and the preceding stimulus. Let xt = (xt−τ , . . . , xt−1 ) represent the vector of preceding stimuli up to but not
including time t. Let yt = (yt−τ , . . . , yt−1 ) be a vector of preceding spike counts up to but not including time t.
t-1

time
xt
yt

t

-0.45 0.02 0.26 0.33 -0.03
0

1

0

0

0

yt

Importantly, note that xt and yt are fixed known variables but yt is an unknown random variable. We posit that yt is
distributed according to a Poisson distribution whose conditional intensity λt is related to the stimulus and previous
spiking history by
λt = f (k · xt + h · yt + µ)

(3)

where k is a stimulus filter of the neuron (i.e. receptive field), h is a post-spike filter to account for spike history
dynamics (e.g. refractoriness, bursting, etc.) and µ is a constant bias to match the neuron’s firing rate. Figure 1a
provides a diagram of Equation 3. Each term of the conditional intensity increases or decreases the average firing
rate depending on the preceding stimulus (term 1) and the spiking history of the neuron (term 2). f is an arbitrary
invertible function termed a link function. One possibilty is to select f = exp(·), termed the natural link function,
because it conveniently simplifies the likelihood calculation and the interpretation of individual terms (McCullagh
and Nelder, 1989). In this case, a stimulus preceding time t that closely matches the filter k increases the average
spike rate by a multiplicative gain exp(k · xt ). Likewise, if a spike occurs just prior to time t, the convolution of the
spike occurrence (a delta function) with the post-spike filter diagrammed in Figure 1a decreases the probability of a
spike by factor exp(h · yt ) to mimic a refractory period.
It is worth noting at this point that the name generalized linear model refers to the requirement that f −1 (λ) is linearly
related to the parameters of interest, θ = {k, h, µ}. Combining equations 2 and 3, the probability of observing the

4
complete spike train Y is
X

L(θ) =

(k · xt + h · yt + µ) − ∆

X

exp(k · xt + h · yt + µ)

(4)

t

t=spike

Importantly, the likelihood is concave everywhere in the parameter space θ = {k, h, µ}, thus no local maxima exist
and ascending the gradient leads to a single, unique global maximum. To calculate the maximum likelihood estimate
of θ one must calculate the gradient and Hessian of the likelihood with respect to each variable,

Gradient
∂
∂k L(θ)
∂
∂h L(θ)
∂
∂µ L(θ)

=
=
=

Hessian
2

∂
∂k2 L(θ)
P
P
∂2
t=spike yt − ∆
t yt λt ∂h2 L(θ)
P
∂2
nsp − ∆ t λt
∂µ2 L(θ)

P

t=spike

xt − ∆

P

t

xt λt

∂2
∂k∂µ L(θ)
P
∂2
−∆ t yt ytT λt ∂h∂µ
L(θ)
P
∂2
−∆ t λt
∂k∂h L(θ)

= −∆
=
=

P

t

xt xTt λt

= −∆

P

xt λt

= −∆

P

yt λt

= −∆

P

xt ytT λt

t

t
t

where nsp denotes the total number of spikes in the spike train. These series of equations are solved by standard
gradient-ascent algorithms (e.g. Newton-Raphson), permitting an estimate of the model parameters for a given spike
train and stimulus movie. In practice, a model with 30 parameters can be fit with roughly 200 spikes.

III. RECEPTIVE FIELDS IN SPACE AND TIME

Generalizing the GLM to multiple spatial dimensions, the conditional intensity retains the same form2
s
X
ki · xi + h · y + µ)
λ = exp(

(5)

i

with the inner product of the stimulus filter and the stimulus summed over all s spatial locations indexed by i. In
practice though the number of parameters in the stimulus filter grows quadratically as the number of spatial locations and samples in the receptive field’s integration period, making any estimate computationally slow and potentially intractable. To exploit matrix algebra in further equations, define K ≡ [k1 k2 . . . ks ]T and X ≡ [x1 x2 . . . xs ]T .
Rather than letting the number of parameters grow quadratically, one simplification is to posit that K is space-time
separable, meaning that K can be written as an outer product,
K = tsT
where the vectors s and t independently capture the spatial and temporal components of the receptive field, respectively. Rewriting the conditional intensity with matrix algebra,
λ = exp(sT Xt + h · y + µ)

2

For the remainder of this manuscript we drop the subscript t on λt , xt and yt for notational simplicity.

(6)

5
the model is not linear, but quadratic, in parameter space θ = {s, t, h, µ} and the likelihood surface is not concave
everywhere (i.e. local maxima exist). In particular, the lack of concavity arises from the new Hessian term measuring
the dependency between the space and time filters.
X
X
X
∂2
L(θ) =
X−∆
Xλt − ∆
XtsT Xλt
∂s∂t
t
t

(7)

t=spike

Ascending the gradient, however, produces a good approximation of the standard model but with a notable reduction in the number of parameters. The gradient formulae are all the same as Section II if one replaces xt with Xt for
the derivatives associated with s, and sX for the derivatives associated with t.

IV. NETWORKS OF NEURONS

The most significant feature of networks of neurons is that the activity of one neuron significantly influences the activity of surrounding neurons. One common example of correlated activity is the observation of synchronized firing,
when two or more neurons fire nearly simultaneously more often then expected by chance. Importantly, correlated
activity can be independent of shared stimulus drive and instead reflect an underlying biophysical mechanisms (e.g.
common input, electrical or synaptic coupling, etc.).
A naive implementation of the GLM (Equation 3) for two neurons with conditional intensities λ1 and λ2 would
fail to capture synchrony or any correlated activity because λ1 and λ2 are independent3 . A simple extension of the
conditional intensity can mimic correlated activity such as synchrony (Pillow et al., 2008). In particular, one can add
post-spike filters hij to permit spikes from one neuron to influence the firing rate of another neuron (see Figure 1b).
The conditional intensity of neuron i is then
λi = exp(ki · x +

X

hij · yj + µi )

(8)

j

where we have added the subscript i throughout to label neuron i. The term

P

j

hij yj sums over post-spike activity

received by neuron i, whether internal dynamics (j = i) or activity from other neurons (j 6= i).
The complete likelihood of a population of neurons is
L(θ) =

X
i,t=spike

(ki · x +

X
j

hij · yj + µi ) − ∆

X
i,t

exp(ki · x +

X

hij · yj + µi )

(9)

j

where θ = {ki , hij , µi } is the parameter set. Note again that the subscript for time t is dropped from the stimulus

3

More precisely, the neurons modeled by λ1 and λ2 would not be independent but conditionally independent (Schneidman et al., 2003). The latter
distinction means that both neurons would be independent given that the stimulus is held fixed. This distinction is important because one
could select a stimulus which drives both neurons and attribute the apparent correlated activity to an underlying mechanism as opposed to the
choice of stimulus.

6
and spike trains for clarity – although the likelihood sums over this implicit subscript. The gradient of the likelihood
of each post-spike coupling term hij must be generalized
X
X
∂
L(θ) =
yj − ∆
yj λi
∂hij
t
t=spike

where care must be taken to pair the appropriate spike train yj with conditional intensity λi . Most of the terms of
the Hessian are preserved from the single neuron case (Section II) except for the terms associated with hij :
X
∂2
L(θ) = −∆
yj yjT λi ,
2
∂hij
t

X
∂2
L(θ) = −∆
xi yjT λi ,
∂ki ∂hij
t

X
∂2
L(θ) = −∆
yjT λi
∂µ∂hij
t

Again, care must be taken to pair the appropriate spike train yj with the appropriate conditional intensity λi . Note
that Hessian terms between different neurons (e.g. ki , kj , for j 6= i) are zero due to the linearity of the conditional intensity. Because the cross terms of the Hessian between neurons are zero, one can fit a parameter set θi = {ki , hij , µi }
for each neuron i independently. In practice, this independence vastly simplifies the fitting procedure because each
neuron can be fit in parallel.

V. CONCLUSIONS

This manuscript provided a very brief summary of analyzing neural data with generalized linear models (GLM’s).
The GLM is a flexible framework that allows one to attribute spiking activity to arbitrary phenomena – whether
network activity, sensory stimuli or other extrinsic factors – and then estimate parameters efficiently using standard
maximum likelihood techniques. Thus far, GLM’s (or related ideas) have been used to model network activity in
the retina (Pillow et al., 2008), motor cortex (Truccolo et al., 2005), visual cortex (Koepsell et al., 2008), hippocampus (Brown et al., 2001) as well as devise new strategies for the design of experiments (Paninski et al., 2007), quantify
stimulus information in spike trains (Pillow and Paninski, 2007) and help a lowly graduate student complete a Ph.D.
thesis (Shlens, 2007).

References

Brown, E., D. Nguyen, L. Frank, M. Wilson, and V. Solo, 2001, Proceedings of the National Academy of Sciences 98(21), 12261.
Dayan, P., and L. Abbott, 2001, Theoretical Neuroscience: Computational and Mathematical Modeling of Neural Systems (The MIT Press).
Kay, S. M., 1993, Fundamentals of Statistical Signal Processing, Volume I: Estimation Theory (Prentice Hall PTR).
Koepsell, K., T. Blanche, N. Swindale, and B. Olshausen, 2008, in Computational and Systems Neuroscience (Cosyne) 2008 (Salt Lake
City, Utah).
McCullagh, P., and J. A. Nelder, 1989, Generalized Linear Models, Second Edition (Monographs on Statistics and Applied Probability)
(Chapman & Hall/CRC).
Moon, T. K., and W. C. Stirling, 1999, Mathematical Methods and Algorithms for Signal Processing (Prentice Hall).
Paninski, L., 2004, Network 15(4), 243.

7
Paninski, L., J. Pillow, and J. Lewi, 2007, in Computational Neuroscience: Theoretical Insights into Brain Function, edited by P. Cisek,
T. Drew, and J. Kalaska (Elsevier, Amsterdam).
Pillow, J., and L. Paninski, 2007, Neural Computation , submitted.
Pillow, J., J. Shlens, L. Paninski, A. Sher, A. Litke, E. Chichilnisky, and E. Simoncelli, 2008, Nature , accepted.
Schneidman, E., W. Bialek, and M. Berry, 2003, J Neurosci 23(37), 11539.
Shlens, J., 2007, Synchrony and concerted activity in the neural code of the retina, Ph.D. thesis, University of California, San Diego.
Truccolo, W., U. Eden, M. Fellows, J. Donoghue, and E. Brown, 2005, J Neurophysiol 93(2), 1074.

