Data

Retinal tissue was obtained from adult (8 weeks old) male Long-Evans rat (Rattus norvegicus) and continuously perfused with Ames Solution (Sigma-Aldrich) and maintained at 32 °C. Ganglion cell spikes were recorded extracellularly from a multi-electrode array with 252 electrodes spaced 60 μm apart (custom fabrication by Innovative Micro Technologies, Santa Barbara, CA). Experiments were performed in accordance with institutional animal care standards. The microelectrode covered a total retinal area of ∼ 1 mm². For the rat this corresponds to 16-17 degrees of visual angle [25]. The spike sorting was performed with an in-house method based on [22].

Visual stimulus

The stimulus movie consisted of randomly moving dark discs (r = 100 μm) against a bright background (100% contrast, 2 ⋅ 10¹² photons/cm²/s). The discs followed mutually avoiding trajectories generated through an Ornstein-Uhlenbeck process of the form: (1) where v_i is the velocity of the disc i, Δt = 0.01 is the integration timestep, τ = 0.8 is the damping time constant, σ = 0.5 is the random force magnitude, dW is a zero-mean, unit-variance Gaussian random variable, and f_i is the hard-core central repulsive force between the discs and between each disc and the frame bounding box, with a decay of ∝ r⁻⁶, where r is the distance between the disc centers or between the disc and closest point of the bounding box. The resulting distribution of disc speeds peaked at v ≈ 0.6 μm/ms and had a width of about σ_v ≈ 0.4 μm/ms. The discs covered the recorded area uniformly to a very good approximation, with occupancy deviations at different encoding sites of ∼ 3%.

The movie was divided in segments of 1, 2, 4 and 10 discs, each 675 s long. Segments with increasing number of discs were presented sequentially and in total 3 segments of each type were shown, amounting to a total experiment time of 135 min. Each segment was regularly interspersed with 18 short (7.5 s) clips of repeated stimulus: in sum, 54 repeated clips were shown for each stimulus with different number of discs. The stimulus was convolved with a bank of 400 spatial symmetric gaussian filters (σ = 66.67 μm) placed in a regular square 20x20 grid with d = 53μm spacing to produce local luminance traces. The filter normalization ensures the resulting traces are bounded in (0,1). The width of the filters was selected in preliminary tests to optimize decoding performance; specifically, in preliminary tests we found that filter widths in range σ ∼ 50 − 100μm minimized the mean-squared-error of L1-regularized linear decoders trained for a subset of decoding sites. The movie stimulus was shown at a refresh rate of 80 Hz. The response spike trains were binned accordingly in bins of 12.5 ms, and time aligned to the stimulus. The spatio-temporal receptive fields of the retinal ganglion cells were obtained through reverse correlation to a flickering checkerboard stimulus. The checkerboard was constructed from squares of 130 μm that were randomly selected to be black or white at a rate of 40 Hz. Retinal spontaneous activity was recorded in full darkness (blackout condition) for 2.5 min.

Linear decoder

Let be a one-dimensional stimulus trace of length N time bins. In the linear decoding framework we assume that an estimate of the stimulus can be obtained from the neural response Σ as , where is a linear filter. In this formulation, the response of the retina is represented by the matrix , where C is the number of cells and ΔT the size in bins of the time window we associate with a single point in (for all analyses ΔT = 61 corresponding to a window stretching from -375 ms to 375 ms around the time bin of interest, i.e., the decoding is performed using “acausal” filters). The extra dimension is a column of ones to account for the bias term in the decoding. Thus, the decoding filter is structured as where is the filter corresponding to cell i and L₀ is the bias term. We learned the filters by minimizing the square error function with L1-regularization (2) To solve the minimization problem computationally we made use of the Lasso algorithm with the routines by Kim et al. [26]. Data was divided into training and testing sets (4.9⋅10⁴ training points, 2.3⋅10⁴ testing points). The filters were obtained from the training set and all measures of performance refer to the testing set. Regularization parameter λ was chosen through 2-fold cross-validation on the training set. The regularization term ensures the sparsity of the filters. Due to this sparsity some cells have negligible filter norms and therefore do not contribute to the decoding. This allows us to establish a hierarchy of cells by sorting them according to their filter norm . “Single-best cell” for every site refers to the cell with the largest norm. “Contributing cells” are the subset of cells with largest norm that jointly account for at least half of the total filter norm .

All the results in the paper use acausal filters. If the decoding filters are restricted to be causal, the decoding performance can be significantly decreased; for test sites where we explored this effect, the cross correlation between true and linearly decoded luminance trace could decrease from ≈ 0.8 to ≈ 0.6. Given that our stimulus is a stochastic process and retinal processing necessarily entails some processing delay, this is not surprising.

Kernelized nonlinear decoder

