Neural Encoding Models Maneesh Sahani Gatsby Computational - PowerPoint PPT Presentation

Spikes, or rate? Most neurons communicate using action potentials — statistically described by a point process: � � spike ∈ [ t , t + dt ) = λ ( t | H ( t ) , stimulus , network activity ) dt P To fully model the response we need to identify λ . In general this depends on spike history H ( t ) and network activity. Three options: ◮ Ignore the history dependence, take network activity as source of “noise” (i.e. assume firing is inhomogeneous Poisson or Cox process, conditioned on the stimulus).

Spikes, or rate? Most neurons communicate using action potentials — statistically described by a point process: � � spike ∈ [ t , t + dt ) = λ ( t | H ( t ) , stimulus , network activity ) dt P To fully model the response we need to identify λ . In general this depends on spike history H ( t ) and network activity. Three options: ◮ Ignore the history dependence, take network activity as source of “noise” (i.e. assume firing is inhomogeneous Poisson or Cox process, conditioned on the stimulus). ◮ Average multiple trials to estimate the mean intensity (or PSTH) � 1 λ ( t , stimulus ) = λ ( t | H n ( t ) , stimulus , network n ) , lim N N →∞ n and try to fit this.

Spikes, or rate? Most neurons communicate using action potentials — statistically described by a point process: � � spike ∈ [ t , t + dt ) = λ ( t | H ( t ) , stimulus , network activity ) dt P To fully model the response we need to identify λ . In general this depends on spike history H ( t ) and network activity. Three options: ◮ Ignore the history dependence, take network activity as source of “noise” (i.e. assume firing is inhomogeneous Poisson or Cox process, conditioned on the stimulus). ◮ Average multiple trials to estimate the mean intensity (or PSTH) � 1 λ ( t , stimulus ) = λ ( t | H n ( t ) , stimulus , network n ) , lim N N →∞ n and try to fit this. ◮ Attempt to capture history and network effects in simple models.

Spike-triggered average Decoding: mean of P ( s | r = 1 )

Spike-triggered average Decoding: mean of P ( s | r = 1 ) Encoding: predictive filter

Linear regression s 1 s 2 s 3 . . . s T s T + 1 . . . � T r ( t ) = s ( t − τ ) w ( τ ) d τ 0

Linear regression s 1 s 2 s 3 . . . s T s T + 1 . . . s 1 s 2 s 3 . . . s T � T � �� � r ( t ) = s ( t − τ ) w ( τ ) d τ 0 w t r T s 1 s 2 s 3 . . . s T + 1 . . . × = w 3 w 2 w 1

Linear regression s 1 s 2 s 3 . . . s T s T + 1 . . . s 1 s 2 s 3 . . . s T � T � �� � s 1 s 2 s 3 . . . s T s T r ( t ) = s ( t − τ ) w ( τ ) d τ � �� � 0 . . . s 1 s 2 s 3 s T + 1 w t r T s 2 s 3 s 4 . . . s T + 1 . r T + 1 . . × = w 3 . . . . . w 2 . w 1

Linear regression s 1 s 2 s 3 . . . s T s T + 1 . . . s 1 s 2 s 3 . . . s T � T � �� � s 1 s 2 s 3 . . . s T s T r ( t ) = s ( t − τ ) w ( τ ) d τ � �� � 0 . . . s 1 s 2 s 3 s T + 1 w t r T s 2 s 3 s 4 . . . s T + 1 . r T + 1 . . × = w 3 . . . . . w 2 . w 1 SW = R

Linear regression s 1 s 2 s 3 . . . s T s T + 1 . . . s 1 s 2 s 3 . . . s T � T � �� � s 1 s 2 s 3 . . . s T s T r ( t ) = s ( t − τ ) w ( τ ) d τ � �� � 0 . . . s 1 s 2 s 3 s T + 1 w t r T s 2 s 3 s 4 . . . s T + 1 . r T + 1 . . × = w 3 . . . . . w 2 . w 1 SW = R − 1 ( S T R ) W ( ω ) = S ( ω ) ∗ R ( ω ) W = ( S T S ) � �� � � �� � | S ( ω ) | 2 Σ SS STA

