11 Encoding and Decoding with Stochastic Neuron models

11.2 Encoding Models in Systems Neuroscience

Generalized integrate-and-fire models have been used not only for spike prediction of neurons in brain slices, but also for measurements in systems neuroscience, i.e., in the intact brain driven by sensory stimuli or engaged in a behavioral task. Traditionally, electrophysiological measurements in vivo have been performed with extracellular electrodes or multi-electrode probes. With extracellular recording devices, the presence of spikes emitted by one or several neurons can be detected, but the membrane potential of the neuron is unknown. Therefore, in the following we aim at predicting spikes in extracellular recordings from a single neuron, or in groups of connected neurons.

A B
Fig. 11.5: The encoding problem in visual neuroscience. A. A stimulus is presented on a screen while a spike train is recorded from an area in visual cortex. B. Models designed to predict the spike train first filter the stimulus 𝐱{\bf x} with a spatial filter 𝐤{\bf k} (linear processing step), pass the result u=𝐤𝐱u={\bf k}\cdot{\bf x} through a nonlinearity ff and then generate spikes stochastically with Poisson statistics. The main difference between an Linear-Nonlinear-Poisson (LNP, top) and a soft-threshold generalized integrate-and-fire model (GLM, bottom) is the presence of spike-triggered currents η~(s)\tilde{\eta}(s) in the latter.

11.2.1 Receptive fields and Linear-Nonlinear Poisson Model

The linear properties of a simple neuron in primary visual cortex can be identified with its receptive field, i.e., the small region of visual space in which the neuron is responsive to stimuli (cf. Chapters 1 and 12). Receptive fields as linear filters have been analyzed in a wide variety of experimental preparations.

Experimentally, the receptive field of a simple cell in visual cortex can be determined by presenting a random sequence of spots of lights on a gray screen while the animal is watching the screen. In a very limited region of the screen, the spot of light leads to an increase in the probability of firing of the cell, in an adjacent small region to a decrease. The spatial arrangement of these regions defines the spatial receptive field of the cell and can be visualized as a two-dimensional spatial linear filter (Fig. 11.5B).

Instead of a two-dimensional notation of screen coordinates, we choose in the following a vector notation where we label all pixels with a single index kk. For example, on a screen with 256×256256\times 256 pixels we have 1kK1\leq k\leq K with K=65536K=65^{\prime}536. A full image corresponds to a vector 𝒙=(x1,xK)\mbox{\boldmath\(x\)}=(x_{1},\dots x_{K}) while a single spot of light corresponds to a vector with all components equal to zero except one (Fig. 11.6A).

The spatial receptive field of a neuron is a vector 𝒌k of the same dimensionality as 𝒙x. The response of the neuron to an arbitrary spatial stimulus 𝒙x depends on the total drive 𝒌𝒙t\mbox{\boldmath\(k\)}\cdot\mbox{\boldmath\(x\)}_{t}, i.e., the similarity between the stimulus and the spatial filter.

More generally, the receptive field filter 𝒌k can be described not only by a spatial component, but also by a temporal component: an input 100ms ago has less influence on the spiking probability now than an input 30ms ago. In other words, the scalar product 𝒌𝒙t\mbox{\boldmath\(k\)}\cdot\mbox{\boldmath\(x\)}_{t} is a short-hand notation for integration over space as well as over time. Such a filter 𝒌k is called a spatio-temporal receptive field.

In the linear-nonlinear-Poisson (LNP) model, one assumes that spike trains are produced by an inhomogeneous Poisson process with rate

ρ(t)=f(𝒌𝒙t)\rho(t)=f(\mbox{\boldmath\(k\)}\cdot\mbox{\boldmath\(x\)}_{t}) (11.5)