If instead of L1-regularization we enforce L2-regularization, the linear decoding filters can be obtained analytically through the normal equation Thus, an estimate of the stimulus for some new data is given by Since this expression only depends on products of spike trains, we can make use of the kernel trick and substitute the usual scalar product by some appropriate nonlinear function k of the spike trains. In this way, we can express our nonlinear decoding problem as (3) (4) (5) where is the ith row of matrix Σ. This is known as Kernel Ridge Regression [27, 28]. For our analyses we have used the Gaussian kernel (6) Before computing the kernel, it is customary to turn the spike trains into smooth traces for the sake of performance [29]. We convolved our spike trains with a Gaussian filter of 3 time bins width. The data was divided into training and testing sets (9.8⋅10³ training points, 2.3⋅10⁴ testing points). The parameters s and λ were obtained through joint 3-fold cross-validation on the training set. Since decoding at different sites is independent, s and λ were chosen separately at each site (likewise, L1 regularization strength for the linear decoder was also selected independently for each site). The performance of the nonlinear decoder depends on the set of cells considered. Contrary to the linear case where L1-regularization can effectively silence cells by setting their filters to zero, this nonlinear framework cannot ignore cells in a similar way. Therefore, including in the analysis non-informative cells can decrease the generalization performance of the decoder. To determine the best subset of cells for decoding we took advantage of the hierarchy of cells established by the linear L1-regularized decoding. We trained nonlinear decoders with progressively more cells (best cell, best two cells, etc.) and selected the subset of minimum decoding error on the training set (S7 Fig). Effectively, we jointly cross-validated the three parameters s, λ, and the subset size.

Deep neural network

We trained a deep neural network on the decoding task. The architecture of the network is as follows: there are 5460 inputs (activity of 91 cells and 60 time bins) and two or three fully connected hidden layers followed by a fully connected linear output layer of 400 cells (corresponding to the 20x20 grid), see S20 Fig. The hidden layers have each 150 units with tanh activation function.

The loss is L2 loss on the regression error and L1 and L2 regularization on the weights (Elastic Net type), more specifically: (7) where is the network output for input x_i and W^(k) is the weight matrix to layer k + 1, . The regularization improves generalization by making the network weights smaller and creating a sparse connection graph (increases robustness to training set variations). The available data of 126360 input-output pairs was split into 79560, 23400, and 23400 points for training, validation and testing respectively (same test set as for other methods). We trained for 2500 epochs (each epoch trains on all training points once). To avoid that early during training many weights become zero because of the regularization we set λ = 0 for the first 100 epochs. We performed model selection by a grid search through the following hyperparameters: regularization constant: λ ∈ {5 ⋅ 10⁻⁷, 7.5 ⋅ 10⁻⁷, 1 ⋅ 10⁻⁶, 2.5 ⋅ 10⁻⁶, 5 ⋅ 10⁻⁶}, number of hidden layers K ∈ {2, 3}, and optimization method ∈ { stochastic gradient decent (sgd) with learning rate 0.01 and momentum of 0.9, Adam optimizer [30] with learning rate 0.005 and ϵ = 0.0001}. The hyperparameter setting with the smallest validation error where selected, resulting in: K = 3, λ = 2.5 ⋅ 10⁻⁶, and sgd.

Interestingly, only around 42 units per hidden layer have non-zero connections after training. Although, if started with only 50 units we observed worse performance. The mean of the square test error (over locations) is 0.01387 with standard deviation 0.00304.

Classifiers

For classification purposes we assign each time bin to one of two classes: “fluctuating” or “constant”. “Fluctuating” corresponds to discs moving over the site of interest and decreasing the light intensity in that site, while “constant” refers to the constant illumination of the site when no discs are present. To label the time bins we use a simple cut-off criterion plus two further correcting steps to account for retinal adaptation effects. First we label as “fluctuating” every bin with stimulus intensity less than 0.99. Then we apply these corrections: i) Every identified “constant” segment shorter than 30 bins (375 ms) is relabelled as “fluctuating,” and ii) The first 30 bins following a “fluctuating” segment are also labelled “fluctuating.” In this way the stimulus at each site is divided in segments of fluctuating and constant intensity. We train both linear and nonlinear Support Vector Machine (SVM) classifiers to determine, from the spike train response, whether a given time bin is labelled as “constant” or “fluctuating”. Similarly to the decoding framework, to classify a given bin we consider a time window of ΔT = 61 bins around it in the response. For the nonlinear SVM we use the same gaussian kernel as in nonlinear decoding and the parameter values obtained when training the decoder. Note that this is not the optimal nonlinear classifier but allows us to evaluate the classifying power of the decoding kernel.

Measures of performance

Given a stimulus intensity trace and the corresponding decoding prediction we define the decoding error as the Mean Squared Error We also make use of the related Fraction of Variance Explained defined as FVE = 1 − (MSE/Var(y)).

