Extreme synergy: Spatiotemporal correlations enable rapid image reconstruction from computer-generated spike trains

Abstract

Over the brief time intervals available for processing retinal output, the number of spikes generated by individual ganglion cells can be quite variable. Here, two examples of extreme synergy are used to illustrate how realistic long-range spatiotemporal correlations can greatly improve the quality of retinal images reconstructed from computer-generated spike trains that are 25–400 ms in duration, approximately the time between saccadic eye movements. Firing probabilities were specified both explicitly: using time-varying waveforms consistent with stimulus-evoked oscillations measured experimentally, and implicitly: by superimposing realistic fixational eye movements on a biophysical model of primate outer retina. Synergistic encoding was investigated across arrays of model neurons up to 32 × 32 in extent, containing over 1 million pairwise correlations. The difficulty of estimating pairwise, spatiotemporal correlations on single trials from only a few events was overcome by using oscillatory, local multiunit activity to weight contributions from all spike pairs. Stimuli were reconstructed using either an independent rate code or the first principal component of the single-trial, pairwise correlation matrix. Spatiotemporal correlations mediated dramatic improvements in signal/noise without eliminating fine spatial detail, demonstrating how extreme synergy can support rapid image reconstruction using far fewer spikes than required by an independent rate code.

Introduction

The conditions under which retinal ganglion cells transmit visual signals synergistically remains a topic of considerable debate (Latham & Nirenberg, 2005; Schneidman, Bialek, & Berry, 2003; Schnitzer & Meister, 2003). Measurements of synergy between small groups of neurons have found evidence for everything from redundancy (Gawne & Richmond, 1993; Puchalla, Schneidman, Harris, & Berry, 2005; Warland, Reinagel, & Meister, 1997) to statistical independence (Nirenberg, Carcieri, Jacobs, & Latham, 2001; Panzeri, Schultz, Treves, & Rolls, 1999; Reich, Mechler, & Victor, 2001; Rolls, Franco, Aggelopoulos, & Reece, 2003) to both modest (Dan, Alonso, Usrey, & Reid, 1998) as well as more substantial levels of synergy (deCharms & Merzenich, 1996; Gat & Tishby, 1999; Hirabayashi & Miyashita, 2005; Riehle, Grun, Diesmann, & Aertsen, 1997; Samonds, Allison, Brown, & Bonds, 2004; Singer, 1999; Vaadia et al., 1995). Several studies have looked specifically for synergy among larger ensembles (Bezzi, Diamond, & Treves, 2002; Frechette et al., 2005; Narayanan, Kimchi, & Laubach, 2005; Stanley, Li, & Dan, 1999), yet it remains an open question as to how a nonlinear code based on spatiotemporal correlations between hundreds of neighboring ganglion cells—representing hundreds of thousands of neuron–neuron pairs—might convey local pixel-by-pixel intensity information that could not be obtained by analyzing the same spike trains individually, especially over the short time scales—approximately 50 to 300 ms—available for interpreting retinal output (Bacon-Mace, Mace, Fabre-Thorpe, & Thorpe, 2005; Kirchner & Thorpe, 2006; Rolls, Tovee, & Panzeri, 1999; Thorpe, Fize, & Marlot, 1996). 

Here, the effectiveness of two extremely synergistic codes, measured relative to an independent rate code, was quantified using the quality of image reconstructions mediated by short sections of massively parallel spike train data. As a baseline, rate-based reconstructions employed only the number of spikes from each model ganglion cell to estimate the local intensity at each pixel. Synergistically encoded images, on the other hand, used Principal Components Analysis (PCA) to estimate local pixel intensities from the nonlinear spatiotemporal correlations among large numbers of ganglion cell pairs. Two complementary methods for generating artificial spike trains with realistic spatiotemporal correlations were considered. 

Common oscillatory input

First, realistic spatiotemporal correlations were generated by simultaneously modulating the instantaneous firing probabilities of all spike trains underneath the stimulated (foreground) cells, using a common, oscillatory waveform whose amplitude and temporal structure was consistent with both optic tract recordings (Doty & Kimura, 1963; Laufer & Verzeano, 1967; Steinberg, 1966) and with single and multiunit correlation functions (Ishikane, Kawana, & Tachibana, 1999; Neuenschwander & Singer, 1996). By specifying the amplitude, time course, and extent of the spatiotemporal correlations explicitly, it was possible to quantify how such factors affected the fidelity with which synergistically encoded images could be reconstructed using only the information available in short sections of simultaneously recorded spike trains. 

Fixational eye movements

To explore extreme synergy in the context of a very different yet still physiologically plausible mechanism, spatiotemporal correlations were instead produced implicitly by superimposing small fixational eye movements on a biophysically realistic model of the primate outer retina. Detailed representations of cone photoreceptors were combined with reduced descriptions of bipolar cells and horizontal cells, the latter modified to include local electrical coupling (van Hateren, 2005). Fixational eye movements were produced in accord with previously published models derived from measurements of drift and tremor in human subjects (Eizenman, Hallett, & Frecker, 1985). Parallel spike trains were obtained by using the membrane potentials of the model bipolar cells to estimate an instantaneous firing probability for the overlying ganglion cells, consistent with the nearly one-to-one mapping of bipolar cells to midget ganglion cells in the primate fovea (Calkins, Schein, Tsukamoto, & Sterling, 1994). 

Both the explicit and implicit correlation models were used to generate short segments of massively parallel computer-generated spike trains, from 25 to 400 ms in duration, simulating a retinal patch containing up to 1000 output neurons (equivalent to 1 million pairwise correlations). All model ganglion cells were assumed to be of the same type and their receptive field centers were arrayed in a square rectilinear grid. Background and stimulated firing rates were chosen so as to encompass much of the measured dynamic range of cat retinal ganglion cells in response to diffuse, sinusoidally varying stimuli (Troy & Enroth-Cugell, 1993). 

Our findings indicate that over a nearly 16-fold range of stimulated firing rates and without eliminating fine spatial detail, realistic pairwise spatiotemporal correlations over extended neighborhoods can support rapid pixel-by-pixel reconstructions of visual stimuli that are qualitatively and quantitatively superior to analogous reconstructions obtained simply by counting the number of spikes from each neuron, even when the reduced trial-to-trial variability of retinal ganglion cells elicited by high-contrast stimuli is taken into account (Kara, Reinagel, & Reid, 2000). The present findings thus imply that realistic spatiotemporal correlations can greatly augment the stimulus information conveyed by a simple rate code, in which local pixel intensity is conveyed solely by the number of spikes produced by the overlying ganglion cell. Neither purely spatial correlations in the number of spikes, nor purely temporal correlations within the individual spike trains, could by themselves account for the superior quality of the image reconstructions obtained by fully exploiting both spatial and temporal correlations simultaneously. Spatiotemporal correlations mediated performance levels on an ON/OFF pixel discrimination task that, if instead mediated by an independent rate code, would have required approximately four times as many spikes to achieve similar accuracy. By distributing local intensity information across an extended neighborhood, the present study suggests that rapid image reconstructions can be accomplished using far fewer spikes than would otherwise be required. Of course, the question of whether extreme synergy exists among retinal ganglion cells, and whether such encoding is utilized by the visual cortex, can only be answered experimentally. The present theoretical study may be useful, however, in motivating additional experimental investigation. 

Methods

Computer-generated spike trains I: Common oscillatory input

Overview

The procedures used here for generating spatiotemporal correlations via common oscillatory input are similar to those previously described (Kenyon, Theiler, George, Travis, & Marshak, 2004; Stephens, Neuenschwander, George, Singer, & Kenyon, 2006). Starting with a power spectrum containing a single Gaussian peak and characterized by three parameters, specifying the amplitude, width, and central frequency of the oscillatory modulation, an oscillatory waveform was obtained by transforming back to the time domain after randomly choosing the phases of the individual Fourier components. Parallel sets of artificial spike trains with realistic spatiotemporal correlations were constructed using an array of binary event generators, employing the same common waveform to simultaneously modulate their instantaneous firing probabilities. Although on any given trial, every cell activated by the stimulus had the same time-dependent firing probability, the spike trains themselves varied considerably from cell to cell due to the underlying stochastic process. 

Mathematical details

An oscillatory time series of duration T and temporal resolution Δ t was constructed by first defining the discrete frequencies f _k:  

$f_k=\frac{k}{T},0\le k<\frac{T}{\Delta t},$

(1)

in terms of which the discrete Fourier coefficients were defined as follows:  

$C_k=\exp \left(\frac{-(f_k-f_0{)}^2}{2{\sigma }^2}\right)\exp (2\pi i{r}_k),$

(2)

where f ₀ is the central oscillation frequency, here set to 80 Hz, σ is the width of the spectral peak in the associated power spectrum, here set to 10 Hz, and r _k is a uniform random deviate between 0 and 1 that randomized the phases of the individual Fourier components (generated by the Matlab intrinsic function RAND). The procedure was then repeated using a different pseudorandom sequence for each stimulus trial. The coefficients, C _k, were used to convert back to the time domain using the discrete inverse Fourier transform (generated by the Matlab intrinsic function IFFT), so that the firing rate R _n ^{( i)} at each time step, t _n = nΔ t, was given by  

$R_n^{\left(i\right)}=I^{\left(i\right)}\frac{1}{N}\sum _{k=1}^{N-1}{C}_k{e}^{-2\pi i{f}_k{t}_n}+F_0^{\left(i\right)},$

(3)

where the R _n ^{( i)} depends only on the real part of the sum on the RHS of Equation 3, I ^{( i)} and F ₀ ^{( i)} are scale factors used to encode the stimulus intensity—denoted by the superscript i, and N = T/Δ t. 

Once determined, the time series defined by R _n ^{( i)} was used to generate oscillatory spike trains via a pseudorandom process:  

$S_n^{\left(i\right)}=\theta (R_n^{\left(i\right)}\Delta t-r),$

(4)

where R _n ^{( i)}Δ t is the probability of a spike in the nth time bin, θ is a step function, θ( x < 0) = 0, θ( x ≥ 0) = 1, and r is again a uniform random deviate. In the limit that R _n ^{( i)}Δ t ≪ 1, the above procedure, which in general produces a rate-modulated binomial distribution, reduces to a rate-modulated Poisson process. The same time series, R _n ^{( i)}, was used to modulate the firing probabilities of each stimulated, or foreground, element contributing to the artificially generated multiunit spike train, thus producing oscillatory spatiotemporal correlations due to common input. 

Foreground vs. background