Linear models So the (whitened) spike-triggered average gives the minimum-squared-error linear model. Issues:

Linear models So the (whitened) spike-triggered average gives the minimum-squared-error linear model. Issues: ◮ overfitting and regularisation

Linear models So the (whitened) spike-triggered average gives the minimum-squared-error linear model. Issues: ◮ overfitting and regularisation ◮ standard methods for regression

Linear models So the (whitened) spike-triggered average gives the minimum-squared-error linear model. Issues: ◮ overfitting and regularisation ◮ standard methods for regression ◮ negative predicted rates

Linear models So the (whitened) spike-triggered average gives the minimum-squared-error linear model. Issues: ◮ overfitting and regularisation ◮ standard methods for regression ◮ negative predicted rates ◮ can model deviations from background

Linear models So the (whitened) spike-triggered average gives the minimum-squared-error linear model. Issues: ◮ overfitting and regularisation ◮ standard methods for regression ◮ negative predicted rates ◮ can model deviations from background ◮ real neurons aren’t linear

Linear models So the (whitened) spike-triggered average gives the minimum-squared-error linear model. Issues: ◮ overfitting and regularisation ◮ standard methods for regression ◮ negative predicted rates ◮ can model deviations from background ◮ real neurons aren’t linear ◮ models are still used extensively

Linear models So the (whitened) spike-triggered average gives the minimum-squared-error linear model. Issues: ◮ overfitting and regularisation ◮ standard methods for regression ◮ negative predicted rates ◮ can model deviations from background ◮ real neurons aren’t linear ◮ models are still used extensively ◮ interpretable suggestions of underlying sensitivity (but see later)

Linear models So the (whitened) spike-triggered average gives the minimum-squared-error linear model. Issues: ◮ overfitting and regularisation ◮ standard methods for regression ◮ negative predicted rates ◮ can model deviations from background ◮ real neurons aren’t linear ◮ models are still used extensively ◮ interpretable suggestions of underlying sensitivity (but see later) ◮ may provide unbiased estimates of cascade filters (see later)

How good are linear predictions? We would like an absolute measure of model performance. Two things make this difficult:

How good are linear predictions? We would like an absolute measure of model performance. Two things make this difficult: Measured responses can never be predicted perfectly, even in principle: ◮ The measurements themselves are noisy.

How good are linear predictions? We would like an absolute measure of model performance. Two things make this difficult: Measured responses can never be predicted perfectly, even in principle: ◮ The measurements themselves are noisy. Even if we can discount this, a model may predict poorly because either: ◮ It is the wrong model. ◮ The parameters are mis-estimated due to noise.

How good are linear predictions? We would like an absolute measure of model performance. Two things make this difficult: Measured responses can never be predicted perfectly, even in principle: ◮ The measurements themselves are noisy. Even if we can discount this, a model may predict poorly because either: ◮ It is the wrong model. ◮ The parameters are mis-estimated due to noise. Approaches: ◮ Compare I ( resp ; pred ) to I ( resp ; stim ) . ◮ mutual information estimators are biased ◮ Compare E ( resp − pred ) to E ( resp − psth ) where psth is gathered over a very large number of trials. ◮ may require impractical amounts of data to estimate the psth ◮ Compare the predictive power to the predicatable power (similar to ANOVA).

Estimating predictable power P signal P noise response = signal + noise � �� � r ( n )   � � 1 P ( r ( n ) ) = P signal + P noise   � N P ( r ( n ) ) − P ( r ( n ) )   P signal = N − 1 ⇒ P ( r ( n ) ) = P signal + 1     N P noise � P noise = P ( r ( n ) ) − � P signal

Testing a model For a perfect prediction � � P ( trial ) − P ( residual ) = P ( signal )

Testing a model For a perfect prediction � � P ( trial ) − P ( residual ) = P ( signal ) Thus, we can judge the performance of a model by the normalized predictive power P ( trial ) − P ( residual ) � P ( signal )