To measure decoding performance from the fully decoded movie we build Receiver Operating Curves (ROC). We threshold the decoded intensity trace at each site. If intensity is below threshold, the presence of a disc in the site is predicted. By comparing the prediction as a function of the threshold to the original stimulus frames (where the center of the site can only be white when no disc is present, or black when the disc is present), we can evaluate the performance of the decoder as a balance between the True Positive (TP) and False Positive (FP) rates

To assess the performance of the SVM classifiers we use the F₁-score measure defined as where P is the Precision and R the Recall given by For the binary classification task, “fluctuating” is defined as the positive class.

Unless otherwise stated, all of the statistical significance tests were performed with the Wilcoxon signed rank test.

ON/OFF ratio bias estimation

For each site s we determine the set of available cells whose receptive field centers were less than 300 μm distant from the center of the site. We call C_s the total number of available cells at site s. In general, C_s is the sum of ON and OFF subtype cells, . If, from the available cells at site s, we pick a random subset of size N = N^on + N^off, the probability of choosing N^off cells is given by the hypergeometric distribution (random draw without replacement) The average probability over all sites considered is

Separately, for each site s we have established a hierarchy of cells from their decoding filter norms. Following the hierarchy we create decoding sets of different size N (the best cell, the best two cells, etc) and we count the number of OFF type cells N^off in them. We summarize this information in the histogram M(N^off, N) that counts the number of sites where the decoding set of size N contains N^off OFF cells. With this histogram we obtain an empirical probability that we can compare with p(N^off|N). In particular, the bias reported in S5 Fig is given by Only sites with N^off, N^on ≥ 2 were considered for the comparison (n = 115).

Encoding model

We build an encoding model for a single cell, based on the standard GLM type model proposed by Pillow et al [15]. The cell spikes stochastically through a Poisson process with a time-dependent firing rate λ(t) given by where is a spatio temporal filter acting on stimulus and is a temporal filter of the past spike history of the cell represented by . The function f(x) is a rectifying nonlinearity of the log-exp form f(x) = a log(b exp(x + c)). The stimulus filter factorizes into separate spatial and temporal filters. The spatial component is given by a balanced difference of gaussians, with widths σ_c = 35μm for the positive and σ_s = 100μm for the negative part, providing a symmetrical center-surround type filter. The temporal part of the filter is given by a single negative lobe of a sin-like function. The filter for the past spike history takes the form This filter inhibits firing after a spike but, depending on the values of the parameters, it can have a positive lobe after the inhibitory part that tends to increase the firing rate. We consider a span of 250 ms (20 bins) for both the past history filter and the temporal part of the stimulus filter. All elements of the filter are fixed except for the rectifying nonlinearity that is changed according to the value of α. Initially, the parameters of the nonlinearity f_{α = 1}(x) are adjusted to provide an average firing rate similar to that observed in real data. The α = 1 model is taken as the ground-truth and every time α changes, the nonlinearity f_α(x) is fitted anew by maximizing the likelihood on α = 1 rasters, in order to reproduce the firing rate trace (PSTH) as closely as possible to the PSTH generated by α = 1. The model neuron is stimulated with real data and the intensity trace at the central site of its receptive field is the stimulus considered for decoding. The model has been implemented using the Nonlinear Input Model toolbox [31].

Decoding setup

We recorded the spiking activity of C = 91 ganglion cells from a 1 mm² patch of the rat retina, while presenting a complex and dynamical stimulus that consisted of 1, 2, 4 or 10 black discs on a bright background (Fig 1A and Methods). The discs followed self-avoiding random motion, generated as described in the Methods section, which (for decoder training) was non-repeated; all decoding results are reported on withheld (test) segments of the stimulus that were not used during training. The stimulus also contained a segment of repeated trajectories that was randomly interspersed into the non-repeated part and used only to assess the role of noise correlations. Our goal was to reconstruct the light intensity as a function of time at a grid of 20 × 20 spatial positions (“sites”) uniformly tiling the stimulus frame. Specifically, at each site, we convolved the original movie with a small Gaussian filter (see Methods), which defined the “luminance trace” at that site, to be decoded. Stimulus features (here, disc size) were smaller than the receptive field center of a typical recorded RGC, making the decoding task non-trivial.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 1. Linear decoding of a complex movie.