In the following, foreground refers to any pixel whose mean firing rate is above baseline, whereas background refers to any pixel whose firing rate is equal to or below baseline. Consistent with the lack of phase-locking between experimentally recorded ganglion cells responding to spatially separated stimuli, background spike trains (i.e., events arising from pixels outside the stimulus) were not subject to oscillatory modulations, except where noted in control experiments. Coherent oscillations in response to diffuse stimuli have been reported in many species, including frog (Ishikane, Gangi, Honda, & Tachibana, 2005; Ishikane et al., 1999), mudpuppy (Wachtmeister & Dowling, 1978), rabbit (Ariel, Daw, & Rader, 1983; Yokoyama, Kaneko, & Nakai, 1964), cat (Doty & Kimura, 1963; Laufer & Verzeano, 1967; Neuenschwander, Castelo-Branco, & Singer, 1999; Neuenschwander & Singer, 1996; Steinberg, 1966), monkey (Frishman et al., 2000; Laufer & Verzeano, 1967), and humans (De Carli et al., 2001; Wachtmeister, 1998). In both frog (Ishikane et al., 1999) and cat (Neuenschwander & Singer, 1996), the relative phases of such oscillations have been shown be sensitive to global stimulus topology; oscillations arising from simply connected regions remain coherent whereas oscillations arising from noncontiguous regions rapidly become uncorrelated as a function of time from stimulus onset. In control experiments, the sharp distinction between foreground and background spikes trains was relaxed, so that all model ganglion cells were subjected to oscillatory modulations whose mean amplitudes were always proportional to the local firing rate and which were either (1) uncorrelated with, or (2) in phase with, the foreground modulations. 

Representation of visual stimuli

For these experiments, visual stimuli were completely determined by specifying the mean and variance of the spike trains arising from the foreground and background pixels. Images of the corresponding input stimuli are presented solely for purposes of display. The mean stimulus intensity and, where present, the amplitude of the common oscillatory modulation were determined by adjusting the two free parameters in Equation 3, namely, the scale factor I ^{( i)} and the constant offset F ₀ ^{( i)}. These two free parameters were determined empirically so that the mean foreground and background firing rates 〈 R ^{( i)}〉 and standard deviations σ _R ^{( i)} were given by the following relations:  

$\begin{array}{c}\langle R^{\left(0\right)}\rangle =F_0^{\left(0\right)},i=0\\ \langle R^{\left(i\right)}\rangle =F_0^{\left(0\right)}(1+2^{(i-3)}),i>0,\end{array}$

(5)

 

$\begin{array}{c}{\sigma }_R^{\left(0\right)}=0,i=0\\ {\sigma }_R^{\left(i\right)}=F_0^{\left(0\right)}{2}^{(i-3)},i>0,\end{array}$

(6)

where i = {0, 1, 2, 3, 4, 5} and the different stimulus intensities are denoted by the set {0%, 25%, 50%, 100%, 200%, 400%}, with the percentages giving the change from baseline (note that not all stimulus intensities are displayed in each figure). At all background pixels, the baseline intensity was given by I ⁽⁰⁾ = 0, F ₀ ⁽⁰⁾ = 25 impulses per second (ips), so that 〈 R ⁽⁰⁾〉 = F ₀ ⁽⁰⁾, σ _R ⁽⁰⁾ = 0. For i > 0, the methodology was complicated by the fact that negative values of R _n ^{( i)}Δ t were truncated at zero, and values of R _n ^{( i)}Δ t > 1 likewise saturated at 1, making it necessary to determine I ^{( i)} and F ₀ ^{( i)} empirically via an iterative procedure. This was accomplished by explicitly calculating 〈 R ^{( i)}〉 and σ _R ^{( i)} after setting all values of R _n ^{( i)}Δ t < 0 to zero and values of R _n ^{( i)}Δ t > 1 to 1, and adjusting I ^{( i)} and F ₀ ^{( i)} so that Equations 5 and 6 were satisfied. Values of R _n ^{( i)}Δ t were then again truncated between zero and one and the process repeated until the discrepancy from the exact equality expressed by Equations 5 and 6 was less than 0.5%. 

The range of oscillatory modulations employed here produced periodic spatiotemporal correlations between co-activated cells that were consistent with electrophysiological recordings across a variety of vertebrate retinas. The multiunit activity (MUA) obtained by combining spikes across all foreground pixels exhibited periodic deflections that increased sharply in response to stimulus intensity ( Figure 1, column MUA) in a manner that was qualitatively consistent with optic nerve recordings in both cats and primates following full-field stimulation (Doty & Kimura, 1963; Laufer & Verzeano, 1967; Steinberg, 1966). Likewise, cross-correlation functions extracted from the computer-generated spike trains (Figure 1, column Cross-Correlation), obtained by averaging over all foreground pixels, were again qualitatively similar, in terms of relative peak amplitude, shape, and width, to corresponding experimental measures from cats (Neuenschwander et al., 1999; Neuenschwander & Singer, 1996) and frogs (Ishikane et al., 2005, 1999) in response to large or looming stimuli, respectively. 

Figure 1

 

Common oscillatory input. MUA: Multiunit activity. Instantaneous firing rate (1-ms bins) combining spike trains from all 16 × 16 foreground pixels (black lines) exhibit oscillatory modulations that increased with stimulus intensity, indicated on each panel as a percentage above baseline (25 Hz). When averaged over 100 independent trials, oscillatory structure in the instantaneous firing rate largely disappeared (gray lines), due to the lack of time locking to stimulus onset. Baseline firing rate indicated by dashed lines. Cross-correlation: Average pairwise cross-correlation between foreground pixels, expressed as a fraction of baseline. Both the increase in the mean firing rate, as well as the amplitude and shape of the coherent oscillatory modulations, fell within the range of published values.

Fano factors and trial-to-trial variability

Retinal ganglion cells responding to high-contrast and/or highly dynamic stimuli often exhibit greatly reduced trial-to-trail variability, as compared to a Poisson process with the same mean spike rate (Berry, Warland, & Meister, 1997; Kara et al., 2000; Passaglia & Troy, 2004; Uzzell & Chichilnisky, 2004). In a number of different retinal preparations, the Fano factor, a standard measure of the trial-to-trial variability in the number of spikes over a fixed interval (given by the mean over the variance), has been shown to decline in approximately inverse proportion to stimulus contrast (Berry et al., 1997; Kara et al., 2000; Passaglia & Troy, 2004; Uzzell & Chichilnisky, 2004). Here, an rth-order Gamma distribution, for which the Fano factor scales as 1/r, was used to model the inverse relationship between trial-to-trial variability and contrast (Stein, 1965). An instantiation of an rth-order Gamma distribution was constructed from a set of binomial distributions by choosing every rth spike (Stein, 1965). In all other respects, spikes trains based on an rth-order Gamma distribution were constructed following the same procedure as described above. 

Computer-generated spike trains II: Fixational eye movements

Overview

To assess the generality of our conclusions, computer-generated spike trains were produced implicitly via an alternative mechanism in which the amplitude, time course, and extent of the spatiotemporal correlations between pairs of ganglions cells arose entirely from the dynamics of the model, rather than being specified explicitly. Simulated fixational eye movements were superimposed on an array of model cone photoreceptors, whose output was transformed by a reduced model of negative feedback from electrically coupled horizontal cells and converted to spike trains via a piecewise linear mapping with a hard cutoff at zero. Spatiotemporal correlations between model ganglion cells were due entirely to the jitter of the visual image across the photoreceptor array. 

Fixational eye movements

The term “fixational eye movements” refers to small, involuntary changes in gaze that occur even when a subject attempts to fixate on a single point in space. Three separate mechanisms contribute to fixational eye movements, micro-saccades, drift, and tremor. Whereas drift and tremor produce very small amplitude motions, typically on the order of the width of a single photoreceptor, micro-saccades are considerably larger but occur at a relative low rate, at no more than a few per second (Martinez-Conde & Macknik, 2008). The present model addresses only small amplitude fixational eye movements due to drift and tremor that occur continuously between both saccades and micro-saccades. 

Power spectrum analysis of fixational eye movements in human subjects indicates that the drift amplitude falls off as 1/ f, where f is the temporal frequency (Eizenman et al., 1985). Drift was therefore modeled by an amplitude spectrum of the following functional form: 

$C_k^{(\mathrm{d}\mathrm{r}\mathrm{i}\mathrm{f}\mathrm{t})}=\frac{C_0^{(\mathrm{d}\mathrm{r}\mathrm{i}\mathrm{f}\mathrm{t})}}{1+f_k/f_0}\exp (2\pi i{r}_k),$

(7)

where C_k^(drift) is the Fourier amplitude of the drift at the discrete temporal frequency f_k, f₀ = 1 Hz denotes the frequency at which the spectral amplitude declines by half, C₀^(drift) = 3.33′ (min of arc) determined the average drift amplitude, and r_k is a uniform random deviate. 

Physiological nystagmus, or tremor, which is a small, high-frequency oscillation in eye position, was similarly modeled as a function of f _k by the following expression:  

$C_k^{(\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{m}\mathrm{o}\mathrm{r})}=C_0^{(\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{m}\mathrm{o}\mathrm{r})}·\exp (-\frac{(f_k-f^{(\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{m}\mathrm{o}\mathrm{r})}{)}^2}{2{\sigma }_{(\mathrm{t}\mathrm{r}\mathrm{e}\mathrm{m}\mathrm{o}\mathrm{r})}^2})\exp (2\pi i{r}_k),$

(8)

where C _k ^(tremor) is the Fourier amplitude of the tremor at the discrete temporal frequency f _k, C ₀ ^(tremor) is the peak tremor amplitude, which is equal to 6.66 s of arc, f ^(tremor) = 70 Hz denotes the central tremor frequency, and σ _(tremor) = 15 Hz is the Gaussian width of the central peak in the tremor amplitude spectrum. 

After summing the contributions from both drift and tremor, the resulting power spectrum ( Figure 2, top panel), expressed in decibels relative to the power at the discrete frequency nearest to f ₀, was qualitatively similar to published power spectra from fixating human subjects (Eizenman et al., 1985). Taking the discrete inverse Fourier transform of the combined amplitude spectrum and centering on zero yielded horizontal fixational eye movements (Figure 2, middle panel) that were also qualitatively similar in size and temporal structure to those recorded from human subjects (Eizenman et al., 1985). Vertical movements (Figure 2, bottom panel) were modeled similarly using an independent set of random deviates. 

Figure 2

 

Fixational eye movements. (Top) Amplitude spectra consisting of both 1/ f drift and 70-Hz tremor components. (Middle) Simulated horizontal fixational eye movements. Typical deflections due to tremor were of order 1′ of arc. Dotted lines indicate model cone diameter (0.75′). (Bottom) Vertical eye movements produced on the same trial.