Testing a model For a perfect prediction � � P ( trial ) − P ( residual ) = P ( signal ) Thus, we can judge the performance of a model by the normalized predictive power P ( trial ) − P ( residual ) � P ( signal ) Similar to coefficient of determination ( r 2 ), but the denominator is the predictable variance.

Predictive performance Training Error Cross−Validation Error 2.5 1 normalised Bayes predictive power 0.5 2 0 1.5 −0.5 1 −1 0.5 −1.5 0 −2 0 0.5 1 1.5 2 2.5 −2 −1 0 1 normalised STA predictive power

Extrapolating the model performance 3 normalized linearly predictive power 2.5 2 1.5 1 0.5 0 −0.5 0 50 100 150 normalized noise power

Extrapolating the model performance 3 normalized linearly predictive power 2.5 2 1.5 1 0.5 0 −0.5 0 50 100 150 normalized noise power

Extrapolating the model performance 3 normalized linearly predictive power 2.5 2 1.5 1 0.5 0 −0.5 0 50 100 150 normalized noise power

Extrapolating the model performance 3 normalized linearly predictive power 2.5 2 1.5 1 0.5 0 −0.5 0 50 100 150 normalized noise power

Jacknife bias correction Estimate bias by extrapolation in data size: T jn = N T − ( N − 1 ) T loo where T is the training error on all data and T loo is the average training error on all sets of N − 1 data. For a linear model we can find this in closed form: � � � ( r i − s i w ML ) 2 T jn = 1 1 − s i ( S T S ) − 1 s T N i i

Jackknifed estimates 3 2.5 Normalized linearly predictive power 2 1.5 1 0.5 0 −0.5 0 50 100 150 Normalized noise power

Extrapolated linearity 2.5 2 1.5 1 1 0.5 0.8 0 Normalized linearly predictive power −0.5 0 50 100 150 0.6 0.4 0.2 0 −0.2 −5 0 5 10 15 20 25 30 Normalized noise power [extrapolated range: (0.19,0.39); mean Jackknife estimate: 0.29]

Simulated (almost) linear data 3 2.5 2 1.5 1.5 1.4 1 0.5 Normalized linearly predictive power 1.3 0 0 50 100 150 1.2 1.1 1 0.9 0.8 0.7 0.6 0.5 −5 0 5 10 15 20 25 30 Normalized noise power [extrapolated range: (0.95,0.97); mean Jackknife estimate: 0.97]

Beyond linearity

Beyond linearity Linear models often fail to predict well. Alternatives? ◮ Wiener/Volterra functional expansions ◮ M-series ◮ Linearised estimation ◮ Kernel formulations ◮ LN (Wiener) cascades ◮ Spike-trigger covariance (STC) methods ◮ “Maximimally informative” dimensions (MID) ⇔ ML nonparametric LNP models ◮ ML Parametric GLM models ◮ NL (Hammerstein) cascades ◮ Multilinear formulations

The Volterra functional expansion A polynomial-like expansion for functionals (or operators). Let y ( t ) = F [ x ( t )] . Then: � �� y ( t ) ≈ k ( 0 ) + dx d τ k ( 1 ) ( τ ) x ( t − τ ) + dx d τ 1 d τ 2 k ( 2 ) ( τ 1 , τ 2 ) x ( t − τ 1 ) x ( t − τ 2 ) ��� dx d τ 1 d τ 2 d τ 3 k ( 3 ) ( τ 1 , τ 2 , τ 3 ) x ( t − τ 1 ) x ( t − τ 2 ) x ( t − τ 3 ) + . . . + or (in discretised time) � � � y t = K ( 0 ) + K ( 1 ) K ( 2 ) K ( 3 ) x t − i + x t − i x t − j + ijk x t − i x t − j x t − k + . . . i ij i ij ijk For finite expansion, the kernels k ( 0 ) , k ( 1 ) ( · ) , k ( 2 ) ( · , · ) , k ( 3 ) ( · , · , · ) , . . . are not straightforwardly related to the functional F . Indeed, values of lower-order kernels change as the maximum order of the expansion is increased. Estimation: model is linear in kernels, so can be estimated just like a linear (first-order) model with expanded “input”. ◮ Kernel trick: polynomial kernel K ( x 1 , x 2 ) = ( 1 + x 1 x 2 ) n . ◮ M-series.