A: An example stimulus frame. At each site (red dots = partially shown 20×20 grid) the stimulus was convolved with a spatial gaussian filter (red circle = 1σ). Typical RGC receptive field center size shown in gray. B: Responses of 91 RGCs with 750 ms decoding window overlaid in blue. C: Three example luminance traces (red) and the linear decoders’ predictions (blue). D: Decoded frame (same as in A) reconstructed from 20×20 separately decoded traces. Disc contours of the original frame shown for reference in green. E: RF centers of the 91 cells (black dots = centers of fitted ellipses). RF centers overlapping a chosen site (red dot) are highlighted in blue. F: Performance of the linear decoders across space, as Fraction of Variance Explained (FVE). Black dots as in E; black contour is the boundary FVE = 0.4. G: Performance of the linear decoders (FVE) across sites as a function of cell coverage (grayscale = conditional histograms, red dots = means, error bars = ± SD). H: Average decoding error across sites (MSE ± SD) of 10-disc-trained decoders, tested on withheld stimuli with different numbers of discs. I: Cells (black dots = RF center positions) contributing to the decoding at two example sites (red circles); decoding filters shown below. For each site, contributing cells (highlighted in red and joined to the site) account for at least half of the total L1 norm. J: Decoding field of a single cell (here, evaluated over a denser 50×50 grid and normalized to unit maximal variance); the cell’s RF center shown in black.

https://doi.org/10.1371/journal.pcbi.1006057.g001

Sparse linear decoding of a complex movie

To estimate the luminance trace at any given time, we trained a separate sparse linear decoder for each site on a 750 ms sliding window of the complete spiking raster, shown in Fig 1B, and represented as spike counts in Δt = 12.5 ms time bins (see Methods). The decoder minimized the square error between the true and estimated luminance trace at each site, using sparse (L1) penalty on decoding weights, as specified by Eq (2). While each decoder in principle had access to all neural responses, the sparse penalty ensured that the majority of the weights corresponding to redundant or non-informative neural responses for each site were zeroed out, yielding interpretable results which we describe in detail below. When trained on the 10-disc stimulus, this procedure predicted well the luminance traces across individual sites on withheld sections of the stimulus (Fig 1C), allowing us to reconstruct the complete movie (Fig 1D).

We expected the performance of our decoder to depend strongly on local coverage, i.e., on the number of recorded cells whose receptive field centers overlap a given site. Coverage amounted to about six cells on average and exhibited substantial spatial heterogeneity, as shown in Fig 1E. The quality of our movie reconstruction, measured locally by “fraction of variance explained” (FVE, see Methods), showed similar spatial variation (Fig 1F) which correlated strongly with coverage (Fig 1G), and saturated at ≥ 6 cells. In what follows, we restrict our analyses to sites with good coverage that pass a threshold of FVE ≥ 0.4. Despite the high dimensionality of this regression problem (decoders have ∼5 ⋅ 10³ parameters per site), sparse regularization ensured uniformly good performance even when tested on out-of-sample stimuli with varying number of discs (Fig 1H).

To analyze how rich stimuli are represented by a population of ganglion cells with densely overlapping receptive fields, we examined the resulting decoding weights in detail. We found that stimulus readout was surprisingly local. As illustrated for two example sites in Fig 1I, only a few cells whose receptive field centers were in close proximity to the respective sites were assigned non-negligible decoding weights. This was true in general: on average 5.4 ± 2.8 cells, whose RF centers were all located within 200 μm of the decoded site, contributed to the luminance trace reconstruction; cells beyond this spatial scale contained no decodable information (S1 and S2 Figs).

Our framework also allowed us to construct a “decoding field” for every cell (Fig 1J). A decoding field represents an impulse response of the decoder, i.e., an additive contribution to the stimulus reconstruction for every spike emitted by a particular cell. While one can reasonably expect that the receptive and decoding fields overlap in location and spatial extent, there is no theoretical guarantee that this must happen, given that neural encoding is strongly nonlinear. We nevertheless confirm this expectation and observe a very good correspondence between the spatial locations and sizes of the decoding and receptive fields for all cells (S3 Fig), with decoding fields also exhibiting a clear center-surround-like structure. We find that decoding filters shapes for all cells are highly stereotyped (S4 Fig). We further find that the readout of retinal responses is local, in the sense that only cells with receptive field centers close to the decoding site contribute to the decoding (Fig 1I, S6 Fig). Lastly, the readout is structured, in the sense that the cells that contribute to decoding at each site have a preferred ON vs OFF composition that favors recruiting OFF cells (S5 Fig), most likely because the visual feature that moves in our stimulus is a dark disc on a bright background.

Taken together, our results suggest that retinal responses to complex stimuli can be read out in a highly stereotyped, structured, and local manner.

Nonlinear decoding outperforms linear decoding

Could nonlinear decoding improve on these results? We considered two nonlinear regression methods that can tractably be applied to our data: kernel ridge regression (KRR) and regression using deep neural networks.

Kernel ridge regression is a well-known extension of linear regression into the nonlinear domain by means of the kernel trick. We used Gaussian kernels whose width was determined using cross-validatation (see Methods) [29], as specified in Eqs (3)–(6). Importantly, the success of this nonlinear decoder crucially depended on the proper selection of local groups of cells relevant for each site, as identified by linear decoding: its sparse (L1) regularization acted as “feature selection” for the nonlinear problem (Methods, S7 Fig). Nonlinear decoder could then make use of higher-order statistical dependencies within and between the selected spike trains to achieve high performance.