Biophysical outer retinal model

A biophysical model (van Hateren, 2005) was used to simulate the responses of three principal neuron types in the outer retina; cone photoreceptors, bipolar cells, and horizontal cells, to rapidly changing visual stimuli, such as are likely to be produced by fixational eye movements. The biophysical cone model included rate kinetic equations for cGMP-gated channels and the light-activated breakdown and calcium-modulated synthesis of cGMP, with an instantaneous nonlinearity representing voltage-gated channels in the inner retina. Synaptic release was governed by an exponential function of membrane voltage to account for the flux of calcium into the cone terminal, with the operating voltage shifted linearly due to negative feedback from horizontal cell dendrites. Using a single set of generic parameters, the complete biophysical model of the three main outer retinal neuron types was able to qualitatively account for the amplitude and time course of measured horizontal cell responses to a wide range of impulsive, step, and sinusoidal inputs over at least a 100-fold range of background luminance values (van Hateren, 2005). The implementation used here was built upon the public domain Octave script developed by Furusawa and Kamiyama (Visiome). 

In order to model light responses across a 32 × 32 array of photoreceptors, toroidal (wrap around) boundary conditions were employed to mitigate edge effects. Electrical coupling between horizontal cell dendrites was accounted for by adding a term directly proportional to the difference between the local and surrounding horizontal cell membrane potentials, the latter obtained by computing a Gaussian-weighted local average with an e ⁻¹ fall-off radius of 4 pixels. For uniform stimuli, the additional term was identically zero (i.e., horizontal cell coupling had no effect), so that the modified model reduced exactly to the original model. For additional details regarding the biophysical outer retinal model, the reader is referred to the original report (van Hateren, 2005). 

Representation of visual stimuli

Because a primary mechanism by which fixational eye movements might produce strong spatiotemporal correlations is through the passage of a long, straight edge back and forth over an array of photoreceptors, input stimuli for these experiments consisted of vertically oriented, odd-symmetric Gabor gratings, with spatial frequency equal to 1/4 cycle per pixel (1 pixel = 1 cone diameter) and an e ⁻¹ fall-off radius equal to 8 pixels. The maximum contrast of the grating was 100% and the background luminance was equivalent to 100 trolands (as represented in the biophysical outer retinal model). Initial values were established by running the model for several seconds with only uniform background illumination. 

The dynamic variable in the biophysical model representing the bipolar cell membrane potential was employed as a surrogate for the membrane potential of the overlying model ganglion cell. Model bipolar cell membrane potentials were inverted, the mean value across all bipolar cells subtracted, and the results normalize by the maximum value across all bipolar cells (i.e., the most negative value prior to inversion). The inverted, mean-subtracted, normalized bipolar cell membrane potentials were then converted to instantaneous firing rates by multiplying by 100 Hz and adding a constant baseline of 25 Hz, corresponding to a maximum firing rate of 125 Hz to a 100% contrast grating. Spike trains were generated using an expression exactly analogous to Equation 4. A different pattern of fixational eye movements was generated for each stimulus trial. 

Estimating pairwise correlations

Method 1: SYNC

Single-trial cross-correlations were first estimated by counting the number of synchronous events between each pair of spike trains, defined as the occurrence of an event in both trains in the same 1-ms time bin. The expected number of synchronous events due to chance, computed separately on each trial, was subtracted from the actual count, yielding a correlation estimate that was positive, on average, for pairs of event generators that tended to fire together and distributed about zero for event generators that fired randomly with respect to each other. Mathematically, the SYNC-based cross-correlation, X _ij, was given by  

$X_{ij}=\sum _{t_k}\left[S_i\right(t_k)-{\bar{S}}_i]\left[S_j\right(t_k)-{\bar{S}}_j],$

(9)

where S _i denotes the spike train of the ith neuron, either 0 or 1 at the corresponding time step, t _k, and

$\bar{S}$

_i is the single-trial average. Autocorrelations were computed identically, such that each event was treated as “synchronous” with itself, thereby preserving rate-coded information along the diagonals and setting the maximum correlation amplitude. 

Method 2: γMUA*

Over very short durations, only a few events are produced by each model ganglion cell, making the estimation of pairwise correlation on single trials particularly difficult. Pairwise correlations were therefore estimated using a second method that used all pairs of spikes, not just synchronous events, by assigning a weight to each spike determined by its time of occurrence relative to the oscillatory component of the local multiunit activity (MUA). First, the MUA was estimated for each 1-ms time step by summing all events in the neighborhood of each target pixel, out to a maximum radius of 4 pixels. The contribution from each event to the MUA was weighted by the inverse of the distance to the center pixel, with distance measured as the radius of the square circumference centered on the target pixel and passing through the remote pixel. This measure of distance ensured that the contribution from each square perimeter was independent of the number of pixels at that radius, preventing the estimate of the local MUA from being dominated by the outermost rings. The local MUA was band-pass filtered using a rectangular window, typically between 60 and 100 Hz, exclusive of the endpoints, yielding the oscillatory component of the local MUA, denoted here as γMUA*, where the asterisk serves as a reminder that this quantity is merely analogous to the band-pass-filtered multiunit activity that would be measured experimentally. The γMUA*, in turn, was used to weight each event, giving positive weight to events occurring at the peaks of the γMUA* and negative weights to events occurring in the troughs. The average weight of randomly distributed events was zero, because the γMUA* itself had zero mean. 

An estimate of the oscillatory cross-correlation was given by the sum over the product of all pairs of spikes in the two trains, with each event weighted by the γMUA* evaluated at the target pixel. Mathematically, the γMUA*-based cross-correlation, Γ _ij, was given by  

${\mathrm{\Gamma }}_{ij}=\sum _{t_k,t_l}{\gamma }_i\left(t_k\right)S_i\left(t_k\right){\gamma }_i\left(t_l\right)S_j\left(t_l\right),$

(10)

where γ _i is a short-hand for γMUA* at the ith location (target pixel) and the sum is over all pairs of events, or equivalently, over all pairs of time steps, t _k and t _l, regardless of their relative timing. Note that if the multiunit measure γ _i is replaced by the single unit measure S _i, then Equation 10 reduces to Equation 9 to within an additive constant plus higher order terms. The γMUA*-based procedure for estimating spatiotemporal correlations described by Equation 10 was motivated in part by encoding schemes in which events are weighted by their time of occurrence relative to an underlying oscillation (Brody & Hopfield, 2003), only here the weights are combined into a nonlinear product. 

Autocorrelations were computed similarly, so that the diagonal terms in the pairwise correlation matrix reflected the square of the number of events produced by each spike train, with each pair of events modulated by the product of their overlaps with the local γMUA*. For Poisson distributed activity, the autocorrelations computed according to Equation 10 would typically be distributed about zero, except for the contribution of each spike train to its own γMUA*, which effectively ensured some residual autocorrelation. 

The main advantage of the γMUA*-based procedure was that it allowed correlation strengths to be estimated from only a few events based on their relative timing, regardless of whether such events were synchronous or not. Attempts to estimate oscillatory correlation strengths using conventional Fourier analysis (Kenyon, Harvey, Stephens, & Theiler, 2004), presented here as a control, were less accurate due to the extremely limited number of events available from very short spike train segments. By comparison, the γMUA* permitted the instantaneous phase of the stimulus-dependent oscillation to be reliably estimated by averaging over a small neighborhood, yielding a meaningful estimate of the correlation strength from as little as a single pair of spikes. 

Stimulus reconstruction

Rate-based reconstructions

Stimuli were first reconstructed by setting the magnitude of each pixel proportional to the logarithm of the number of events generated by the corresponding spike train, normalized relative to the expected number of spikes due to baseline activity. The maximum reconstructed pixel intensity was determined by the maximum number of events produced in response to the largest stimulus intensity across all trials. The number of trials was 100 unless otherwise noted. The minimum number of spikes was truncated at the expected baseline level, so that the minimum pixel intensity was equal to zero (i.e., log of unity). From the distributions of the number of spikes produced by stimulated vs. unstimulated event generators, an ideal threshold was determined for each stimulus intensity, allowing optimal classification of pixels as either ON or OFF (see Ideal observer classification below). Subthreshold pixel values were set to zero, which had the effect of suppressing background clutter, making the reconstructed images easier to interpret. 

SYNC- and γMUA*-based reconstructions

Stimuli were also reconstructed using the logarithm of the normalized first principal component, or largest eigenvector (eigenimage), of the pairwise correlation matrix computed using either the SYNC or γMUA* method described above. To establish an absolute scale for comparison across stimulus intensities, each eigenimage was multiplied by its corresponding eigenvalue, this being the first term in a complete orthonormal expansion. To the extent that pairwise correlations were determined solely by the common oscillatory modulation, eigenimages depended only on foreground pixels, allowing the original stimulus to be estimated via PCA. To obtain single-trial estimates that scaled linearly with stimulus intensity, the square roots of the individual pixel values were used for all correlation-based reconstructions. Pixel values were further normalized by the average reconstruction in response to baseline activity. The first principal component was obtained using the Matlab intrinsic function EIGS, with the pairwise correlation matrix, either X _ij or Γ _ij, replaced with the explicitly symmetrical construction

${\overleftrightarrow{X}}^T\overleftrightarrow{X}$

or

${\overleftrightarrow{\mathrm{\Gamma }}}^T\overleftrightarrow{\mathrm{\Gamma }}$

(the double arrow denotes a matrix) and taking the square root of the largest eigenvalue. Because only the first principal component was used in the reconstruction, the increase in condition number produced by squaring made no difference. In all runs at finite stimulus intensity, it was verified that the largest eigenvalue was at least twice the magnitude of the next largest eigenvalue, consistent with previous studies showing that the number of large eigenvalues of the pairwise correlation matrix is equal the number of distinct image segments (Kenyon, Harvey et al., 2004). To resolve an overall sign ambiguity, the average value of all stimulated pixels was required to be positive—all pixel values were multiplied by −1 if this condition was violated—except when the intensity was zero. The maximum pixel value was determined by the largest value obtained in response to the largest intensity across all stimulus trials. Background activity was again suppressed by using the distributions of stimulated and unstimulated pixel values to compute an ideal threshold and setting all subthreshold pixels to zero. An identical reconstruction procedure was used for all stimulus intensities. 

As a control, stimuli were also reconstructed using only the diagonal terms of the γMUA*-based correlation matrix, preserving information encoded in the temporal alignment of individual spike trains with the γMUA* but discarding terms that depended on the relative timing between spike pairs from different pixels. Stimuli were also reconstructed after replacing the γMUA*-based weighting of each spike by a uniform weight of one, thereby preserving spatial correlations in the total number of spikes while discarding information regarding relative spike timing or oscillatory temporal structure. 