Wiener Expansion The Wiener expansion gives functionals of different orders that are orthogonal for white noise input x ( t ) . G 0 [ x ( t ); h ( 0 ) ] = h ( 0 ) � G 1 [ x ( t ); h ( 1 ) ] = dx d τ h ( 1 ) ( τ ) x ( t − τ ) �� � G 2 [ x ( t ); h ( 2 ) ] = dx d τ 1 d τ 2 h ( 2 ) ( τ 1 , τ 2 ) x ( t − τ 1 ) x ( t − τ 2 ) − P dx d τ 1 h ( 2 ) ( τ 1 , τ 1 ) ��� G 3 [ x ( t ); h ( 3 ) ] = dx d τ 1 d τ 2 d τ 3 h ( 3 ) ( τ 1 , τ 2 , τ 3 ) x ( t − τ 1 ) x ( t − τ 2 ) x ( t − τ 3 ) �� dx d τ 1 d τ 2 h ( 3 ) ( τ 1 , τ 2 , τ 2 ) x ( t − τ 1 ) − 3 P Easy to verify that E [ G i [ x ( t )] G j [ x ( t )]] = 0 for i � = j . Thus, these kernels can be estimated independently. But, they depend on the stimulus.

Cascade models The LNP (Wiener) cascade n k ◮ Rectification addresses negative firing rates. ◮ Loose biophysical correspondance.

LNP estimation – the Spike-triggered ensemble

Single linear filter n k ◮ STA is unbiased estimate of filter for spherical input distribution. (Bussgang’s theorem) ◮ Elliptically-distributed data can be whitened ⇒ linear regression weights are unbiased. ◮ Linear weights are not necessarily maximum-likelihood (or otherwise optimal), even for spherical/elliptical stimulus distributions. ◮ Linear weights may be biased for general stimuli (binary/uniform or natural).

Multiple filters Distribution changes along relevant directions (and, usually, along all linear combinations of relevant directions). Proxies to measure change in distribution: ◮ mean: STA (can only reveal a single direction) ◮ variance: STC ◮ binned (or kernel) KL divergence: MID “maximally informative directions” (equivalent to ML in LNP model with binned nonlinearity)

STC Project out STA: X T � � X T diag ( Y ) � � X X � X = X − ( X k sta ) k T sta ; C prior = N ; C spike = N spike Choose directions with greatest change in variance: v T ( C prior − C spike ) v k- argmax � v � = 1 ⇒ find eigenvectors of ( C prior − C spike ) with large (absolute) eigvals.

STC Reconstruct nonlinearity (may assume separability)

Biases STC (obviously) requires that the nonlinearity alter variance. If so, subspace is unbiased provided distribution is ◮ radially (elliptically) symmetric ◮ AND independent ⇒ Gaussian. May be possible to correct for non-Gaussian stimulus by transformation, subsampling or weighting (latter two at cost of variance).

More LNP methods ◮ Non-parametric non-linearities: “Maximally informative dimensions” (MID) ⇔ “non-parametric” maximum likelihood. ◮ Intuitively, extends the variance difference idea to arbitrary differences between marginal and spike-conditioned stimulus distributions. k MID = argmax KL [ P ( k · x ) � P ( k · x | spike )] k ◮ Measuring KL requires binning or smoothing—turns out to be equivalent to fitting a non-parametric nonlinearity by binning or smoothing. ◮ Difficult to use for high-dimensional LNP models (but ML viewpoint suggests separable or “cylindrical” basis functions). ◮ Parametric non-linearities: the “generalised linear model” (GLM).

Generalised linear models LN models with specified nonlinearities and exponential-family noise. In general (for monotonic g ): y ∼ ExpFamily [ µ ( x )]; g ( µ ) = β x For our purposes easier to write y ∼ ExpFamily [ f ( β x )] (Continuous time) point process likelihood with GLM-like dependence of λ on covariates is approached in limit of bins → 0 by either Poisson or Bernoulli GLM. Mark Berman and T. Rolf Turner (1992) Approximating Point Process Likelihoods with GLIM Journal of the Royal Statistical Society. Series C (Applied Statistics), 41(1):31-38.