We compared these results to regression using neural networks. An architecture that achieved good performance consisted of an input layer (that received spiking rasters of the same dimension as the linear regression problem), followed by three fully-connected hidden layers with 150 sigmoidal neurons each, followed by a 20 × 20 output layer whose units correspond to the decoding sites of our movie; this architecture is schematized in S20 Fig. The network was trained by minimizing the squared reconstruction error, Eq (7), using standard deep learning tools (see Methods).

Fig 2A shows a luminance trace at one of the example sites, together with its linear and nonlinear reconstruction. Nonlinear decoders track better the detailed structure of luminance troughs, which occur when discs cross the site, as well as exhibiting smaller fluctuations when no discs are crossing the site and the true luminance trace is therefore constant. This is reflected in a substantial overall increase in fraction of variance explained (FVE) across different sites, shown in Fig 2B. A kernelized nonlinear decoder using only two best cells per site outperforms, on average, the best sparse linear decoder constructed from the entire population; nonlinear performance saturates quickly with the number of cells and peaks when decoding from local ∼8-cell groups. The neural network decoder, which we train on the complete neural population, reaches a comparable performance to the best kernelized decoder.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Nonlinear decoding outperforms linear decoding.

A: Luminance trace (red) with linear (blue) and nonlinear KRR (green) and neural network (grey) predictions. B: Average decoder performance (± SD across sites), achievable using increasing numbers of cells with highest L1 filter norm. For nonlinear decoding, “All” is the optimal subset that maximizes performance (S7 Fig). Since the neural network (grey point with an error bar) simultaneously decodes the movie at all sites, it only makes sense to train it using “All” cells. C: Average ROC across all testing movie frames. D: Fractional improvement (average ± SEM across sites) of nonlinear KRR versus linear decoders for test stimuli with different numbers of discs. All decoders were trained only on the 10-disc stimulus. E: Decoding error (MSE; average ± SEM across sites) in fluctuating and constant epochs is significantly larger for linear decoders (p<0.001) relative to nonlinear KRR and the neural network.

https://doi.org/10.1371/journal.pcbi.1006057.g002

An alternative way to compare decoding performance is to threshold the sequence of decoded movie frames (see S8 Fig and S1 Video), thereby assigning each site to a decoded dark disc (“below threshold”) or to the bright background (“above threshold”). Decoded movie frames can then be compared to ground truth (i.e., the original movie frames which can only be either black or white at every location) at each threshold using the receiver operator characteristic (ROC curve), shown for all decoders in Fig 2C. In this metric, the performance of the kernelized and neural network decoders are nearly indistinguishable, and consistently outperform linear decoders. Excess nonlinear performance of between 15 and 20% of FVE was maintained even when decoders were trained on 10-disc stimulus and tested on stimuli with smaller number of discs (Fig 2D). Excess nonlinear performance was also observed when decoding from a cell mosaic of a single functional type (S9 Fig) and on a repeat experiment (S10 Fig).

We note the surprising consistency between the kernelized decoding and neural network results. Despite the fundamental differences in the nature and application of these two regression methods—neural networks are universal approximators, use different regularization from the kernelized decoders, and have been trained on all cells simultaneously to decode at all sites simultaneously, in contrast to the kernelized decoders—their numerical measures of performance appear quantitatively consistent. While it is impossible to exclude the possibility that yet another type of decoder could yield much higher performance, it is also possible that both nonlinear decoders we tried managed to extract all available information about local luminance traces from the recorded spike trains.

Another particularly striking feature of our results was the difficulty of the linear decoder to match the true (constant) luminance trace when no disc was crossing the corresponding site. Rat retinal ganglion cells are continuously active even when there are no coincident on-center luminance changes, with the activity likely resulting from stimulus changes in the surround, from long-lasting sustained responses to previous stimuli, from effective network coupling to cells that do experience varying input, or from true spontaneous excitation that would take place even in complete absence of stimuli [32–35]. Either way, activity of cells at constant local luminance presents a confound that is difficult for a generic linear mechanism to eliminate, which results in decoder fluctuations, or “hallucinations,” of sizeable variance. To quantify this effect, we partitioned the luminance traces at every site into constant and fluctuating epochs by means of a simple threshold (see Methods), and examined decoding errors separately during both epochs. The decrease in decoding error by using nonlinear decoders was similar in absolute terms in both epochs, but represented a much larger fractional decrease during constant epochs, suggesting that nonlinear decoders might specifically be better at suppressing their responses to spontaneous-like neural activity (Fig 2E).

We reasoned that this improvement comes, in part, from the ability of both nonlinear methods to recognize whether there are any on-site luminance fluctuations or not, from the spike trains alone. To test this idea, we trained linear and nonlinear kernelized classifiers, operating on identical inputs and with the same kernel parameters as the decoders, to best separate constant from fluctuating activity. Consistent with our expectations, nonlinear classifiers outperformed linear at every site, irrespectively of whether their input were the rasters of all local cells that contribute to the decoding, as shown in S11A Fig, or the raster of a single best cell at every site, as shown in S11B Fig.