Ideal threshold determination

An ideal discrimination threshold was determined by comparing distributions of reconstructed pixel values at different stimulus intensities. Theoretically, the maximum percentage of pixels that could be correctly discriminated is inversely related to the degree of overlap between the two distributions (Duda, Hart, & Stork, 2001). If the distributions overlap completely, the maximal theoretical performance on an intensity discrimination task would be no better than chance (50% correct). On the other hand, if the distributions were entirely nonoverlapping, the maximum theoretical performance on a discrimination task would be perfect (100% correct). Between these two extremes, corresponding to distributions that partially overlap, maximum theoretical performance on the intensity discrimination task, P, expressed as a fraction of pixels correctly classified, is given by the following formula: 

$P=(2-A_{\mathrm{o}\mathrm{v}\mathrm{e}\mathrm{r}\mathrm{l}\mathrm{a}\mathrm{p}})/2,$

(11)

where A_overlap denotes the total area of the overlap between the two distributions and the maximum value of A_overlap is normalized to one. Error bars on the estimated values of P were determined by assuming the pixel values to either side of the Bayes discriminator obeyed binomial statistics. Error bars were always negligible and therefore omitted on semilogarithmic plots. 

Results

Due to trial-to-trial variability, an independent rate code can be quite noisy, even under circumstances where the evoked activity is a factor of two or more above unstimulated baseline levels (Barlow & Levick, 1969). The intrinsically low signal/noise conveyed by an independent rate code can be illustrated by examining idealized responses to a series of simple visual stimuli presented at various intensities (Figure 3, column IMAGE). 

Figure 3

 

Single-trial reconstructions. Computer-generated spike trains, 100-ms duration, baseline rate = 25 ips, used to simulate firing activity in a 32 × 32 retinal patch. IMAGE: Stimuli were 16 × 16 uniform square spots. Intensity, indicated to the left of each row, denoted the increase in firing rate relative to baseline and, if present, the RMS amplitude of the oscillatory modulation. RATE: Representative reconstructions based on the logarithm of the number of spikes. SYNC: Reconstructions based on the logarithm of the largest principal component of the (32 × 32) × (32 × 32) correlation matrix computed from the number of synchronous events between each pair of spike trains relative to chance. γMUA*: Reconstructions computed as for SYNC but with correlations estimated by weighting each event pair based on the oscillatory component of the local multiunit activity at each target pixel ( γMUA*) and summing over all weighted event pairs. For all three reconstruction methods, the ideal threshold for classifying individual pixels as either ON or OFF was determined for each intensity level and subthreshold pixels were set to zero. Performance of optimal classifier indicated at bottom right of each panel. Fano factors indicated at bottom left. Nonlinear correlations dramatically improved stimulus reconstructions at all intensities.

Computer-generated spike trains, 100 ms in duration and derived from a binomial distribution operating close to the Poisson limit (step size, 1 ms), were used to model both background activity (25 impulses per second (ips)), as well as stimulated activity, which ranged from 25% to 400% above baseline levels (note that the intensity range plotted in most figures is less). These firing parameters encompassed both measured baseline rates (Troy & Robson, 1992) as well as much of the dynamic range of Y ganglion cells in the cat retina in response to diffuse drifting gratings presented at up to 100% luminance contrast and temporally modulated at 4 Hz (Troy & Enroth-Cugell, 1993). Although retinal neurons can fire in a highly reproducible, non-Poisson manner, especially in response to temporally complex or high-contrast stimuli (Passaglia & Troy, 2004; Reich, Victor, Knight, Ozaki, & Kaplan, 1997; Reinagel & Reid, 2000; Uzzell & Chichilnisky, 2004), a binomial distribution operating close to the Poisson limit is used here to establish a nonconfounded baseline against which the role of extreme synergy can be measured. Intrinsic reliability is addressed in separate control experiments. 

Single-trial, rate-based reconstructions were obtained by setting the intensity of each pixel to the logarithm of the number of spikes generated by the corresponding neuron, normalized relative to the average number of spikes due to baseline activity ( Figure 3, column RATE). A logarithmic scale was used to maintain sensitivity over the full 32-fold dynamic range. Background activity was suppressed by choosing a discrimination threshold that maximized the percentage of pixels correctly classified as either ON or OFF (Duda et al., 2001) and setting all subthreshold pixel values to zero. If the computed discrimination threshold was less than zero, a value of zero was used instead. The percentage of correct ON/OFF classifications, averaged across all trials, is shown at the bottom right of each reconstruction. Representative examples, here and in subsequent panels, corresponded to individual trials on which the percentage of correct classifications was closest to the multitrial average. 

Doubling the firing rate relative to baseline, from 25 ips to 50 ips over a period of 100 ms, supported performance levels of only 73% for the ON/OFF classification task, compared to an expected performance level of 50% for purely random assignments. Such poor performance can be readily explained, since doubling the firing rate produced on average only 2.5 extra spikes per pixel, corresponding to a signal/noise of only 1.6, consistent with the trial-to-trial variability exhibited by retinal ganglion cells over similar intervals in response to low-contrast stimuli (Berry et al., 1997; Kara et al., 2000; Passaglia & Troy, 2004; Uzzell & Chichilnisky, 2004). Fano factors, shown at the bottom left of each panel, declined slightly from 1.0 as a function of stimulus intensity due to deviations from a Poisson process at higher firing rates. 

Nonlinear synergistic encoding schemes, in which an oscillatory modulation was used to simultaneously vary the firing rates of all foreground spike trains, mediated substantial improvements relative to an independent rate code. To mimic the contrast dependence of the coherent oscillations between ganglion cells measured experimentally, which increase in strength with luminance contrast (Neuenschwander et al., 1999), the RMS amplitude of the oscillatory modulation, expressed as a fraction of the average response to the stimulus, was set equal to the fraction by which the mean stimulated firing rate exceeded the baseline rate (see Figure 1). 

Correlation strengths were evaluated in one of two ways. First, the number of simultaneous events between each pair of cells was summed over the 100-ms analysis interval and expressed as the difference from the expected number of synchronous events due to chance. Second, a measure that was more sensitive to periodic structure and which took account of temporal correlations at all relative delays was used. 

Briefly, an estimate of the local multiunit activity (MUA) at each target pixel was obtained by summing all spikes out to a given radius, typically 4 pixels (less for smaller retinal patches), such that the influence of each event decreased as one over the distance to the target pixel. The MUA was band-pass filtered using a rectangular window from 60 Hz and 100 Hz (or nearest discrete frequencies), yielding an estimate of the oscillatory component of the local multiunit activity, denoted as γMUA*. The γMUA* was then used to weight each spike, assigning positive weights to pairs aligned with the peaks of the local oscillation and negative weights to pairs in which one of the component spikes fell in a trough. The sum over the product of all weighted spike pairs provided an estimate of the oscillatory correlations between any two spike trains that was much more sensitive than conventional Fourier analysis when only a few events were available. 

Depending on how pairwise correlations were estimated, either using synchrony (SYNC) or the oscillatory component of the local multiunit activity ( γMUA*), stimuli were reconstructed by computing the first principal component of the corresponding pairwise correlation matrix (specifically, as the product of the largest eigenvector and eigenvalue, representing the leading term in an orthonormal expansion). Using synchrony to measure pairwise correlation strengths ( Figure 3, column SYNC), the first principal component yielded a noticeable improvement over the reconstructions mediated by an independent rate code but only at higher stimulus intensities. 

If the γMUA* was instead used to estimate pairwise correlation strengths ( Figure 3, column γMUA*), the product of the largest eigenvector and eigenvalue, or first principal component, yielded dramatic improvements in stimulus reconstruction over an independent rate code across a range of intensities spanning nearly 4 log ₂ units (16-fold). This improvement was quantified by measuring relative performance on the ON/OFF pixel classification task. Whereas a doubling of the firing rate supported only 73% correct classifications using a binomial distribution to model an independent rate code, an oscillatory correlation code, based on a rate-modulated binomial distribution, supported an average of 92% correct classifications for the same average increase in firing activity. Reliability also increased with stimulus intensity, with Fano factors falling to a minimum value of 0.38. However, separate control experiments showed that this reduction in trial-to-trial variability did not fully account for the improvements in γMUA*-based reconstructions ( Figures 10, top panel, and 12). 

Spatiotemporal correlations supported not only improved performance on the ON/OFF pixel discrimination task but also mediated improved discrimination between different stimulus intensities, here separated by factors of 2. As a function of stimulus strength, the distributions of foreground pixel values were much better separated when periodic spatiotemporal correlations were used to reconstruct the input image as opposed to the rate-based reconstructions ( Figure 4, top row). The black dotted line indicates the distribution of foreground pixel values in the absence of stimulation. Progressively lighter shades corresponded to increasing stimulus strengths, except for the highest intensity, which was denoted by a solid black line. The x-axis was scaled logarithmically to better separate the distributions at low stimulus intensities, plotted here in normalized units. Vertical dashed lines indicate the discrimination thresholds used in the ON/OFF classification task (note some thresholds were zero and thus off scale). All distributions were normalized to unity (the apparent area was distorted by the logarithmic x-scale). 

Figure 4

 

Discrimination between different stimulus intensities. Reconstruction method indicated at the top of each column. (Top) Distribution of foreground pixel values (semilogarithmic scale). Dotted line denotes the baseline distribution. Lighter shades correspond to increasing stimulus intensity (from 25% to 200%). The maximum intensity (400%) was denoted by a solid black line. Vertical dashed lines indicate discrimination thresholds used in the ON/OFF pixel class task, with shading matched to corresponding intensity. Distributions from γMUA*-based reconstructions are more widely separated than distributions for SYNC- or RATE-based reconstructions. (Middle) Percentage of correct pixel discriminations as a function of intensity difference, plotted in log ₂ units. The reference intensity, which increased from the bottom rightmost to the top leftmost curve on each plot, was matched by an upward shift in the percent correct, reflecting improved signal/noise at higher firing rates. Performance was markedly superior for the γMUA*-based reconstructions. (Bottom) Intensity discrimination threshold in log ₂ units, as a function of the percentage of correct classifications. Discrimination thresholds were 1.5 to 3 log ₂ units lower for the γMUA*-based reconstructions (curves shown for the three lowest reference intensities, which increased from top to bottom).

Distributions of foreground pixel values obtained from the rate-based reconstructions overlapped considerably, especially between lower stimulus strengths, implying poor intensity discrimination. At larger stimulus strengths, the distributions obtained from the SYNC-based reconstructions were reasonably well separated. However, only the distributions derived from γMUA*-based reconstructions were relatively well separated at nearly all intensities. 