given by a cascade of two simple steps (Fig. 11.5B). The linear stage, 𝒌𝒙t\mbox{\boldmath\(k\)}\cdot\mbox{\boldmath\(x\)}_{t}, is a linear projection of 𝒙t\mbox{\boldmath\(x\)}_{t}, the (vector) stimulus at time tt, onto the receptive field 𝒌k; this linear stage is then followed by a simple scalar nonlinearity f(.)f(.) which shapes the output (and in particular enforces the non-negativity of the output firing rate ρ(t)\rho(t)). A great deal of the systems neuroscience literature concerns the quantification of the receptive field parameters 𝒌k.

Note that the LNP model neglects the spike history effects that are the hallmark of the SRM and the GLM – otherwise the two models are surprisingly similar; cf. Fig. 11.5B. Therefore, a LNP model cannot account for refractoriness or adaptation, while a GLM in the form of a generalized soft-threshold integrate-and-fire model does. The question arises whether a model with spike-history effects yields a better performance than the standard LNP model.

A B
Fig. 11.6: Spatial receptive field measurement. A. While the animal focuses on the center (star), light dots are presented at random positions on a gray screen; in the present trial, pixel 19 lights up. The input is denoted as a vector 𝒙t\mbox{\boldmath\(x\)}_{t}. B. Schematic of input matrix XX. Matrix representing a sparse input, such as a single spot of light. Rows of the matrix correspond to different trials, marked by the observation time tt.

Both models, LNP and GLM, can be fitted using the methods discussed in Ch. 10. For example, the two models have been compared on a data set where retinal ganglion cells have been driven by full-field light stimulus, i.e., the stimulus did not have any spatial structure (398). Prediction performance ranged in similar values as for cortical neurons driven by intracellular current injection, with up to 90% of the PSTH variance predicted in some cases. LNP models in this context have significantly worse prediction accuracy; in particular, LNP models greatly overestimate the variance of the predicted spiking responses. See Fig. 11.7 for an example.

Example: Detour on reverse correlation for receptive field estimation

Reverse correlation measurements are an experimental procedure based on spike-triggered averaging (113; 94). Stimuli 𝒙x are drawn from some statistical ensemble and presented one after the other. Each time the neuron elicits a spike, the stimulus 𝒙x presented just before the firing is recorded. The reverse correlation filter is the mean of all inputs that have triggered a spike

𝒙RevCorr=𝒙spike=tnt𝒙ttnt,\mbox{\boldmath\(x\)}_{\rm RevCorr}=\langle\mbox{\boldmath\(x\)}\rangle_{\rm spike% }=\frac{\sum_{t}n_{t}\mbox{\boldmath\(x\)}_{t}}{\sum_{t}n_{t}}\,, (11.6)

where ntn_{t} is the spike count in trial tt. Loosely speaking, the reverse correlation technique finds the typical stimulus that causes a spike. In order to make our intuitions more precise, we proceed in two steps.

First, let us consider an ensemble p(𝒙)p(\mbox{\boldmath\(x\)}) of stimuli 𝒙x with a ‘power’ constraint |𝒙|2<c|\mbox{\boldmath\(x\)}|^{2}<c. Intuitively, the power constraint means that the maximal light intensity across the whole screen is limited. In this case, the stimulus that is most likely to generate a spike under the linear receptive field model (11.10), is the one which is aligned with the receptive field

𝒙opt𝒌;\mbox{\boldmath\(x\)}_{\rm opt}\propto\mbox{\boldmath\(k\)}\,; (11.7)

cf. Exercises. Thus the receptive field vector 𝒌k can be interpreted as the optimal stimulus to cause a spike.

Second, let us consider an ensemble of stimuli 𝒙x with a radially-symmetric distribution, where the probability of a possibly multidimensional 𝒙x is equal to the probability of observing its norm |𝒙||\mbox{\boldmath\(x\)}| : p(𝒙)=pc(|𝒙|)p(\mbox{\boldmath\(x\)})=p_{c}(|\mbox{\boldmath\(x\)}|). Examples include the standard Gaussian distribution, or the uniform distribution with power constraint p(𝒙)=p0p(\mbox{\boldmath\(x\)})=p_{0} for |𝒙|2<c|\mbox{\boldmath\(x\)}|^{2}<c and zero otherwise. We assume that spikes are generated with the LNP model of Eq. (11.5). An important result is that the experimental reverse correlation technique yields an unbiased estimator of the filter 𝒌k, i.e.,