Nonlinear decoders make use of spike-history dependencies in individual spike trains

Next, we attempted to identify the statistics of spike trains that are necessary to explain the excess performance of nonlinear decoders. Our starting point was the following observation: the simplest nonlinear decoders that used a single best cell for each site, when interrogated with a test-set epoch of pure spontaneous activity (i.e., neural responses to a completely blank screen), yielded luminance traces with significantly smaller variance than their linear counterparts (S12 Fig). Since the only structure in spike trains during spontaneous activity is, by definition, due to “noise correlations”—pairwise or higher-order dependencies between spikes within an individual spike train or across different spike trains—we hypothesized that certain noise correlations could be used by nonlinear decoders also during stimulus presentation to boost their decoding performance.

To test this hypothesis, we made use of many identical repeats of a particular stimulus fragment embedded in our disc movie (these repeats were used neither for training nor testing). Using the same decoders as above, we decoded the original response rasters corresponding to the repeated fragment, as well as rasters in which we shuffled the spikes to remove spike-history dependencies, or to remove cell-cell noise correlations, as shown in Fig 3A, to assess how decoding is impacted by the removal of certain types of correlation in the spike trains. Note that these manipulations left the firing rates of all cells intact, and thus preserved all correlations in the spike trains that are due to the neurons responding to a spatio-temporally correlated stimulus. Crucially, for our analysis we did not retrain our decoders on the shuffled spike trains, because we wanted to ask whether the same decoders that we trained to read the real (unshuffled) neural code can also read the modified neural code lacking various components of the noise correlations. If the decoder performance were unaffected by such removal, then noise correlations are not crucial for our decoder; in contrast, a drop in decoder performance would suggest that noise correlations may be necessary. Alternatively, if we were to retrain our decoders on the shuffled spike trains, we would be answering a different question: Is there any decoder that can read the shuffled neural code? While interesting, (i) it is unclear what statements about the actual neural code such an analysis would provide (since these decoders would be trained on synthetic, shuffled codes that only exist in our computer); (ii) technically, the number of distinct training samples would be drastically too small to train such decoders, since the experiments impose a hard trade-off between the number of repeats and the duration of the repeated fragment, on which the decoder would have to be trained. For these scientific and methodological reasons, we performed the following analysis using decoders trained on actual (unshuffled) responses to unrepeated stimuli.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Spike-history dependencies affect decoding performance.

A: Shuffles of responses to repeated stimulus presentations remove different types of correlations, but preserve average locking to the stimulus (PSTH), and thus stimulus-induced correlations. B: A repeated stimulus fragment (red trace), nonlinear kernelized decoder predictions using real responses (green), and using responses without different types of correlations (gray); shown is the prediction mean ± SD over repeats. C: Increase in decoding error (MSE) when spike-history dependencies or noise correlations are removed (average ± SEM across sites); percentages report fractional differences relative to the original performance. D: Spike count distributions for a single example cell. Removing spike-history dependencies broadens the distributions, in particular in constant epochs. Dashed line = expectation for a fully randomized spike train with a matched firing rate. E: Variance-to-mean ratio F of spike count distributions for spike trains with and without spike-history dependencies. Each point is a cell that contributes most to decoding at a particular site (when the same cell contributes to multiple sites, average ± SD across sites is shown).

https://doi.org/10.1371/journal.pcbi.1006057.g003

Fig 3B shows a stimulus reconstruction at an example site by the nonlinear kernelized decoder, for original rasters as well as rasters with removed spike-history dependencies or cell-cell noise correlations. Removing cell-cell noise correlations leads to a small increase in the variance of the reconstructions across stimulus repeats, with only marginal differences in the mean reconstructed trace, compared to decoding from intact rasters. Surprisingly, removing spike-history dependencies leads to much worse reconstructions, whose mean is strongly biased and variance increased; as a result, the dynamic range of the decoded trace is substantially lower compared to decoding from intact rasters. These observations are summarized across sites in Fig 3C, which shows the increase in decoding error when spike-history dependencies or cell-cell noise correlations are removed. Removal of cell-cell noise correlations leads to small increases in error, roughly of the same magnitude for both linear and nonlinear decoders; in contrast, while removal of spike-history dependencies leads to increases in error for both decoders, the effect is two-to-three-fold larger for the two nonlinear methods. We emphasize that kernelized decoders and neural network are two fundamentally different regression methods, yet the removal of spike-history dependences strongly decreases the performance of both, suggesting that our observations are likely not a consequence of choosing a particular decoder type. Qualitatively similar conclusions hold for the classifiers trained to separate constant from fluctuating input epochs (S13 Fig), as well as for decoders and classifiers trained on the single best cell per site (S14 Fig).