For each pair of reconstructed images, obtained in response to different stimulus strengths, a discrimination threshold was chosen that allowed optimal classification of each foreground pixel with respect to the two intensities being compared. Performance on the intensity discrimination task, measured as the percentage of correctly classified pixels, was plotted as a function of the relative intensity difference between the two stimuli, measured in log ₂ units ( Figure 4, middle row). 

For the rate-based reconstructions, performance improved monotonically both as a function of the relative intensity difference and as the reference intensity—or weaker of the two intensities—was increased ( Figure 4, RATE, reference intensity was smallest for the bottom rightmost curve, largest for upper leftmost curve). The improved performance with increasing reference intensity reflected the fact that higher firing rates supported larger signal/noise. 

For the correlation-based reconstructions, performance on the intensity discrimination task was substantially better for most combinations of relative and reference stimulus strengths ( Figure 4, columns SYNC and γMUA*). This was especially true for the more sensitive γMUA*-based correlation measure, for which performance levels exceeded 90% for intensity differences between 1 and 3 log ₂ units, whereas an intensity difference of 3 or more log ₂ units was required to ensure the same level of performance using the rate-based reconstructions. 

Similar results were obtained when the data were reanalyzed to extract the discrimination threshold, or minimum intensity difference, required to support a given level of performance on the intensity discrimination task ( Figure 4, bottom row). Compared to the rate-based reconstructions, discrimination thresholds were uniformly lower for the correlation-based reconstructions. At the 90% performance level, discrimination thresholds for the γMUA*-based reconstructions ranged from 1 to 2.5 log ₂ units lower than for the rate-based reconstructions. Thus, on a pixel-by-pixel basis, spatiotemporal correlations permitted greatly superior discrimination between different stimulus intensity levels. 

The dramatic improvements in signal/noise mediated by spatiotemporal correlations presumably reflected the aggregate contribution from hundreds of cells, or tens of thousands of stimulated cell pairs. It follows that smaller patches, containing smaller numbers of spike train pairs, should support proportionately smaller improvements in signal/noise compared to an independent rate code. 

To test this prediction, reconstructions were performed using progressively smaller retinal patches ( Figure 5). The percentage of correct classifications on the ON/OFF pixel discrimination task was plotted as a function of the total size of the retinal patch for the 4 largest stimulus intensities employed. For the rate-based reconstructions (solid line), the percentage of correct classifications was insensitive to the size of the retinal patch, as expected for independent encoders. 

Figure 5

 

Dependence on retinal patch size. Performance on the ON/OFF pixel classification task was plotted as a function retinal patch diameter, with the correlation matrix constructed from all cell pairs. Stimulus intensity indicated on each panel. The precision of the rate-based reconstructions (solid lines) was insensitive to patch size. At small to moderate stimulus intensities, performance mediated by the γMUA*-based reconstructions (dashed lines) rose steeply as a function of the total number of cells, saturating for patch sizes of approximately 24 × 24. At the highest intensity, superior γMUA*-based reconstructions required only a few cells. The quality of the SYNC-based reconstructions (dotted lines) exhibited a more complex dependence on patch size, due to noise in the underlying estimates of pairwise correlation strength, which was the dominant factor at low stimulus intensities. Number of trials: {100, 100, 200, 300, 400, 600, 800, 1000} for patch diameters {32, 28, 24, 20, 16, 12, 8, 4}, respectively.

For the γMUA*-based reconstructions (dashed line), performance increased monotonically with patch size, saturating at approximately 24 × 24 total pixels, or at 12 × 12 stimulated pixels. This latter number is in reasonable agreement with estimates of the size of redundant cell neighborhoods in the salamander retina (Puchalla et al., 2005). 

Paradoxically, for SYNC-based reconstructions at lower stimulus intensities, increasing patch size led to lower performance on the ON/OFF pixel classification task, presumably because the first principal component was dominated by noise in the off-diagonal terms of the pairwise correlation matrix, swamping rate-coded information distributed along the diagonal. On the other hand, even for very small patches, the performance mediated by the γMUA*-based reconstructions was either equivalent to, or for higher stimulus intensities measurably better than, the performance mediated by the rate-based reconstructions. 

Overall, these results support the hypothesis that for weak to moderate pairwise correlation strengths, the improvements in signal/noise provided by the γMUA*-based reconstructions required simultaneous processing of hundreds of spike trains, corresponding to hundreds of thousands of pairwise correlations. 

Remarkably, the above nonlinear reconstruction methods do not eliminate fine spatial detail, as each term in the pairwise correlation matrix, which is based on a sum over the product of event pairs, was dependent on local firing activity. The preservation of substantial fine spatial detail was demonstrated implicitly in the ability to reproduce sharp edges in the correlation-based reconstructions of uniform square spots (e.g., Figure 3). When individual pixels were randomly deleted from the stimulus, thus creating an arbitrary complex shape, γMUA*-based reconstructions continued to preserve substantial spatial structure at the level of individual pixels ( Figure 6). 

Figure 6

 

Preservation of substantial fine spatial detail. Explanation of panels as in Figure 3 with column SYNC omitted. Stimuli were again 16 × 16 square spots but with 10% of the pixels randomly deleted. The percentage of correct ON/OFF pixel classifications was essentially identical as for uniform square spots. The largest principal component of the pairwise correlation matrix mediated greatly improved signal/noise without eliminating fine spatial detail. Similar improvements in signal/noise via conventional spatial averaging using predefined templates would have required precise foreknowledge of which pixels had been deleted.

The reconstruction of images with randomly deleted pixels is not physiologically realistic, as the present model ignores lateral inhibition, which would normally suppress the firing of an isolated pixel surrounded by stimulated pixels. Nonetheless, the analysis of a random pattern of deletions demonstrates how nonlinear spatiotemporal correlations allow information to be distributed across many pixels while maintaining a high fidelity representation of the pixel-by-pixel intensity. The loss of some spatial detail is revealed by the fact that there were very few mistakes among either ON foreground pixels or OFF background pixels outside the central square, whereas mistakes among deleted pixels were much more common. Nonetheless, to obtain comparable improvements in signal/noise via conventional spatial averaging—without eliminating fine spatial detail—would have required a separate template for every combination of deleted pixels. Reconstruction based on γMUA*, on the other hand, required only the information available on single trials from the spike trains themselves. The present results demonstrate that even over very short time scales, from tens to hundreds of milliseconds, enough information may be available in the pairwise correlation matrix, especially when referenced to the local multiunit activity (i.e., γMUA*), to construct nonlinear spatial filters on the fly, thereby conferring the main advantages of spatial averaging (i.e., improved signal/noise) without eliminating fine spatial detail. 

Both psychophysical and electrophysiological data suggest that 50 to 300 ms are required to process visual scenes (Bacon-Mace et al., 2005; Kirchner & Thorpe, 2006; Rolls et al., 1999; Thorpe et al., 1996). The lowest estimates reflect the latency of image-specific responses recorded from individual neurons in IT compared to latencies recorded in V1 (Rolls et al., 1999). The longest estimates are based on the upper range of mean saccade reaction times for an image classification task in which only correct trials are considered (Kirchner & Thorpe, 2006) and falls within the range of intersaccade intervals measured in humans and other primates (Martinez-Conde, Macknik, & Hubel, 2004). 

Consistent with the observed rapidity of visual processing, oscillatory spatiotemporal correlations were found to support improved stimulus reconstructions over time scales as short as 25 ms ( Figure 7). Superior reconstructions were obtained from 25-ms spike trains even though the individual frequency components used in estimating the γMUA* at each target cell were quantized at 40-Hz intervals, so that periodic structure could at best be only poorly resolved. Correlation-based reconstructions were particularly accurate at higher intensities, with the first principal component becoming nearly perfect as the stimulus strength approached maximum values. 

Figure 7

 

Reconstructions from 25-ms spike trains. Explanation of panels as in Figure 3. The γMUA*-based reconstructions were again greatly superior to those mediated by an independent rate code, despite the necessarily low resolution of the underlying frequency components. SYNC-based reconstructions were also superior at very high stimulus intensities.

Retinal ganglion cells typically fire most strongly during the transient phase of the response immediately following stimulus onset, settling to a lower plateau level of activity after several tens of milliseconds (Cleland, Levick, & Sanderson, 1973; Creutzfeldt, Sakmann, Scheich, & Korn, 1970). While the largest intensities modeled here are consistent with the highest levels of activity measured in cat Y ganglion cells in response to diffuse, sinusoidally modulated gratings (Troy & Enroth-Cugell, 1993), abrupt light steps produce transient peaks of activity that equal or exceed the maximum response amplitudes considered here (Cleland et al., 1973; Creutzfeldt et al., 1970). When assessing the quality of rate- vs. correlation-based reconstructions extracted from very short time spike train segments—on the order of 25 ms—the results derived here from higher stimulus intensities may therefore be more relevant, with the initial transient response component potentially serving the same function as performed by coherent oscillations in the present study. 

To characterize how the quality of the rate-, SYNC- and γMUA*-based reconstructions depended on the length of the analysis window, optimal theoretical performance on the ON/OFF pixel classification task was plotted as a function of spike train duration, ranging from 25 to 400 ms ( Figure 8). As expected, performance for all three reconstruction methods declined as the length of the analysis window was decreased, yet even for the shortest spike trains tested, the γMUA*-based reconstructions mediated substantially superior ON/OFF pixel discrimination across a nearly 16-fold range of intensities. 

Figure 8

 

Reconstructions vs. spike train duration. Reconstruction quality assessed as the percentage of pixels correctly classified as either ON or OFF using an optimal discrimination threshold. Stimulus intensity indicated on each panel. Note logarithmic time scale. Rate-based reconstructions (solid lines) took substantially longer than γMUA*- (dashed lines) and, at high intensities, SYNC-based reconstructions (dotted lines), to achieve comparable accuracy.

At a stimulus intensity of 100%, representing a doubling in the firing rate and an RMS oscillatory amplitude that was two times baseline, γMUA*-based reconstructions supported performance levels exceeding 90% in less than 100 ms, whereas nearly 400 ms were required to achieve the same level of performance using rate-based reconstruction methods. Thus, γMUA*-based reconstructions required approximately 1/4 the number of events to achieve accuracy comparable to an independent rate code. Given that the signal/noise in response to slowly varying stimuli presented at low to moderate contrast generally improves as the square root of the number of events (Barlow & Levick, 1969), the present findings imply that utilizing spatiotemporal correlations approximately doubles the potential signal/noise, an improvement comparable to quadrupling the number of spikes. 