𝒙RevCorr=𝒌.\langle\mbox{\boldmath\(x\)}_{\rm RevCorr}\rangle=\mbox{\boldmath\(k\)}\,. (11.8)

The proof (84; 478) follows from the fact that each arbitrary input vector 𝒙t\mbox{\boldmath\(x\)}_{t} can be separated into a component parallel to 𝒌k and one orthogonal to it. Since we are free to choose the scale of the filter 𝒌k we can impose |𝒌|=1|\mbox{\boldmath\(k\)}|=1 and write

𝒙t=(𝒌𝒙t)𝒌+(𝐞𝒙t)𝐞\mbox{\boldmath\(x\)}_{t}=(\mbox{\boldmath\(k\)}\cdot\mbox{\boldmath\(x\)}_{t}% )\,\mbox{\boldmath\(k\)}+({\bf e}\cdot\mbox{\boldmath\(x\)}_{t})\,{\bf e} (11.9)

where 𝐞{\bf e} is a unit vector in the subspace orthogonal to 𝒌k. For firing only the component parallel to 𝒌k matters. The symmetry of the distribution p(𝒙)p(\mbox{\boldmath\(x\)}) guarantees that spike-triggered averaging is insensitive to the component orthogonal to 𝒌k; cf. Exercises.

In summary, reverse correlations are an experimental technique to determine the receptive field properties of a sensory neuron under an LNP model. The success of the reverse correlation technique as an experimental approach is intimately linked to its interpretability in terms of the LNP model.

Reverse correlations in the LNP model can also be analyzed in a statistical framework. To keep the arguments simple, we focus on the linear case and set

f(𝒌𝒙t)=ρ0+𝒌𝒙t.f(\mbox{\boldmath\(k\)}\cdot\mbox{\boldmath\(x\)}_{t})=\rho_{0}+\mbox{% \boldmath\(k\)}\cdot\mbox{\boldmath\(x\)}_{t}\,. (11.10)

The parameters minimizing the squared error between model firing rate ρ0+𝒌\rho_{0}+\mbox{\boldmath\(k\)}\cdot and observed firing rate ntn_{t} are then

𝒌opt=(XTX)-1(tnt𝒙t)/   d   t.\mbox{\boldmath\(k\)}_{opt}=(X^{T}X)^{-1}\left(\sum_{t}n_{t}\mbox{\boldmath\(x% \)}_{t}\right)/{\text{d}}t\,. (11.11)

We note that for short observation intervals Δ\Delta, the spike count ntn_{t} is either zero or one. Therefore the term in parenthesis on the right-hand side is proportional to the classical spike-triggered average; cf. Eq. (11.6). The factor in front of the parenthesis, XTXX^{T}X, is a scaled estimate of the covariance of the inputs 𝒙x. For stimuli consisting of uncorrelated white noise or light dots at random positions, the covariance structure is particularly simple.

See (383) for further connections between reverse correlation and likelihood-based estimates of the parameters in the LNP model.

Fig. 11.7: Example predictions of retinal ganglion ON-cell activity using the generalized linear encoding model with and without spike history terms. A: Recorded responses to repeated full-field light stimulus (top) of true ON-cell (‘RGC’), simulated LNP model (no spike history terms; ‘LNP’), and generalized linear model including spike-history terms (‘GLM’). Each row corresponds to the response during a single stimulus presentation. B. Magnified sections of rasters, with rows sorted in order of first spike time within the window in order to show spike timing details. Note that the predictions of the model including spike history terms are in each case more accurate than those of the Poisson (LNP) model. (PSTH variance accounted for: 91%, compared to 39% for the LNP model); C. Time-dependent firing rate plotted as a PSTH. D. Variance of the time-dependent firing rate. All data shown here are cross-validated ‘test’ data (that is, the estimated model parameters were in each case computed based on a non-overlapping ‘training’ data set not shown here). From (380) based on data from (528).