Generalised linear models Poisson distribution ⇒ f = exp () is canonical ( natural params = β x ). Canonical link functions give concave likelihoods ⇒ unique maxima. Generalises (for Poisson) to any f which is convex and log-concave: log-likelihood = c − f ( β x ) + y log f ( β x ) Includes: ◮ threshold-linear ◮ threshold-polynomial f ( z ) = [ z 3 ] + ◮ “soft-threshold” f ( z ) = α − 1 log ( 1 + e α z ) . f ( z ) f ( z ) = log ( 1 + e z ) f ( z ) = 1 3 log ( 1 + e 3 z ) f ( z ) = [ z ] + z

Generalised linear models ML parameters found by ◮ gradient ascent ◮ IRLS Regularisation by L 2 (quadratic) or L 1 (absolute value – sparse) penalties (MAP with Gaussian/Laplacian priors) preserves concavity.

Linear-Nonlinear-Poisson (GLM) Poisson point stimulus filter nonlinearity spiking stimulus ( t ) k

GLM with history-dependence (Truccolo et al 04) Poisson exponential stimulus filter nonlinearity spiking + stimulus ( t ) k post-spike filter h conditional intensity (spike rate) • rate is a product of stim- and spike-history dependent terms • output no longer a Poisson process • also known as “soft-threshold” Integrate-and-Fire model

GLM with history-dependence Poisson exponential stimulus filter nonlinearity spiking + stimulus ! ( t ) k post-spike filter h traditional IF “soft-threshold” IF spike rate “hard threshold” filter output filter output • “soft-threshold” approximation to Integrate-and-Fire model

GLM dynamic behaviors stimulus x(t) post-spike waveform regular spiking stim-induced 0 0 spike-history induced 0 50 100 0 100 200 300 400 500 time after spike time (ms)

GLM dynamic behaviors stimulus x(t) post-spike waveform regular spiking stim-filter output 0 0 spike-history filter output irregular spiking 0 0 0 10 20 0 100 200 300 400 500 time after spike time (ms)

GLM dynamic behaviors stimulus x(t) post-spike waveform bursting 0 0 adaptation 0 0 -10 0 20 40 0 100 200 300 400 500 time after spike time (ms)

Generalized Linear Model (GLM) exponential probabilistic stimulus filter nonlinearity spiking + stimulus post-spike filter

multi-neuron GLM exponential probabilistic stimulus filter nonlinearity spiking + post-spike filter neuron 1 stimulus + neuron 2

multi-neuron GLM exponential probabilistic stimulus filter nonlinearity spiking + post-spike filter neuron 1 stimulus coupling filters + neuron 2

GLM equivalent diagram: ... time t conditional intensity (spike rate)

Non-LN models? The idea of responses depending on one or a few linear stimulus projections has been dominant, but cannot capture all non-linearities. ◮ Contrast sensitivity might require normalisation by � s � . ◮ Linear weighting may depend on units of stimulus measurement: amplitude? energy? logarithms? thresholds? (NL models – Hammerstein cascades) ◮ Neurons, particularly in the auditory system are known to be sensitive to combinations of inputs: forward suppression; spectral patterns (Young); time-frequency interactions (Sadogopan and Wang). ◮ Experiments with realistic stimuli reveal nonlinear sensivity to parts/whole (Bar-Yosef and Nelken). Many of these questions can be tackled using a multilinear (cartesian tensor) framework.

Input nonlinearities The basic linear model (for sounds): � w tf � r ( i ) = s ( i − j , k ) , jk ���� ���� � �� � jk predicted rate STRF weights stimulus power

Input nonlinearities The basic linear model (for sounds): � w tf � r ( i ) = s ( i − j , k ) , jk ���� ���� � �� � jk predicted rate STRF weights stimulus power How to measure s ? (pressure, intensity, dB, thresholded, . . . )

Input nonlinearities The basic linear model (for sounds): � w tf � r ( i ) = s ( i − j , k ) , jk ���� ���� � �� � jk predicted rate STRF weights stimulus power How to measure s ? (pressure, intensity, dB, thresholded, . . . ) We can learn an optimal representation g ( . ) : � w tf ˆ r ( i ) = jk g ( s ( i − j , k )) . jk