Pairwise spatiotemporal correlations could support reasonably accurate reconstructions even when the firing rate was unchanged from baseline ( Figure 9). Although signal/noise was degraded in the absence of rate-coded signals, strong coherent oscillations alone could support reasonably accurate reconstructions, allowing the intensity of a large, contiguous region to be discriminated even at locations where the firing rate remained at baseline levels. This control demonstrates that pairwise correlations are capable of encoding precise spatial information that is independent of, and thus can augment, rate-coded signals. 

Figure 9

 

Reconstructions in the absence of rate-coded information. Explanation of panels as in Figure 6. Firing rates were held fixed at baseline levels. Stimulus intensities indicate the RMS amplitude of the oscillatory modulation relative to baseline. Rate-based reconstructions were indistinguishable from chance. Oscillatory spatiotemporal correlations supported levels of accuracy in the range of 80% to 90% at higher stimulus intensities.

To parse the relative factors underlying the ability of oscillatory spatiotemporal correlations to mediate superior reconstructions, a series of experiments was performed to isolate the contributions from various sources of information, such as purely temporal modulations, spatial correlations, changes in the mean firing rate, and the common reference signal provided by the γMUA* ( Figure 10). Reconstruction quality was assessed as optimum performance on the ON/OFF pixel classification task, plotted as a function of stimulus intensity. 

Figure 10

 

Relative contributions of spatial and temporal factors. Reconstruction quality plotted as optimal ON/OFF pixel classification performance vs. stimulus intensity. Solid and dashed black lines in all three panels give performance of the previously described rate- and γMUA*-based reconstructions, respectively. (Top) Rate-based reconstructions using spike trains subjected to γ-band temporal modulations (gray-solid line) were only slightly improved relative to Poisson-like event trains, despite reduced trial-to-trial variability. Alternatively, γMUA*-based reconstructions were noticeably degraded in the absence of proportionate firing rate increases (gray-dashed line), although spatiotemporal correlations alone supported performance similar to that mediated by a pure rate code. (Middle) Incorporating spatial correlations in the total number of spikes into rate-based reconstructions, but ignoring temporal correlations, yielded no improvement (gray-solid line), as spatial fluctuations in spike counts were not stimulus related. Ignoring off-diagonal terms in the pairwise correlation matrix degraded γMUA*-based reconstructions (gray-dashed line), highlighting the importance of spatiotemporal correlations. (Bottom) When oscillatory pairwise correlations were assessed using standard Fourier analysis, the corresponding reconstructions were degraded (gray-dashed line), underlining the relative efficiency of the γMUA*-based procedure. The performance mediated by SYNC-based reconstructions is shown for comparison (gray-dotted line).

At low levels of activity, Fano factors measured from cat retinal ganglion cells approach 1.0, consistent with a Poisson process, but in response to optimal, high-contrast stimuli, Fano factors are much smaller than 1.0, ranging from 0.05 to 0.52 with a mean value of 0.15 (Kara et al., 2000). The increased reliability of strongly stimulated retinal ganglion cells was partially mimicked in the present study by the reduction in trial-to-trial variability produced by large amplitude oscillatory modulations, which reduced the average Fano factors from a value of 1.0 at baseline firing rates to a value of 0.38 ± 0.05 at the highest stimulus intensity. For a stationary Poisson process, the Fano factor always equals 1.0 regardless of mean firing rate. However, spike trains in the present study are derived from binomial distributions, which only reduce to rate-matched Poisson processes in the limit that the probability for an event in any given bin is much less than one. For very strong oscillatory modulations, this assumption was violated. 

To isolate the contribution from reduced variability, a comparison was made between (1) rate-based reconstructions that used temporally modulated spike trains, generated via the same mathematical process used to produce γ-band oscillations ( Figure 10, top panel, gray-solid line) and (2) rate-based reconstructions using nonmodulated spike trains (black-solid line). Oscillatory temporal structure, by reducing the variability in the number of spikes generated over the 100-ms trial, yielded only a small improvement in signal/noise. Although retinal spike trains can exhibit significant temporal structure (Reich et al., 1997; Rodieck, 1967; Troy, Schweitzer-Tong, & Enroth-Cugell, 1995), these findings indicate that the oscillatory temporal structure employed here cannot, by itself, account for the superior quality of the γMUA*-based reconstructions. 

The mean firing rate did contribute substantially to the quality of γMUA*-based reconstructions ( Figure 10, top panel, dashed gray line). Although γMUA*-based reconstructions were adversely affected by the loss of rate-coded information, oscillatory spatiotemporal correlations—even with the firing rate held at baseline—supported performance levels roughly equivalent to those mediated by a stationary rate code (gray-dashed vs. black-solid lines). 

A related question is whether purely spatial correlations, independent of the contribution from oscillatory temporal structure, could have contributed meaningfully to the superior quality of the γMUA*-based reconstructions. A new correlation matrix was computed, which ignored precise temporal structure while preserving spatial correlations in the total number of spikes. Foreground pixels were still subjected to oscillatory modulations but the recomputed correlation matrix discarded all dependence on relative spike timing ( Figure 10, middle panel, gray-solid vs. black-solid lines). 

The reconstructions based on purely spatial correlations were nearly identical to the rate-based reconstructions. Thus, spatial correlations in the number of events, even between cells subjected to a common oscillatory modulation, did not contribute meaningfully to the superior quality of the γMUA*-based reconstructions. This finding is reasonable given that spatial correlations in the total number of events were not explicitly tied to the stimulus intensity. 

To assess the contribution from precise temporal structure alone, modified γMUA*-based reconstructions were obtained using only the diagonal terms in the pairwise correlation matrix, thus preserving the autocorrelation of each spike train with the local phase of the multiunit activity but discarding all cross-correlations. Such diagonal-only constructions were markedly inferior to the standard γMUA*-based reconstructions that utilized off-diagonal correlations as well ( Figure 10, middle panel, gray-dashed vs. black-dashed lines), except at high stimulus intensities where differences may have been masked by saturation effects. Spatiotemporal correlations thus make a critical contribution to the superior quality of γMUA*-based reconstructions, but only when both spatial correlations and precise temporal structure are simultaneously taken into account. 

Experiments were also conducted to evaluate the importance of the γMUA*-based reconstruction technique, which allowed pairwise correlations to be estimated from only a limited number of events. A new pairwise correlation matrix was computed by performing a conventional Fourier analysis of each spike train segment. Correlation strengths were estimated by the peak in the cross-power spectra between 60 Hz and 100 Hz (given by product of the Fourier amplitudes from the two spike trains) multiplied by the cosine of the difference in relative phase at the peak frequency. The dependence on the relative phase ensured that synchronous oscillatory activity corresponded to a positive correlation, whereas activity that was out of phase corresponded to a negative pairwise correlation. Taking the average of all such products over the interval 60 Hz to 100 Hz, roughly corresponding to the area under the peak, yielded equivalent results (data not shown). Reconstructions based on a pairwise cross-power spectra analysis were greatly inferior to the γMUA*-based reconstructions ( Figure 10, bottom panel, gray-dashed vs. black-dashed lines) but were slightly better than SYNC-based reconstructions (gray-dotted curve). 

The poor reconstructions derived from a conventional Fourier analysis reflect the difficulty of estimating oscillatory correlations from a very small number of events. The γMUA*, by averaging over a local neighborhood containing multiple cells, yielded an estimate of the instantaneous phase of the common oscillation that allowed a meaningful weight to be assigned to every interspike interval. Thus, by taking into account temporal structure linking cells responding to the same stimulus, the oscillatory component of the local multiunit activity permitted an estimate of the correlation strength to be obtained for each pair of events. 

So far, it has been assumed that the background pixels remained completely free of oscillatory modulations regardless of the strength of the foreground modulation, an assumption that is physiologically untenable. In practice, especially for high-contrast stimuli, foreground/background pixels would be separated by an inhibited buffer region. Nonetheless, the physiological mechanisms responsible for coherent high-frequency retinal oscillations among contiguously activated foreground pixels would likely have some effect on nearby background pixels as well, potentially producing strong spatiotemporal correlations both within and between foreground and background regions. To investigate whether oscillatory modulations of the background firing rate would obviate the proposed synergistic encoding scheme, event generators corresponding to background pixels were subject to oscillatory modulations that were of the same central frequency and relative magnitude as the foreground modulations (i.e., at a stimulus intensity of 100%, the RMS amplitude of the oscillatory modulations applied to the background firing rates was equal to 25 Hz; Figure 11). 

Figure 11

 

Effect of background oscillations. Columns 1–3: Input stimuli, rate- and γMUA*-based reconstructions for unmodulated background firing rates of 25 Hz (panels replotted from Figure 3 for ease of comparison). Columns 4–5: Background pixels subject to oscillatory modulations of the same relative magnitude as the foreground modulations using either uncorrelated ( γMUA*-I) or identical ( γMUA*-II) waveforms. The proposed synergistic encoding scheme was robust with respect to relatively strong background oscillations except at low stimulus intensities.

Whether the background modulations were governed by a waveform that was statistically independent of the foreground modulations ( γMUA*-I) or whether the same waveform (rescaled) was used for both regions ( γMUA*-II), reconstructions based on oscillatory spatiotemporal correlations were superior to rate-based reconstructions, except at low stimulus intensities. At higher stimulus intensities, the difference in pairwise correlation strength due to the difference in absolute modulation amplitude was sufficient to permit reasonable discrimination between foreground and background pixels, despite the presence of relatively strong oscillatory modulations across the entire image, even when such modulations were in phase with the foreground modulations. 

Although strong oscillatory modulations by themselves can substantially reduce trial-to-trial variability, here lowering Fano factors from 1.0 to 0.38, there remains the question of how intrinsic, stimulus-related increases in the reliability of retinal output, independent of increases in pairwise correlation strength, might affect the performance of a rate- vs. correlation-based reconstruction. Moreover, Fano factors measured from cat retinal ganglion cells in responses to high-contrast, drifting sinusoidal gratings have an average value of 0.15 and range as low as 0.05 for some neurons, which is near the theoretically derived minimum (Kara et al., 2000). Thus, there remains the related question of what additional pixel-by-pixel intensity information, if any, might be conveyed by extreme synergy when the underlying spike trains span the full range of measured Fano factors. To addresses these questions, an rth-order Gamma distribution, in which the Fano factors scale as 1/r (Stein, 1965), was used to generate artificial spike trains whose trial-to-trial variability could be specified independently of pairwise correlation strength (Figure 12). 

Figure 12

 