11.2.2 Multiple Neurons

Using multi-electrode arrays, Pillow et al. (2008) recorded from multiple ganglion cells in the retina provided with spatiotemporal white noise stimuli. This stimulation reaches the ganglion cells after being transducted by photoreceptors and inter-neurons of the retina. It is assumed that the effect of light stimulation can be taken into account by a linear filter of the spatio-temporally structured stimulus. An SRM-like model of the membrane potential of a neuron ii surrounded by nn other neurons is

ui(t)=fηi(t-ti(f))+𝐤i𝐱(t)+jifϵij(t-tj(f))+urest.u_{i}(t)=\sum_{f}\eta_{i}(t-t_{i}^{(f)})+\mathbf{k}_{i}\cdot\mathbf{x}(t)+\sum% _{j\neq i}\sum_{f}\epsilon_{ij}(t-t_{j}^{(f)})+u_{\rm rest}. (11.12)

The light stimulus 𝐱(t)\mathbf{x}(t) filtered by the receptive field of neuron ii, 𝐤i\mathbf{k}_{i}, replaces the artificially injected external current in Eq. (11.1). The spike afterpotential ηi(t)\eta_{i}(t) affects the membrane potential as a function of the neuron’s own spikes. Spikes from a neighboring neuron jj modify the membrane potential of neuron ii according to the coupling function ϵij(t)\epsilon_{ij}(t).

The extracellular electrodes used by Pillow et al. (399) did not probe the membrane potential. Nonetheless, by comparing the spike times with the conditional firing intensity ρ(t|{S})=1τ0exp(u(t)-θ0)\rho(t|\{S\})={1\over\tau_{0}}\,\exp\left(u(t)-\theta_{0}\right) we can maximize the likelihood of observing the set of spike trains {S}\{S\} (Ch. 10). This way we can identify the spatio-temporal receptive field 𝐤i\mathbf{k}_{i}, the spike afterpotential η(t)\eta(t) and the coupling functions ϵij(t)\epsilon_{ij}(t).

The fitted functions 𝐤i\mathbf{k}_{i} showed two types of receptive fields (399). The ON cells were sensitive to recent increase in luminosity while the OFF cells were sensitive to recent decrease. The coupling functions also reflect the two different neuron types. The coupling from ON cells to ON cells is excitatory and the coupling from ON cells to OFF cells is inhibitory, and conversely for couplings from OFF cells.

How accurate are the predictions of the multi-neuron model? Figure 11.8 describes the prediction performance. The spike trains and PSTHs of the real and modeled neurons are similar. The spike train likelihood reaches 2 bits per spike and the PSTH is predicted with 80-93% accuracy. Overall, the coupled model appears as a valid description of neurons embedded in a network.

Pillow et al. (399) also asked about the relevance of coupling between neurons. Are the coupling functions an essential part of the model or can the activity be accurately predicted without them? Optimizing the model with and without the coupling function independently, they found that the prediction of PSTH variance was unaffected. The spike prediction performance, however, showed a consistent improvement for the coupled model. (See (536) for further analysis using a model incorporating unobserved common noise effects.) Inter-neuron coupling played a greater role in decoding as we will see in the next section.

Fig. 11.8: Spike-train prediction of a retinal ganglion cell within its network. A, Raster of responses of a retinal ganglion cell (RGC; top) to 25 repetitions 1-s stimulus, and responses of the fully coupled model (Full model; bottom) to the same stimulus. B. PSTH of the RGC (full black line) and the fully coupled model (dashed black line). C. PSTH prediction for the fully coupled model of different cells plotted against the the PSTH prediction of a model fitted without inter-neuron coupling. D. Log-likelihood (Eq. (10.41)) for the fully coupled model of different cells plotted against the the log-likelihood of a model fitted without inter-neuron coupling. Modified from Pillow et al. (399).