Neuronal Dynamics (11)
Encoding and Decoding with Stochastic Neuron models
The online version of this chapter:
Chapter 11 Encoding and Decoding with Stochastic Neuron models https://neuronaldynamics.epfl.ch/online/Ch11.html
Encoding Models for Intracellular Recordings
We focus on GLMs with escape noise, also called soft-threshold integrate-and-fire models.
Predicting Membrane Potential
The SRM model \[ u(t)=\sum_{f}^{} \eta(t-t^{(f)})+\int_{0}^{\infty} \kappa(s)I^{ext}(t-s) \mathrm{d}s +u_{rest}. \tag{11.1} \]
For both the main type of excitatory neurons and the main type on inhibitory neurons, the membrane filter \(\kappa(t)\) is well described by a single exponential. Different cell types have different amplitudes and time constants. The inhibitory neurons are typically faster, with a smaller time constant than the excitatory neurons, suggesting we could discriminate between excitatory and inhibitory neurons in terms of the shape of \(\kappa(t)\). When we take into account the spike-afterpotential, discrimination of cell types is much improved. The shape of \(\eta(t)\) in inhibitory cells is very different than that in excitatory ones.
While the spike-afterpotential is a monotonically decreasing function in the excitatory cells, in the inhibitory cells the function \(\eta(t)\) is better fitted by two exponentials of opposite polarity. Spike afterpotential of inhibitory neurons has an oscillatory component.
If we set \(\kappa(s)=(1/C)\exp (-s/\tau_m)\), we can take the derivative of (11.1) and write it in the form of \[ C\frac{\mathrm{d}u(t)}{\mathrm{d}t}=-\frac{1}{R}(u-u_{rest})+\sum_{f}^{} \tilde{\eta}(t-t^{(f)})+I^{ext}(t) \tag{11.2} \] where \(\tilde{\eta}\) is the time course of the net current triggered after a spike.
After every spike the total current that can charge the membrane capacitance is \[ C\frac{\mathrm{d}u}{\mathrm{d}t} \propto \eta_{C}(t-\hat{t})(u-E_{rev}) \tag{11.3} \] where \(\eta_{C}\) is the spike triggered change in conductance and \(E_{rev}\) its reversal potential. The prediction performance based on conductance change instead of spike after-effects in terms of potential is not significantly improved.
Predicting Spikes
Assuming a moving threshold that can undergo a stereotypical change at every spike \(\theta(t)=\theta_0+\sum_{f}^{} \theta_1(t-t^{(f)})\) we can model the conditional firing intensity as (c.f., 9.27) \[ p(t|S)=\frac{1}{\tau_0}\exp \left[ \beta\left(u(t)-\theta_0-\sum_{t^{(f)}\in S}^{} \theta_1(t-t^{(f)}) \right)\right] \tag{11.4} \]
Since the parameters regulating \(u(t)\) were optimized using the subthreshold membrane potential in 11.1.1, the only free parameters left are those of the threshold, that is, \(\theta_0\),\(\beta\) and the function \(\theta_1(t)\). Once the function \(\theta_1\) is expanded in a linear combination of basis functions, maximizing the likelihood (10.40) can be done through a convex gradient descent.
Dynamic threshold is not necessary for some neurons. The threshold in those cells is constant in time. However, the excitatory cells have a strongly moving threshold which is characterized by at least two decay time constants. Inactivation of sodium channels is a likely candidate for the biophysical causes of a moving threshold.
GLMs predict more than 80 percent of the 'predictable' spikes. Some cells were predicted better than others such that the \(M\) reached up to 95%. Other optimization methods but with similar models could improve the spike timing prediction of inhibitory neurons, reaching up to \(M\)=100% for some cells. However, optimizing a GLM model with refractory effects but no adaptation reduces the prediction performance by 20-30%, for both the excitatory and inhibitory cortical cells. A single spike has a measurable effect more than 10 seconds after the action potential has occurred. Thus, adaptation is not characterized by a single time scale and shows up as a power-law decay in both spike-triggered current and threshold.
Encoding Models in Systems Neuronscience
Receptive fields and Linear-Nonlinear Poisson Model
For a two-dimensional image, we label all pixels with a single index \(k\). A full image corresponds to a vector \(\mathbf{x}=(x_1,\cdots ,x_{K})\) while a single spot of light corresponds to a vector with all components equal to zero except one.
The spatial receptive field of a neuron is a vector \(\mathbf{k}\) of the same dimensionality as \(\mathbf{x}\). The response of the neuron to an arbitrary spatial stimulus \(\mathbf{x}\) depends on the total drive \(\mathbf{k}\cdot \mathbf{x}_t\), i.e., the similarity between the stimulus and the spatial filter.
More generally, the receptive field filter \(\mathbf{k}\) can be described not only by a spatial component, but also by a temporal component. The scalar product \(\mathbf{k}\cdot \mathbf{x}_t\) is a shorthand notation for integration over space as well as over time. Such a filter \(\mathbf{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 \[ \rho(t)=f(\mathbf{k}\cdot \mathbf{x}_t) \tag{11.5} \] 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 suprisingly similar.
Example: Detour on reverse correlation for receptive field estimation
Reverse correlation measurements are an experimental procedure based on spike-triggered averaging. Stimuli \(\mathbf{x}\) are drawn from some statistical ensemble and presented on after the other. Each time the neuron elicits a spike, the stimulus \(\mathbf{x}\) presented just before the firing is recorded. The reverse correlation filter is the mean of all inputs that have triggered a spike \[ \mathbf{x}_{RevCorr}=\langle \mathbf{x}\rangle _{spike}=\frac{\sum_{t}^{} n_t \mathbf{x}_t}{\sum_{t}^{} n_t},\tag{11.6} \] where \(n_t\) is the spike count in trial \(t\). The reverse correlation technique finds the typical stimulus that causes a spike.
Consider an ensemble \(p(\mathbf{x})\) of stimuli \(\mathbf{x}\) with a 'power' constraint \(\lvert \mathbf{x} \rvert ^{2}<c\). In this case, the stimulus that is most likely to generate a spike under the linear receptive field model is the one which is aligned with the receptive field \[ \mathbf{x}_{opt}\propto \mathbf{k} \tag{11.7} \] The receptive field vector \(\mathbf{k}\) can be interpreted as the optimal stimulus to cause a spike.
Then, consider an ensemble of stimuli \(\mathbf{x}\) with a radially-symmetric distribution, where the probability of a possibly multidimensional \(\mathbf{x}\) is equal to the probability of observing its norm \(\lvert \mathbf{x} \rvert \colon p(\mathbf{x})=p_c(\lvert \mathbf{x} \rvert )\). An important result is that the experimental reverse correlation technique yields an unbiased estimator of the filter \(\mathbf{k}\), i.e., \[ \langle \mathbf{x}_{RevCorr}\rangle =\mathbf{k}. \tag{11.8} \]
Multiple Neurons
An SRM-like model of the membrane potential of a neuron \(i\) surrounded by \(n\) other neurons is \[ u_i(t)=\sum_{f}^{} \eta_i(t-t_i^{(f)})+\mathbf{k}_i \cdot \mathbf{x}(t)+\sum_{j \neq i}^{} \sum_{f}^{} \varepsilon_{ij}(t-t_j^{(f)})+u_{rest}. \tag{11.12} \]
Decoding
In this section, we apply 'Bayes' rule to obtain the posterior probability of the stimulus, conditional on the observed response: \[ p(\mathbf{x}|D)\propto p(D|\mathbf{x})p(\mathbf{x}), \tag{11.13} \]
Maximum a posteriori decoding
The Maximum A Posteriori (MAP) estimate is the stimulus \(\mathbf{x}\) that is most probable given the observed spike response \(D\), i.e., the \(\mathbf{x}\) that maximizes \(p(\mathbf{x}|D)\).
The log-posterior, \[ \log p(\mathbf{x}|D)=\log p(D|\mathbf{x})+\log p(\mathbf{x})+c \tag{11.14} \] is concave as long as the stimulus log-prior \(\log p(\mathbf{x})\) is itself a concave function of \(\mathbf{x}\) (e.g. \(p\) is Gaussian). In this case, again, we may easily compute \(\hat{x}_{MAP}\) by numerically ascending the function \(\log p(\mathbf{x}|D)\).
The MAP estimate of the stimulus is, in general, a nonlinear function of the observed spiking data \(D\).
Example: Linear stimulus Reconstruction
We predict the stimulus \(\mathbf{x}_t\) by linear filtering of the observed spike times \(t^{1},t^{2},\cdots ,t^{F}<t\), \[ x(t)=x_0+\sum_{f}^{} k(t-t^{f}) \tag{11.15} \]
The aim is to find the shape of the filter \(k\), i.e., the optimal linear estimator (OLE) of the stimulus. It can be obtained using standard least-squares regression of the spiking data onto the stimulus \(\mathbf{x}\).
Assessing decoding uncertainty
In addition to providing a reliable estimate of the stimulus underlying a set of spike responses, computing the MAP estimate \(\hat{x}_{MAP}\) gives us easy access to the variance of the posterior distribution around \(\hat{x}_{MAP}\). It tells us something about which stimulus features are best encoded by the response \(D\).
For example, along stimulus axes where the posterior has small variance (i.e. the posterior declines rapidly as we move away from \(\hat{x}_{MAP}\)), we have relatively high certainty that the true \(\mathbf{x}\) is close to \(\hat{x}_{MAP}\). Conversely, we have relatively low certainty about any feature axis along which the posterior variance is large.
We measure the curvature at \(\hat{x}_{MAP}\) by computing the 'Hessian' matrix \(A\) of second-derivatives of the log-posterior, \[ A_{ij}=-\frac{\partial ^{2}}{\partial x_i \partial x_j}\log p(\mathbf{x}|D). \tag{11.16} \]
Moreover, the eigendecomposition of this matrix \(A\) tells us exactly which axes of stimulus space correspond to the 'best' and 'worst' encoded feature of the neural response: small eigenvalues of \(A\) correspond to directions of small curvature, where the observed data \(D\) poorly constrains the posterior distribution \(p(\mathbf{x}|D)\) (and therefore the posterior variance will be relatively large in this direction), while conversely large eigenvalues in \(A\) imply relatively precise knowledge of \(\mathbf{x}\), i.e., small posterior variance (for this reason the Hessian of the log-likelihood \(p(D|x)\) is referred to as the 'observed Fisher information matrix' in the statistic literature).
We can use this Hessian to construct a useful approximation to the posterior \(p(\mathbf{x}|D)\). The idea is simply to approximate this log-concave bump with a Gaussian function, where the parameters of the Gaussian are chosen to exactly match the peak and curvature of the true posterior. \[ p(\mathbf{x}|D)\thickapprox (2\pi)^{-d/2}\lvert A \rvert ^{1/2}\mathrm{e}^{-(\mathbf{x}-\hat{x}_{MAP})^{\mathsf{T}}A(\mathbf{x}-\hat{x}_{MAP})^{2}} , \tag{11.17} \] with \(d=\text{dim}(\mathbf{x})\). The approximate posterior entropy of variance of \(x_i\) is \(var(x_i|D)\thickapprox [A^{-1}]_{ii}\).