Effect of intrinsic reliability. RATE: Rate-based reconstruction derived from a stationary binomial distribution; Fano factors (indicated at bottom left of each panel) decline slightly with intensity due to deviations from a Poisson process. RATE-II: Rate-based reconstruction derived from a stationary rth-order Gamma distribution; from top to bottom, Fano factors scale as 1/ r, r = {3, 4, 5, 6}. RATE-III: Rate-based reconstruction derived from rth-order Gamma distribution modulated by a common oscillatory waveform; Fano factors approach theoretical minimum values at the highest stimulus intensity. γMUA*-III: Correlation-based reconstructions obtained from the same spike trains used in previous column. Despite increased intrinsic reliability, oscillatory pairwise correlations supported improved stimulus reconstructions at all intensities.

Based on a several studies showing that Fano factors from retinal ganglions fall in approximately inverse proportion to stimulus contrast (Berry et al., 1997; Kara et al., 2000; Passaglia & Troy, 2004; Uzzell & Chichilnisky, 2004), r-values were set equal to the superscript i denoting stimulus intensity in Equations 3–6, yielding Fano factors that ranged from 1.0 down to 0.18. Compared to spike trains derived from a stationary binomial distribution, reductions in trial-to-trial variability resulting from using an rth-order Gamma distribution yielded quantitative improvements in rate-based reconstructions at all stimulus intensities tested (Figure 12, RATE-II vs. RATE). When a common oscillatory modulation was superimposed on an rth-order Gamma distribution, further improvements in rate-based reconstructions were again observed (Figure 12, RATE-III vs. RATE-I). Here, an rth-order, oscillatory Gamma distribution was constructed from a rate-modulated binomial process by taking every rth spike, after scaling the initial instantaneous firing rate by a factor of r. Despite improvements in rate-based reconstructions due to increased reliability of the underlying spike trains, reconstructions based on oscillatory pairwise correlations remained superior at all stimulus intensities (Figure 12, γMUA*-III). The superiority of the γMUA*-III-based reconstructions demonstrates that intrinsic reliability can work in concert with extensive pairwise correlations to yield greater reductions in trial-to-trial variability than would result from either mechanism alone. 

The above results illustrate how, for a fixed level of trial-to-trial variability, stimuli incident on the retina can be encoded using fewer spikes by distributing local, pixel-by-pixel information across hundreds of neurons in the form of oscillatory spatiotemporal correlations. Specifically, the above example of an extremely synergistic code based on synchronous oscillations allowed pixel-by-pixel information to be read out over physiologically relevant time scales (i.e., 25–400 ms) employing only the spike trains themselves, in the absence of any foreknowledge of the input stimulus. 

The above results, however, were based on an explicit mathematical model in which spatiotemporal correlations were specified a priori. Although specifying the spatiotemporal correlations explicitly permitted a more detailed analysis, it also reduces the generality of the results. To investigate whether the above findings were dependent on the details of a particular set of mathematical assumptions, a final set of experiments was conducted in which spatiotemporal correlations were produced implicitly through the action of a completely different yet well-established physiological mechanism, namely, the rapid translation of images across the retina due to fixational eye movements. 

Here, intersaccadic fixational eye movements, which excluded micro-saccades as well as large amplitude saccades, were modeled using simplified mathematical expressions to approximate the frequency spectrum due to drift and tremor recorded from fixating human subjects (Eizenman et al., 1985). A stimulus consisting of a vertically oriented Gabor grating, with a spatial frequency of approximately 24 cycles/degree—well below the Nyquest frequency for central vision (Hirsch & Curcio, 1989; Williams, 1986), was presented to a model of the primate outer retina (van Hateren, 2005), modified to incorporate short-range electrical coupling between horizontal cell dendrites. Assuming zero divergence in the foveal P-pathway (Calkins et al., 1994), bipolar cell membrane potentials in the biophysical model were used as surrogates for midget ganglion cell firing probabilities, after a linear transformation to match baseline and maximum firing rates of 25 and 125 ips, respectively. In order to reproduce the natural jitter of the retinal image due to drift and tremor, the input stimulus was simultaneously translated in both the vertical and horizontal directions in a manner that approximated the amplitude and frequency of human intersaccadic fixational eye movements. Assuming a conversion factor for human retina of 291 μm/degree and a foveal cone center-to-center spacing of 3.0 μm (Hirsch & Curcio, 1989; Williams, 1986), the full 32 × 32 model cone array corresponded to a patch of central vision approximately 1/3 degree on a side. 

Spatiotemporal correlations produced entirely by intersaccadic fixational eye movements, as assessed by the number of synchronous events above baseline measured on single trials, either 128 ms or 256 ms in duration, yielded substantial improvements in the quality of reconstructed images as compared to an independent rate code ( Figure 13). Image regions with positive contrast were labeled foreground pixels, while image regions with negative contrast were labeled background. Rate-coded reconstructions were degraded due to the trial-to-trial variability in the number of spikes generated at each pixel and by the cumulative blurring that resulted from ongoing fixational eye movements. Reconstructions based on pairwise synchrony, on the other hand, were generally consistent with the original image and improved with increased recording duration over the range of intersaccadic intervals tested. Whereas simply averaging the number of spikes from each pixel must invariably be confounded by retinal image motion, detailed spatial information can be preserved in the form of nonlinear pairwise correlations. Reconstructions based on the γMUA* were relatively poor, as there was little narrow-band oscillatory energy in the model spike trains (results not shown). Performance on the ON/OFF pixel discrimination task was effectively at chance for the rate-based reconstructions; in part due to the use of graded as opposed to binary images. Performance on the same ON/OFF pixel discrimination task was substantially better using the reconstructions based on pairwise synchrony, which reflected the preservation of fine spatial detail. Previously, it had been demonstrated that fixational eye movements can improve visual perception (Rucci, Iovin, Poletti, & Santini, 2007b) and associated modeling results suggest a critical role for spatiotemporal correlations between retinal neurons (Poletti & Rucci, 2008). Consistent with previous studies, the present findings suggest that the spatiotemporal correlations resulting from intersaccadic fixational eye movement encode useful pixel-by-pixel information about visual stimuli not present in the time-averaged response at each location. 

Figure 13

 

Fixational eye movements. Spatiotemporal correlations were generated by simulating the small amplitude, rapid shifts in retinal image location due to ocular drift and tremor. Transformed bipolar cell membrane potentials derived from a modified biophysical model of the outer primate retina were employed as surrogates for ganglion cell firing probabilities (baseline and max firing rates, 25 and 125 ips, respectively). IMAGE: Vertically aligned Gabor gratings (100% contrast, spatial frequency ∼24 cycles/degree) presented for (top) 128 ms and (bottom) 256 ms. Red traces: Examples of fixational eye movements of indicated duration. RATE: Rate-coded reconstructions were degraded, due to variability in the number of spikes and to cumulative blurring that resulted from fixational eye movements. SYNC: Reconstructions based on pairwise synchrony were not degraded by image motion, but rather improved with longer spike train durations. Performance on the ON/OFF pixel discrimination task (percent correct shown at bottom right of each panel) was made more difficult by the use of graded as opposed to binary images, yet the spatiotemporal correlations due to fixational eye movements yielded clear improvements relative to an independent rate code.

Discussion

As independent encoders, retinal ganglion cells appear poorly designed to convey precise pixel-by-pixel intensities on short time scales. In response to sustained or smoothly varying stimulation at low to moderate contrast, retinal spike trains can be quite noisy, exhibiting considerable trial-to-trial variability (Barlow & Levick, 1969) with Fano factors close to 1.0 at low luminance contrast (Berry et al., 1997; Kara et al., 2000; Passaglia & Troy, 2004; Uzzell & Chichilnisky, 2004). The problem of extracting information from retinal spike trains is especially critical given that subjects require only 50–300 ms to process visual scenes (Bacon-Mace et al., 2005; Hung, Kreiman, Poggio, & DiCarlo, 2005; Kirchner & Thorpe, 2006; Rolls et al., 1999; Thorpe et al., 1996), a time interval over which the standard deviation in the number of spikes generated by a typical retinal ganglion cell is often a substantial fraction of the mean (Barlow & Levick, 1969). How then, in the face of such pronounced trial-to-trial variability, can the perceived acuity of visual experience be accounted for? 

The findings presented here suggest an explanation based on extreme synergy, in which the messages conveyed by individual retinal ganglion cells can only be properly interpreted within the context of the concurrent firing activity across extended neighborhoods encompassing hundreds of cells, roughly commensurate in size to the classical antagonistic surround. Using computer-generated spike trains consistent with recordings from retinal ganglion cells, the present findings demonstrate that spatiotemporal correlations can support pixel-by-pixel estimates of the input stimulus intensity that in many cases are dramatically superior to those mediated by an independent rate code, even when the reduction in trial-to-trial variability in the total spike count produced at high luminance contrast is taken into account. The proposed synergistic encoding mechanism is intrinsically nonlinear and thus represents an alternative to reconstruction strategies based on optimal linear filters (Stanley et al., 1999; Warland et al., 1997), in which the contribution from each spike is incorporated independently of any other spikes. 

Here, biologically plausible, extremely synergistic spike trains were constructed in which neurons conveyed information through multiple channels; in their mean firing rates, as temporally modulated firing patterns, and via pairwise spatiotemporal correlations. The present findings illustrated how these additional information channels could be exploited to greatly abet the process of reconstructing retinal images from parallel spikes trains, using fewer events than would otherwise be required if only rate-coded signals were considered. As a general decoding strategy, Principal Components Analysis (PCA) was used to extract the largest eigenimage, defined as the image that in an orthogonal expansion best accounted for the measured pairwise correlations. The pairwise correlation matrices considered here included spike-count information along the diagonals, so that in the degenerate case in which all off-diagonal correlations were negligible, PCA was no worse than a rate-based reconstruction. By design, there was only one principal component (one large eigenvalue), which was proportional to the vector of foreground pixels. In general, the same PCA analysis can be used to extract multiple eigenimages, each corresponding to a different principal component of a visual scene (Kenyon, Harvey et al., 2004). 

Pairwise correlations were first estimated by counting the number of synchronous events between each pair of cells over the recording interval. When the expected number of synchronous events was too small, however, estimates of pairwise correlations were very noisy, severely compromising the resulting reconstructions. Using longer time windows to define synchronous events did not improve the correlation estimates, as signal/noise was always largest at zero lag. Alternatively, when modulated by a common oscillatory input, single-trial estimates of the pairwise correlation matrix could be greatly improved by weighting each pair of spikes by their time of occurrence relative to the oscillatory multiunit activity of the surrounding cells. This improvement resulted from the fact that spike pairs whose relative timing was aligned with the oscillatory component of the local multiunit activity at the target pixel were more likely to represent the stimulus, whereas nonaligned spikes were more likely to reflect background activity or noise. In the extreme (unphysiological) limit in which the retina produced only synchronous foreground spikes and randomly distributed background spikes, the optimal decoding strategy would be to estimate local stimulus intensities using only synchronous spikes. Although the present reconstruction algorithm, based on PCA, did not involve such a clean division between foreground and background spikes, the underlying principle is analogous. 