Input nonlinearities The basic linear model (for sounds): � w tf � r ( i ) = s ( i − j , k ) , jk ���� ���� � �� � jk predicted rate STRF weights stimulus power How to measure s ? (pressure, intensity, dB, thresholded, . . . ) We can learn an optimal representation g ( . ) : � w tf ˆ r ( i ) = jk g ( s ( i − j , k )) . jk Define: basis functions { g l } such that g ( s ) = � l w l l g l ( s ) and stimulus array M ijkl = g l ( s ( i − j , k )) . Now the model is � w tf jk w l ˆ r ( i ) = l M ijkl j

Input nonlinearities The basic linear model (for sounds): � w tf � r ( i ) = s ( i − j , k ) , jk ���� ���� � �� � jk predicted rate STRF weights stimulus power How to measure s ? (pressure, intensity, dB, thresholded, . . . ) We can learn an optimal representation g ( . ) : � w tf ˆ r ( i ) = jk g ( s ( i − j , k )) . jk Define: basis functions { g l } such that g ( s ) = � l w l l g l ( s ) and stimulus array M ijkl = g l ( s ( i − j , k )) . Now the model is � r = ( w tf ⊗ w l ) • M . w tf jk w l r ( i ) = ˆ � l M ijkl or j

Multilinear models Multilinear forms are straightforward to optimise by alternating least squares. Cost function: � � 2 � r − ( w tf ⊗ w l ) • M � � E = � Minimise iteratively, defining matrices B = w l • M A = w tf • M and and updating w tf = ( B T B ) − 1 B T r w l = ( A T A ) − 1 A T r . and Each linear regression step can be regularised by evidence optimisation (suboptimal), with uncertainty propagated approximately using variational methods.

Some input non-linearities w l 0 25 40 55 70 l (dB−SPL)

Parameter grouping Separable models: (time) ⊗ (frequency). The input nonlinearity model is separable in another sense: (time, frequency) ⊗ (sound level). frequency frequency weight time intensity intensity time Other separations: r = ( w tl ⊗ w f ) • M , ◮ (time, sound level) ⊗ (frequency): � r = ( w fl ⊗ w t ) • M , ◮ (frequency, sound level) ⊗ (time): � r = ( w l ⊗ w f ⊗ w l ) • M . ◮ (time) ⊗ (frequency) ⊗ (sound level): �

Some examples 32 16 f (kHz) w l 8 (time, frequency) ⊗ (sound level): 0 4 0 2 180 120 60 0 25 40 55 70 t (ms) l (dB−SPL) 70 l (dB−SPL) 55 w f (time, sound level) ⊗ (frequency): 40 0 0 25 180 120 60 0 25 50 100 t (ms) f (kHz) 70 l (dB−SPL) 0 55 w t (frequency, sound level) ⊗ (time): 40 0 25 2 4 8 16 32 180 120 60 0 f (kHz) t (ms)

Variable (combination-dependent) input gain ◮ Sensitivities to different points in sensory space are not independent.

Variable (combination-dependent) input gain ◮ Sensitivities to different points in sensory space are not independent. ◮ Rather, the sensitivity at one point depends on other elements of the stimulus that create a local sensory context.

Variable (combination-dependent) input gain ◮ Sensitivities to different points in sensory space are not independent. ◮ Rather, the sensitivity at one point depends on other elements of the stimulus that create a local sensory context. ◮ This context adjusts the input gain of the cell from moment to moment, dynamically refining the shape of the weighted receptive field.

A context-sensitive model s ( i , k ) r ( i )

A context-sensitive model � J � K w tf ˆ r ( i ) = c + j + 1 , k s ( i − j , k ) j = 0 k = 1

A context-sensitive model � � � J � K � M � N w tf w τφ ˆ r ( i ) = c + j + 1 , k s ( i − j , k ) 1 + m + 1 , n + N + 1 s ( i − j − m , k + n ) j = 0 k = 1 m = 0 n = − N