Having established that spike-history dependencies are crucial to the performance of the nonlinear decoders, we looked at the detailed statistical structure of individual spike trains. For each neuron that best decoded the luminance trace at a specified site, we focused on 250 ms (20 time bins) response sequences and constructed a distribution over the number of occupied time bins (“spike counts”), separately for epochs where the luminance trace was fluctuating or where it was constant. As shown in Fig 3D, these distributions differed significantly: the count distribution was much tighter in constant epochs, while the mean firing rate between the epochs did not change much. During fluctuating-input epochs, observing more spikes in a 250 ms window was more likely than at constant input, but—perhaps surprisingly—patterns with very low numbers of spikes (e.g., zero or one) were also more likely during fluctuating-input epochs. The count distribution at fluctuating light was very similar to binomial (and, at this temporal resolution, Poisson), while it was tighter at constant light. These changes could be summarized by a simple statistic, the variance-to-mean ratio F = (variance in spike count)/(mean spike count). Note that unlike the standard Fano factor, our F is not computed across the repetitions of the same stimulus and thus measures the total variability in the response, which includes variance due to the changing stimulus. When we removed spike-history dependencies, the variance-to-mean ratio F increased for both distributions and they became harder to distinguish from each other. Fig 3E shows that this behavior was consistent across all sites, highlighting the very high regularity of neural spiking that resulted in sub-Poisson variance (F substantially below 1) during epochs of constant luminance.

How could spike-history dependencies help in stimulus decoding? A possible scenario would involve the situation where decoders should “sum” multiple spikes from the same neuron in the recent past super-linearly, to optimally reconstruct the stimulus. In this case, without spike-history dependencies that are responsible for precise firing with sub-Poisson variance, the Poisson spiking in the absence of spike-history effects would cause large (compared to linear decoder) spurious variance in the decoder output. Adding spike-history dependencies would, in contrast, tighten the number of emitted spikes, giving the nonlinear decoder a reliable option to sum spike effects super-linearly without being swamped by spiking noise. We emphasize that this is only the simplest scenario we could think of as an example where spike-history dependencies could be beneficial; there are likely many others.

Taken together, our results show that: (i), spike-history dependencies within individual spike trains are crucial for nonlinear decoder performance; (ii), these dependencies shape the distribution of spike counts on timescales relevant for decoding; (iii), during constant local luminance, spiking activity is very regular (and statistically similar to true spontaneous activity, see S15 Fig); (iv), a simple statistic, which summarizes the effects of spike-history dependencies in different epochs and their changes when the spike trains are shuffled, is the spike variance-to-mean ratio F. This does not imply that nonlinear decoders actually compute some version of a local estimate for F: they could be sensitive to other statistics, e.g., the interspike interval distribution, which also differs substantially between the epochs, see S16 Fig. Because nonlinear decoders we use have no well-defined set of sufficient statistics, it is impossible to claim which precise statistic of the spike train they are sensitive to, beyond stating that they clearly are sensitive to the removal of spike-history dependencies. Note further that we can only establish clearly that nonlinear decoders that we trained are sensitive to the removal of spike-history dependencies; we, however, cannot exclude the option that there exist nonlinear decoders of the same class that reach similar performance as ours but are at the same time robust to the removal of spike-history dependencies. Nevertheless, subsequent analyses on synthetic data that we provide below, as well as the robustness of our observations with respect to the nonlinear method (kernelized decoder and the neural network) suggest that crucial decoding information really is present in the spike-history dependencies, and that the underlying reason for nonlinear decoder performance is its ability to recognize high regularity of spiking during epochs of constant local luminance.

A simple neural encoding model can recapitulate spike train statistics crucial for nonlinear decoding

Can the observed spike-history dependencies, which enable successful nonlinear decoding, be generated by simple and generic neural encoding models? To address this question, we made use of generalized linear models (GLMs) [36, 37], probabilistic functional models of spiking neurons that extend the paradigmatic linear-nonlinear (LN) framework by incorporating the recurrent feedback from neuron’s past spiking, as schematized in Fig 4A. Previously, GLMs have been successfully applied to responses of the mammalian retina [15, 20] and in the cortex [38, 39], and also reproduced well the firing rates of cells recorded in our experiment on the repeated stimulus fragment (S17 Fig).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Spike-history dependencies of intermediate strength facilitate nonlinear decoding in simple models of neural processing.