The extreme synergy hypothesis was motivated in large part by extrapolations of marginal entropy measurements to large ensembles of weakly correlated ganglion cells stimulated by movie clips, which indicate that the information generated by the responses of individual retinal neurons is almost wholly accounted for by the surrounding activity (Schneidman, Berry, Segev, & Bialek, 2006). Marginal entropy measurements in simultaneously recorded mosaics of primate parasol ganglion cells do not exhibit similar scaling with neighborhood size in response to either uniform gray backgrounds or checker board white noise (Shlens et al., 2009), consistent with the assertion that large, contiguous stimuli are necessary for producing long-range correlations over distances beyond nearest neighbors. Here, the use of PCA permitted the information distributed across large neural ensembles to be recombined at each pixel, rapidly accomplishing significant improvements in signal/noise without sacrificing fine spatial detail. 

Previous attempts to measure synergy between pairs of neurons in the visual pathway have obtained contradictory results (Dan et al., 1998; Gawne & Richmond, 1993; Hirabayashi & Miyashita, 2005; Nirenberg et al., 2001; Puchalla et al., 2005; Reich et al., 2001; Rolls et al., 2003; Samonds et al., 2004; Singer, 1999). No previous experimental study, however, has directly considered the precise pixel-by-pixel information encoded across hundreds of thousands of correlated cell pairs (corresponding to hundreds of individual cells) recorded simultaneously; leaving open the possibility that extreme synergy may only be apparent at sufficiently large spatial scales. Consistent with measurements of synergy between pairs of ganglion cells responding to natural images (Nirenberg et al., 2001; Puchalla et al., 2005), the present theoretical analysis found relatively little synergy across small retinal patches containing only a few ganglion cells, unless the pairwise correlations were both very strong and explicitly stimulus-dependent. Rather, it was only when patch size increased sufficiently so that large numbers of pairwise correlations could be processed simultaneously that significant amounts of additional information became available. 

Other experimental studies have shown that pairwise correlations among large numbers of retinal ganglion cells do not improve estimates of global stimulus parameters based on population-coded averages (Frechette et al., 2005). The present theoretical study, however, relates only to the precise pixel-by-pixel information distributed across simultaneously recorded neurons. When estimating the global stimulus properties of a large object—such as speed—using a predefined spatiotemporal filter, the improvements in signal/noise obtained by averaging over many cells are so substantial that pairwise spatiotemporal correlations are unlikely to yield detectable improvements—and indeed may limit the advantages of spatial averaging (Shadlen & Newsome, 1994), although the stimulus information encoded by the correlations themselves can compensate for this effect (Kenyon, Theiler et al., 2004). However, to average over cells in a manner that preserves fine spatial detail requires prior knowledge of the precise dimensions of the stimulus, whereas PCA required only the pairwise correlation matrix extracted on single trials from the spike trains themselves without any prior knowledge of the input image. 

It has also been proposed that rapid, pixel-by-pixel intensity information is transmitted to the brain via spike latency codes (Van Rullen & Thorpe, 2001), and reconstructions based on the timing of first spikes measured relative to the population response exhibit similar improvements—compared to a rate code—as illustrated here (Gollisch & Meister, 2008). In principle, latency codes could be multiplexed with long-rage spatiotemporal correlations, as correlations per se do not impose any particular order of firing within any given cycle. Similarly, spatiotemporal correlations are not inconsistent with evidence that synchronous spikes can encode higher spatial resolution (Schnitzer & Meister, 2003), as the short-range synchrony postulated to underlie spatial hyperacuity may involve different, although possibly overlapping, synaptic mechanisms (Meister & Berry, 1999) to those responsible for long-range synchrony (Kenyon et al., 2003). Finally, optimal predefined filters, especially those that take account of overlapping ganglion cell surrounds (Stanley et al., 1999; Warland et al., 1997), would likely provide additional intensity information to the explicitly nonlinear, synergistic encoding scheme examined here. Although extreme synergy could not produce reliable reconstructions at lower stimulus intensities, there remains the intriguing possibility that a combination of retinal mechanisms, including fixational eye movements, long-range oscillations, latency codes, and specialized filters, could collectively explain the vividness of visual perception even at low contrast. 

The present findings utilized two entirely different mechanisms for generating realistic spatiotemporal correlations. The primary mechanism considered here, coherent oscillations between retinal neurons, has been observed in many vertebrate species (Ariel et al., 1983; De Carli et al., 2001; Doty & Kimura, 1963; Ishikane et al., 1999; Laufer & Verzeano, 1967; Neuenschwander et al., 1999; Steinberg, 1966; Wachtmeister & Dowling, 1978; Yokoyama et al., 1964) and are potentially of clinical relevance (Adams & Dawson, 1971; Fortune, Bearse, Cioffi, & Johnson, 2002; Frishman et al., 2000; Rangaswamy, Hood, & Frishman, 2003; Wachtmeister, 1998). Experimental observation of coherent oscillations among retinal neurons may be confounded, however, by the lack of strong phase-locking to the stimulus onset and the possibility of wash-out in multitrial averages due to random fluctuations in the underlying Fourier components (Stephens et al., 2006). Species differences, adaptation state, and the level and type of anesthesia have also been reported to affect the amplitude of retinal oscillations (Doty & Kimura, 1963; Steinberg, 1966) and thus may confound experimental observations. Moreover, high-frequency retinal oscillations are strongly size dependent (Ishikane et al., 1999; Neuenschwander et al., 1999; Stephens et al., 2006), precluding the use of standard visual stimuli that are matched to the preferred spatial frequency of the cells under study. In vitro preparations also lack long-range spatiotemporal correlations likely to result from high-frequency fixational eye movements (Martinez-Conde et al., 2004), which theoretical models indicate may tap into high-frequency resonant circuitry (Miller, Denning, George, Marshak, & Kenyon, 2006). Even in the absence of a contribution from resonance circuitry in the inner retina, the present findings suggest that the spatiotemporal correlations generated by fixational eye movements are by themselves sufficient to convey precise pixel-by-pixel intensity information in an extremely synergistic manner. A fundamental assumption of the extreme synergy hypothesis is that retinal mechanisms exist to produce useful, long-range correlations in response to natural stimuli. 

Separate from the issue of whether long-range spatiotemporal correlations can be evoked under physiological conditions and, if so, what additional information they might encode, there remains the question of how such information might be utilized by downstream elements. Both theoretical (Abeles, 1982; Kenyon, Fetz, & Puff, 1990; Kenyon, Puff, & Fetz, 1992) and experimental (Usrey, Alonso, & Reid, 2000) studies indicate that downstream neurons could operate as coincidence detectors, suggesting that spatiotemporal correlations of retinal origin directly influence the responses of targets elsewhere in the CNS (Castelo-Branco, Neuenschwander, & Singer, 1998; Doty & Kimura, 1963). Unsupervised, activity-dependent Hebbian plasticity rules, operating at individual synapses, have been used to train nodes in artificial neural networks to respond selectively to principal components (Hyvarinen, Karhunen, & Oja, 2001). Combined with nonlinear dendritic subunits (Gasparini & Magee, 2006; Oesch, Euler, & Taylor, 2005; Polsky, Mel, & Schiller, 2004), which may allow selective responses to pairs of synapses activated in close spatial and temporal proximity, there exists at least the possibility of extracting principal components defined by the spatiotemporal correlations present in short segments of spike train data. Sensitivity to oscillatory input, by exploiting resonances at the synaptic (Beierlein & Connors, 2002; Wespatat, Tennigkeit, & Singer, 2004), cellular (Vigh, Solessio, Morgans, & Lasater, 2003), and network levels (Miller et al., 2006; Rager & Singer, 1998), could extend the sensitivity to synchronous inputs so as to emphasize periodic temporal structure. 

The psychophysical literature does not appear to include any studies that directly address the plausibility of a retinal code based on extreme synergy, although studies have indicated a role for fixational eye movements (Rucci, Iovin, Poletti, & Santini, 2007a). Fast displays have been used to manipulate the spatiotemporal correlations between pixels, and by implication, the spatiotemporal correlations between associated visual neurons, in order to study the neural mechanisms underlying the perceptual grouping of distributed image features (Greene, 2007; Usher & Donnelly, 1998). Analogous procedures might be adapted to test the predictions of the extreme synergy hypothesis by measuring changes in psychophysical thresholds as a function of externally imposed spatiotemporal correlations, although fixational eye movements would likely have to be precisely accounted for to permit a robust interpretation (Rucci et al., 2007a). 

Finally, it is reasonable to speculate as to why the retina might have evolved a coding strategy based on extreme synergy in which local intensity information is distributed across many surrounding neurons. A possible clue is that to encode similar amounts of information using independent firing rates would require much higher levels of activity than are observed physiologically. This line of reasoning suggests that the primary advantage of extreme synergy is to convey large amounts of information using a minimum number of spikes. Metabolic constraints have been previously suggested to influence the retinal code (Laughlin, de Ruyter van Steveninck, & Anderson, 1998). Given that the 2-point correlation function of natural scenes extends out to several degrees of visual angle (Reinagel & Zador, 1999), well beyond the central excitatory portion of a typical ganglion cell's central receptive field, it follows that individual retinal neurons are usually stimulated as part of larger visual features. The statistics of natural scenes therefore may permit the retina to distribute information among multiple neurons representing instantaneously defined ensembles of functionally related cells. 

Conclusions

The present findings suggest that the retina employs a distributed spatiotemporal code that utilizes extreme synergy, in which the firing events from any one neuron can only be properly interpreted in the context of multiple surrounding events. Such a code may have evolved to conserve energy by conveying maximal information using a minimal number of spikes. 

Acknowledgments

This research was supported by LDRD program at the Los Alamos National Laboratory and by the National Science Foundation (#0749348). The author is deeply grateful to David Marshak, Ilya Nemenman, Alex Nugent, and John Mosher for their comments and suggestions and the editor and reviewers for their detailed contributions to revised versions of the manuscript. 

Commercial relationships: none. 

Corresponding author: Garrett T. Kenyon. 

Email: [e-mail]. 

Address: Physics Division, Los Alamos National Laboratory, P-21, P.O. Box 1663, MS D454, Los Alamos, NM 87545, USA.