A: Schematic of a single-cell Generalized Linear Model (see Methods). The neuron’s sensitivity to the stimulus is determined by a radially symmetric difference-of-Gaussians spatial filter that has a monophasic timecourse (), and combines additively with the neuron’s sensitivity to its own past spiking, given by filter (with strong refractoriness followed by weak facilitation). Importantly, shapes spike-history dependencies in the resulting spike trains. A nonlinear function f(⋅) (here, threshold-linear) of the combined sensitivities gives the neuron’s instantaneous firing rate that can be used to generate individual spike train instances. Shapes, as well as the temporal and spatial scales of the filters, were realistic for our data. B: Example rasters (50 repeats) generated with the encoding model for a given intensity trace and different magnitudes (α) of spiking history filter . The rasters are matched in PSTH (bottom) but differ in temporal noise correlations. C: Average spike count variance-to-mean ratio, F, (± SD) of the model as a function of α in fluctuating and constant epochs. D: Decoding error as a function of α. Decoders are trained for each separate α and tested on withheld stimuli; shade = SD over 10 spike train realizations.

https://doi.org/10.1371/journal.pcbi.1006057.g004

To link encoding models and decodability in a way that would generalize beyond the specifics of our dataset, we created the simplest stereotyped model cell, shown in Fig 4A. Crucially, we parametrized the magnitude of the self-coupling filter with α: α = 0 thus corresponded to a pure LN model, while increasing values of α made neural spike trains non-Poisson, progressively enforcing dependence on past spiking and consequently increasing the magnitude of the resulting temporal correlations.

With this model in hand, we generated a “baseline” raster of repeated responses to a randomly moving disc stimulus at an initial value of α = 1, which corresponds to the strength of spike-history dependence inferred from our data, as shown in Fig 4B. The average firing rate was chosen to be the typical rate of our recorded ganglion cells. We then systematically changed the value of α and, for each value, refitted the nonlinearity by maximizing the likelihood to the baseline raster at α = 1 (see Methods). This procedure generated synthetic rasters that, to an excellent approximation, were matched in their peri-stimulus time histograms (PSTH) and stimulus preference, yet differed in the strength of spike-history dependencies.

Following our previous analyses, we partitioned the luminance trace into constant and fluctuating epochs, and looked at the spiking statistics in 250 ms (20 time bin) windows. Spike count variance-to-mean ratio F in constant epochs decreased as a function of α and dropped substantially below 1; in contrast, when on-center luminance was fluctuating, F behaved non-monotonically (Fig 4C). In line with expectations and behavior observed in our data, F at constant luminance was always below F at fluctuating luminance. Having ensured that the statistics of synthetic rasters qualitatively agreed with the data for the range of α we examined, we asked about the performance of linear and nonlinear decoders, trained and tested at different values of α. Fig 4D plots the decoding error as a function of α; see S19 Fig for the separation into error in fluctuating vs constant epochs as a function of α. Overall, the error levels are in range of those observed for real data (cf. Fig 1H), with nonlinear decoders outperforming linear by ∼10 − 30%. Interestingly, the minimal error for both decoders is achieved at an intermediate value of α* ≈ 0.4, which also corresponds to the point where nonlinear decoders maximally outperform their linear counterparts. At α = 0, where the encoding models are effectively LN neurons, the decoders differ only marginally in performance (analogous results hold for the classifiers separating fluctuating from constant epochs, see S18 Fig).

How close are real retinal ganglion cells to the value of α that permits best nonlinear reconstruction? This question cannot be answered precisely with the toy models we use. Within the class of generalized linear models (GLMs) considered here, we can reliably show that nonlinear decoding performance significantly outperforms linear decoding performance for a broad range of α values that includes both α* ≈ 0.4 as well as α = 1 (which, by our definition, corresponds to the best fit of GLM model to our data), and this effect is robustly true for all the cells that we examined. It is, however, likely that GLM models are too simplistic for the cells we are considering (realistic models for rat ganglion cells may require two nonlinear stages of stimulus processing, i.e., LNLN models [40]), if we wanted to make a quantitative statement about how close real cells are to the value of α that permits optimal nonlinear reconstruction. In these simulations we also haven’t modeled cell-to-cell noise correlations; further, it is likely that GLM does not capture all spike-history dependencies; and we decoded only the central pixel of the model’s receptive field. These differences between the simulations and the real experiment are likely responsible for the fact that the difference between nonlinear decoding performance from GLM-generated spike trains and same spike trains with shuffle-removed spike-history dependencies are much smaller than what we see in real data. Nevertheless, while the quantitative match between real neural data and GLM simulations is beyond the scope of this paper, we have shown qualitatively that in a generic class of encoding models that have been widely applied to both peripheral as well as central neural processing, there exists a non-trivial strength of spike-history dependence that facilitates nonlinear stimulus reconstruction. Intuitively, the existence of optimal α* > 0 can be explained as a trade-off between ensuring regularity of spiking during constant epochs, which the nonlinear decoder can make use of, while not impeding stimulus encoding during fluctuating epochs; during these epochs, stimulus-driven term should dominate over sensitivity to past spiking, otherwise excessive dependence on spiking history (e.g., α ≥ 1 in Fig 4B) could perturb reliable locking to the stimulus.