Inter-mosaic coordination of retinal receptive fields.

The output of the retina is organized into many detector grids, called 'mosaics', that signal different features of visual scenes to the brain1-4. Each mosaic comprises a single type of retinal ganglion cell (RGC), whose receptive fields tile visual space. Many mosaics arise as pairs, signalling increments (ON) and decrements (OFF), respectively, of a particular visual feature5. Here we use a model of efficient coding6 to determine how such mosaic pairs should be arranged to optimize the encoding of natural scenes. We find that information is maximized when these mosaic pairs are anti-aligned, meaning that the distances between the receptive field centres across mosaics are greater than expected by chance. We tested this prediction across multiple receptive field mosaics acquired using large-scale measurements of the light responses of rat and primate RGCs. ON and OFF RGC pairs with similar feature selectivity had anti-aligned receptive field mosaics, consistent with this prediction. ON and OFF RGC types that encode distinct features have independent mosaics. These results extend efficient coding theory beyond individual cells to predict how populations of diverse types of RGC are spatially arranged.

Main

Mosaics are a fundamental organizing feature of retinal output[2,7]. This output is formed by RGCs, which transmit signals to the brain down their axons. Most, if not all, RGC types form mosaics, meaning the RFs of each type are regularly spaced and approximately tile the retina. This produces a complete, uniform and efficient encoding of visual scenes[1-3,6,8]. Furthermore, many RGC types come in functionally similar pairs, with one member signaling increments of a particular visual feature and the other signaling decrements. These functionally related pairs also serve to efficiently encode visual stimuli[9,10], and frequently form ‘pathways’ through the brain that converge downstream. For example, ON and OFF parasol RGCs in primates, which form the magnocellular pathway and ultimately converge in primary visual cortex[5]. The convergence of ON and OFF signals strongly shapes cortical processing, and possibly perception[11-14]. This organization of visual signaling raises the question, how should the mosaics of different RGC types be arranged to efficiently encode natural scenes? Answering this question is important because visual processing ultimately depends not on the activity of a single cell or cell type, but on the activity of large populations of diverse cell types.

There are many ways in which two mosaics can be spatially arranged, however we focus on three broad categories: independent, aligned and anti-aligned. For a pair of independent mosaics, the RF locations in one mosaic provide no information about RF locations in the other mosaic. This organization is commonly assumed based on anatomical measurements of retinal cell bodies[15-17], but cell body locations weakly constrain RF locations. For aligned mosaics, the RF locations across mosaics would tend to be closer than expected by chance. Finally, for anti-aligned mosaics, the RF centers across mosaics would tend to be farther apart than expected by chance. Each arrangement has potential advantages. Alignment puts two detector grids in spatial registration, potentially facilitating the encoding or decoding of stimuli. Anti-alignment may provide higher spatial resolution or make the system more robust to noise[18]. Independence may strike a balance between these two states, or simply be straightforward to achieve developmentally.

Efficient Coding Predicts Anti-Alignment

We begin by examining the predictions of a simple model of retinal processing that is optimized according to ECT, based on previous work[6]. The model consists of an ensemble of linear filters (representing RFs) and rectifying nonlinearities (representing the transformation from filtered visual input to spike rates). Input and output noise are included to model noise in phototransduction and spike generation, respectively. Finally, the model includes a cost function for generating spikes (Fig. 1a). We optimized the linear filters and nonlinearities to maximize the mutual information between a library of 10,000 natural images and simulated spike counts of the model retina. Consistent with previous results[6], the optimized model converged on center-surround RFs organized into two mosaics: a mosaic of ‘ON’ RFs and a mosaic of ‘OFF’ RFs (Fig. 1b, c, d). This division of ON and OFF RFs reduces the spiking output of the system[6,9]. Other than the polarity difference, the RFs and nonlinearities were similar across the two mosaics, indicating these mosaics mimic a pair of cell types that encode similar visual features, but with opposite polarities.

Figure 1:

Efficient coding predicts anti-alignment between ON and OFF RGCs with similar feature selectivity.

(a) The schematized model is trained on ~10,000 natural images to maximize the mutual information (MI) between input and collection of responses, (r1, r2,…, rj) under noise and metabolic constraints (see Methods). (b) Noisy spatial filters at the beginning of training (left) become center-surround filters by end of training (right). The two learned filters exemplify ON (positive contrast) and OFF (negative contrast) RFs. Change in MI with training (center). (c) All 100 filters after optimization. Each square patch is 18×18 pixels (orange box). (d) The center of mass (filled circles) and 1 s.d. contour around the center of mass (solid line) are shown for n=46 ON-center (green) and n=54 OFF-center (magenta) filters, on the 18 × 18 pixel grid. These define the ON and OFF mosaics. (e) Center points from all ON-center (green) and OFF-center (magenta) model neurons. (f, g) 2-D z-scored IMCE map (f) and z-scored radial average IMCE (g) for the ON and OFF mosaic pairs in (e) show that the mosaics generated from the efficient coding model are anti-aligned (see Fig. 2). Shaded area is s.d.

We next analyzed the relative spatial arrangements of RFs across the pair of optimized mosaics to determine if they were independent, aligned or anti-aligned. To distinguish between these possibilities, each RF was reduced to its center-of-mass (COM), thus reducing the mosaics to two grids of points (Fig. 1e). The relative spacing between points across the two mosaics was quantified by the sum over the squared inverse distances between all pairs of heterotypic points, which was restricted to a maximum value for points closer than a distance r (see Methods). This distilled the inter-mosaic relationship down to a single number, that we call the ‘inter-mosaic coordination energy’ (IMCE). The IMCE is motivated by the idea of a potential energy function that weights local interactions more strongly than distal interactions[11].

We then rigidly shifted one mosaic with respect to the other and recomputed the IMCE for each Δx and Δy shift (Fig. 2e), to compare the observed arrangement with alternative arrangements at different shifts. The IMCE across all shifts are then z-scored (normalized by the s.d. of observed IMCEs). The distribution of IMCE approximately followed a normal distribution, confirmed by kurtosis, skewness and percentile measurements (Extended Data Table 2).

Figure 2:

Analysis framework for measuring mosaic coordination.

(a-c) Pairs of ON (green) and OFF (magenta) mosaics generated from modified PIPP model with specific constraints on homotypic (ON-ON, OFF-OFF) and heterotypic (ON-OFF) interactions (see Methods) exhibit anti-alignment (a), alignment (b) and independence (c) respectively. Each point represents the center of mass of the spatial RF. X and Y axis units are arbitrary. Number of points for each of ON and OFF mosaic: n=50. (d) The repulsive energy between two points as a function of their separation ‘r’. The energy is set to a constant value for r≤rmin (see Eq. 3). Inset: Inverse cube repulsive force, F, and potential energy, E, between two heterotypic points (green and magenta). (e) ON mosaic (green) rigidly shifted with respect to OFF mosaic (magenta), by different amounts (length of arrow) and along different directions (direction of arrow), for estimating 2-D ‘inter-mosaic coordination energy’ (IMCE) map. (f, g, h) The 2-D z-scored IMCE maps for anti-aligned (f), aligned (g) and independent (h) mosaic pairs (n=1). X- and Y-axis have units of normalized shift. The z-scored IMCE is averaged over an annulus (shown in f) at different radial distances to obtain the radial z-scored IMCE curve. (i, j, k) An average over n=100 radial z-scored IMCE curves, each obtained by radially averaging the 2-D z-scored IMCE map for anti-aligned (i), aligned (j) and independent (k) mosaic pairs. The radial z-scored IMCE curve for the 2D z-scored IMCE map (h) is shown in blue (k). Shaded areas are s.e.m.

Extended Data Table 2 |

The distribution of inter-mosaic coordination energy (IMCE) values is approximately normal.

Mosaic pairs	16th percentile	(μ − σ) †	84th percentile	(μ + σ)
PIPP* simulateci mosaic pairs (Fig. 2a)	0.643	0.638	0.884	0.894
Mosaic pairs from efficient coding model (Fig. 1e)	0.283	0.241	0.685	0.686
Brisk transient mosaic pairs from rat retina (average± s.e.m. over 5 retinas: Fig. 3g)	0.025±0.001	0.025±0.001	0.033±0.002	0.034±0.002
Brisk sustained mosaic pairs from rat retina (average± s.e.m. over 3 retinas: Fig. 3i)	0.033±0.001	0.032±0.001	0.063±0.009	0.064±0.010
Parasol mosaic pairs from primate retina (average± s.e.m. over 3 retinas: Fig. 3k)	0.049±0.001	0.049±0.002	0.056±0.002	0.056±0.002

PIPP: Pairwise interaction point process

μ, σ: mean and s.d. of IMCE values

To gain an intuition about coordination from the IMCE, we used a modified PIPP model[19-21] to generate two sets of points (types 1 and 2) with specific homotypic and heterotypic interactions (see Methods). This allowed homotypic pairs of points to maintain a regular spacing observed in mosaics, while generating mosaic pairs that were anti-aligned (Fig. 2a), aligned (Fig. 2b), or spatially independent (Fig. 2c). If two mosaics are anti-aligned (Fig. 2a), the IMCE tends to increase as one mosaic is rigidly shifted with respect to the other, yielding a lower IMCE for arrangements (Fig. 2f). The radial average of the two-dimensional IMCE map yields a function that has a trough near zero (Fig. 2i). Conversely, if two mosaics are aligned, the IMCE tends to decrease as a function of shift distance (Fig. 2b, g, j). For independent mosaics (Fig. 2c), the IMCE map doesn’t exhibit a consistent topographical structure (Fig. 2h), and hence the radial average is relatively flat as a function of shift distance (Fig. 2k). Note, we only use this ‘energy’ analysis to examine coordination between mosaics; we do not postulate that an actual energy exists between mosaics or was involved in their development.

Applying this analysis to the ON and OFF mosaics optimized according to ECT indicated that they were anti-aligned (Fig. 1e, f, g). The robustness of this result was confirmed by (1) training the model on a set of independently collected natural images (see Methods), (2) using more or less severe cost functions on spiking, (3) changing the number of filters, and (4) forcing the number of filters with ON and OFF centers respectively to be unequal (Extended Data Fig. 1), including vastly different cell densities in the simulated mosaics (Extended Data Fig. 2). The optimized mosaics from the efficient coding model remained anti-aligned under these manipulations.

Extended Data Figure 1 |

Anti-alignment predicted by efficient coding theory is conserved in mosaics with different densities and boundary conditions.

(a-d) Optimal spatial filters of 48 ON units and 52 OFF units, each on a 18×18 pixel grid (orange box) (a). The COMs of optimal filters forming the ON (green) and OFF (magenta) mosaics. Training was performed using a circular mask over the images to reduce edge artifacts (b). The 2-D z-scored IMCE map (c) and radial average z-scored IMCE (d) for the mosaic pairs shown in (b). (e-h) Same as in (a-d) with equal cell density (n=50) for ON and OFF mosaics. (i-l) Same as in (a-d) with number of ON and OFF units fixed at n=45 and n=55 respectively. (m-p) Same as in (a-d) with n=49 ON units and n=51 OFF units, however training was performed without a circular mask. Shaded areas are s.d.

Extended Data Figure 2 |

Mosaic coordination can persist under widely diverging RF densities.

(a, b, c) Bivariate point pattern (type 1: green, type 2: magenta) generated by modified pairwise interaction point process model (PIPP; see Methods) with interaction terms for anti-alignment (a), alignment (b) and independence (c). The density of type 2 points is four times higher than the density of type 1 points. (d, e, f) The 2-D z-scored IMCE maps corresponding to (a), (b) and (c) respectively. (g, h, i) The radial average z-scored IMCE averaged over n=100 mosaic pairs that are aligned (g), anti-aligned (h) and independent (i). Shaded areas are s.e.m.

ON and OFF RGC Mosaics are Anti-Aligned

Following these predictions from ECT, we next tested if measured RGC mosaics exhibited this predicted spatial organization. To measure the spatial relationships across RGC mosaics, we used a large-scale multi-electrode array (MEA) to map the RF locations of diverse RGC types[22,23]. A spatiotemporal white noise stimulus was presented to ex vivo rat retina and the spike-triggered average stimulus was calculated for each RGC to estimate its spatial and temporal RF[24]. A typical experiment yielded spikes from ~300–400 RGCs that were functionally classified according to their spatiotemporal RFs and other response properties (see Methods)[23,25]. RF locations were not utilized in the classification (see Methods), yet the RFs of each clustered RGC type exhibited a mosaic-like arrangement (Fig. 3a, d) - confirming functionally distinct and irreducible RGC types[1,23,26,27].

Figure 3:

Receptive field mosaics of functionally paired ON and OFF RGC types are anti-aligned.

(a) Representative ON (left: light blue shaded) and OFF (right: gray shaded) RF mosaics of brisk transient RGCs from rat (blue). RF contours are at 61% of peak (~1-s.d. for a Gaussian RF). Solid black circles show center-of-mass (COM). Inset: spatial RF of an individual RGC with the contour and COM. (b) The COMs of RFs of all ON (blue) and OFF (gray) brisk transient RGCs. Black contour shows region of interest (ROI) used in analysis (see Methods). (c) Left: 2-D z-scored IMCE map as a function of shift between ON and OFF brisk transient mosaics (a). Mosaic shift is normalized to the mean homotypic nearest-neighbor distance. Right: radial average z-scored IMCE. (d-f) Same as (a-c) but for ON and OFF brisk sustained RGCs (magenta). (g) RF mosaics of ON and OFF brisk transient RGCs from n=5 retinas; each column from the same retina. (h) Radial average z-scored IMCE as a function of mosaic shift for 5 mosaic pairs shown in (g). (i, j) Same as (g, h) for n=3 pairs of ON and OFF brisk sustained RGCs, each from a different retina. (k) RF mosaics of ON (top) and OFF (bottom) parasol RGCs from n=3 primate retinas (green). (l) Left: 2-D z-scored IMCE map of a representative primate retina (indicated by box in k). Right: radial average z-scored IMCE for the 3 primate retinas as a function of mosaic shift. Dashed curve corresponds to the 2D z-scored IMCE map on the left. Shaded areas are s.d.

We focused on four RGC types for which we frequently measured nearly complete mosaics (Fig. 3a, d, g, i). These four types comprised of two pairs of ON and OFF cells that were otherwise functionally similar: brisk transient and brisk sustained RGCs[23]. Thus, they were analogous to the pair of types generated by the efficient coding model (Fig. 1) and allowed us to examine the relationships between multiple pairs of real RGC mosaics that either encode similar (up to a polarity flip) or distinct visual features.

First, we examined ON and OFF brisk transient RGCs: these types both encode temporally bandpass stimuli and exhibit strongly rectified contrast responses[23]. The IMCE maps were less periodic than those generated from ECT (Fig. 1) because the biological mosaics were more disordered. Nevertheless, the IMCE between these RGC mosaics rose with shift distance (Fig. 3b, c), indicating that they were anti-aligned. Next, we examined ON and OFF brisk sustained RGCs, which exhibit relatively low-pass temporal filtering and more linear contrast responses. These mosaics also exhibited an IMCE that rose with shift distance (Fig. 3e, f), indicating anti-alignment. These results were consistent across multiple mosaics measured from different rat retinas (Fig. 3g, h, i, j). Thus, ON and OFF RGC types that form functional pairs exhibit anti-aligned mosaics, consistent with predictions from ECT (Fig. 1f, g).

If anti-alignment between ON and OFF RGC mosaics is driven by the efficient encoding of natural scenes, it should persist in other terrestrial species. Thus, we examined the spatial relationships between mosaics of ON and OFF parasol cells in the primate retina[28]. An analysis of IMCE also revealed anti-alignment between these functionally similar (but opposite polarity) RGC types (Fig. 3k, l). Thus, anti-alignment of ON-OFF mosaics generalizes across RGC types in at least two species.

To test the robustness of these results to mosaic imperfections such as gaps caused by missing RFs, we examined the extent to which artificially adding or removing RFs from measured mosaics impacted these results (Extended Data Fig. 3a). We found that the radial dependence of the IMCE was robust up to a ~20% reduction in the number of RFs in the mosaic (Extended Data Fig. 3b, c). Artificially adding RFs to ‘gaps’ in the measured mosaics reduced variability in the IMCE estimates but did not alter the qualitative results for mosaic anti-alignment observed in rat or primate retinas (Extended Data Fig. 3b, c).

Extended Data Figure 3 |

Mosaic coordination estimates are robust to RF subsampling.

(a) RF mosaics illustrating three different cases. ‘Measured’: no RFs are removed or added; ‘depleted’: a fraction of randomly selected RFs removed (dashed ellipses); and ‘filled’: RFs artificially added to fill mosaic gaps (thick solid ellipses). The gradient (bottom) illustrates the percentage of RFs remaining after removing or adding RFs. (b-c) The radial average z-scored IMCE is shown for different amounts of subsampling and filling of ON and OFF bt (rat) mosaics (b; blue), ON and OFF bs (rat) mosaics (b; purple), and ON and OFF parasol (primate) mosaics (c; green). bt: brisk transient, bs: brisk sustained. Each curve corresponds to an individual mosaic pair. Results are representative of n=5 retinas for ON-OFF bt, n=3 retinas for ON-OFF bs, and n=3 retinas for ON-OFF parasol RGCs. The percentage of RFs relative to measured (100%), is indicated by ‘f’. Shaded areas are s.d.

We next assessed the likelihood that the observed anti-alignment between mosaics occurred by chance. We constructed a null distribution for mosaic coordination by combining ON and OFF RGC mosaics measured from different retinas of rats and primates. These ‘pseudo-pairs’ of mosaics are statistically independent, and therefore should yield radial IMCE curves distinct from those of aligned or anti-aligned mosaics. We spatially registered each pseudo-pair of mosaics based on their center of mass (COM, see Methods). The distribution of IMCE across pseudo-pairs remained largely unaltered for small offsets between the COM. To quantify anti-alignment, we calculated the weighted area under the radial IMCE profile curve (see Methods) for both real and pseudo-pairs (Fig. 4a). The real pairs exhibited values that were outside the confidence bounds of the null distribution generated from the pseudo-pairs (Fig. 4b, c; Extended Data Table 1), indicating a low probability that the observed anti-alignment in the measured mosaics arose by chance.

Figure 4:

Mosaic anti-alignment is unlikely to have arisen by chance.

(a) Example mosaics of ON (left) and OFF (right) parasol RGCs, from two retinas. Mosaics from the same retina (real pairs) are indicated by red arrow; mosaics from different retinas (pseudo-pairs) are indicated by black arrow. Inset: radial average z-scored IMCE (open circles) for real (red) and pseudo (gray) mosaic pairs. Parametric fit (solid line) and area under curve (shaded region) are used to estimate coordination index (CI; see Methods). X-axis unit is normalized radial shift. (b) Sampling distribution from bootstrap estimates of mean CI for pseudo pairs (gray) and mean CI for real pairs (filled colored circles) for ON-OFF brisk transient: bt (left), and ON-OFF brisk sustained: bs (right). Number of real, pseudo pairs = 5, 40 (left); 3, 12 (right). CI mean±s.d. = 0.39±0.13 (left), 0.33±0.19 (right). The gray shaded region to the right of the vertical dashed line indicates value exceeding 95% confidence interval based on one-sample two-sided t-test statistic. (c) Same as (b) for ON-OFF parasol mosaics. Number of real, pseudo pairs = 3, 12. CI mean±s.d.=0.68±0.13. P-values and Cohen’s ds are: 3.03×10−8 and 1.61, 3.72×10−3 and 1.23, and 8.78×10−48 and 1.19, for ON-OFF bt, ON-OFF bs and ON-OFF parasol respectively. Results are representative of n=5 retinas for ON-OFF bt (rat), n=3 retinas for ON-OFF bs (rat), and n=3 retinas for ON-OFF parasol RGC types (primate).

Extended Data Table 1 |

Coordination Index values for real and pseudo mosaic pairs.

Species	RGC types *	Coordination	Real pairs (μ ± σ) †	Pseudo pairs (μ ± σ) †
	ON bt - OFF bt	Anti-aligned	0.393 ± 0.134	0.002 ± 0.385
	ON bs - OFF bs		0.326 ± 0.186	− 0.054 ± 0.360
Rat	ON bt - ON bs		0.070 ± 0.437	− 0.051 ± 0.411
	OFF bt - OFF bsON bt - OFF bs	Independent	− 0.076 ± 0.715 − 0.124 ± 0.330	− 0.078 ± 0.305 − 0.211 ± 0.324
	ON bs - OFF bt		0.099 ± 0.360	0.179 ± 0.366
Primate	ON parasol - OFF parasol	Anti-aligned	0.683 ± 0.131	0.049 ± 0.365

bt: brisk transient, bs: brisk sustained

μ, σ: mean and s.d. of Coordination Index

Mosaics are Independent Across Pathways

Finally, we examined whether inter-mosaic coordination was present between RGC types encoding distinct visual features (e.g. ON brisk transient and OFF brisk sustained or ON brisk transient and ON brisk sustained). Mosaic pairs that differed in their feature selectivity appeared statistically independent (Extended Data Fig. 4a, b, f, d, h; Extended Data Table 1); this was confirmed by the failure of the coordination index of these pairs to fall outside the null distribution established from pseudo-pairs (Extended Data Fig. 4c, g, e, i). Thus, mosaic anti-alignment appears restricted to pairs of RGC types that form pathways encoding similar visual features with ON and OFF responses.

Extended Data Figure 4 |

Mosaics encoding distinct visual features appear to be independent.

(a) Example mosaics of ON and OFF brisk transient and brisk sustained RGC types. Coordination was tested across cell type (orange), and across cell type and polarity (green). (b, d, f, h) 2D z-scored IMCE map of a representative pair (left) and radial average z-scored IMCE of all pairs (right), for ONbt-ONbs (b, n=3), OFFbt-OFFbs (d, n=3), ONbt-OFFbs (f, n=3), ONbs-OFFbt (h, n=3) mosaic combinations; bt: brisk transient, bs: brisk sustained. The radial average z-scored IMCE corresponding to the 2D z-scored IMCE map (left) is shown by dashed curve. Shaded areas are s.d. (c, e, g, i) Sampling distribution from bootstrap estimates of mean coordination index for pseudo pairs (gray) and real pairs (orange/green filled circle, arrow). Number of pseudo pairs: n=12 (c), n=12 (e), n=12 (g), n=12 (i). The gray shaded region to the right of the vertical dashed line indicates value exceeding 95% confidence interval based on one-sample two-sided t-test statistic: P value = 0.33, 0.98, 0.37 and 0.46 respectively for (c), (e), (g) and (i); (n.s., not significant). Cohen’s d = 0.36, −0.006, 0.28 and −0.25 respectively for (c), (e), (g) and (i).

Discussion

We have demonstrated a novel prediction from ECT: ON and OFF mosaics with similar feature selectivity should be anti-aligned to optimize the encoding of natural scenes (Fig. 1). We showed that measured RGC mosaics conform with this prediction across cell types and species (Figs. 3 and 4). This indicates that instead of each RGC type representing an independent and parallel pathway of visual signaling, at least some pairs of RGC types have spatially coordinated RF mosaics, and this coordination likely supports the efficient encoding of natural scenes.

We speculate that anti-alignment is preferred because it mitigates the effects of noise. When input and output noise in the efficient coding model were small (e.g. the variability was less than the mean) the optimized mosaics were aligned. We think the higher noise state is more likely to be representative of retinal processing because of many Poisson-like noise sources in the retina. For example, at the input stage of retinal processing, the absorption of photons is Poisson. Within photoreceptors, the spontaneous activation of phototransduction molecules along with Poisson-like variability in vesicular glutamate release, all add significant amounts of noise[29-32]. After the photoreceptors, noise is further added by bipolar cells and in spike generation by RGCs[33,34]. The input and output noise in the efficient coding model are not easily mapped onto these biological noise sources. However, retinal noise is substantial and is thought to strongly influence mammalian vision[6,35-37]. Furthermore, noise is likely to dominate signaling under conditions that are most survival relevant, such as detecting camouflaged predators or prey. Thus, anti-alignment may be a strategy to increase the reliability of coding when noise dominates signal.

How much do anti-aligned mosaics improve the encoding of natural scenes? To gain some insight to the question, we compared the mutual information between a test set of natural scenes and two mosaic pairs: one ON-OFF pair that was anti-aligned, and a second pair that was aligned. The aligned pair encoded ~4% less information than the anti-aligned pair (see Methods). For context, eliminating the RF surrounds in two mosaics of anti-aligned ON and OFF RGCs resulted in ~19% decrease in mutual information. Thus, inter-mosaic coordination represents ~20% of the encoding benefit provided by RF surrounds, a feature of retinal encoding that is broadly considered fundamental to early visual processing.

One important caveat to these results is that RF structure can depend on light level and other stimulus variables[38,39]. Thus, anti-alignment could also depend on stimulus conditions. In one experiment, we tested the persistence of anti-alignment from the cone-mediated light levels examined here (Fig 4) to rod-mediated light levels (Extended Data Fig. 5). The anti-alignment among ON and OFF BT RGCs persisted: ON and OFF BS mosaics were insufficiently complete in this experiment for a conclusive test. More experiments are needed that examine a wider range of stimulus manipulations, but these data suggest mosaic anti-alignment persists across light levels.

Extended Data Figure 5 |

Anti-alignment between ON brisk transient RF mosaics persists across light levels.

(a) RF mosaics of ON (left) and OFF (right) brisk transient RGCs, measured at photopic light level (10,000 Rh*/rod/s). The COM of RFs are indicated by black filled circles. (b-c) The 2D z-scored IMCE map (b) and radial average z-scored IMCE (right) for the mosaic pair shown in (a). (d) RF mosaics of ON (left) and OFF (right) brisk transient RGC types, at scotopic light level (1.0 R*/rod/s), from the same retina as in (a). (e-f) Same as in (b-c) for the mosaic pairs shown in (d). (g) Change in RF COMs of ON brisk transient RGCs from photopic to scotopic light level (black filled circles). Homotypic nearest-neighbor (NN) distance between RFs estimated at photopic and scotopic light levels is shown by solid and dashed red lines, respectively. (h) Distribution of fractional change in RF position of ON brisk transient RGCs across light levels expressed as a fraction of the mean NN homotypic distance at photopic light level. The smooth curve is kernel density estimate. (i-j) Same as in (g-h) for OFF brisk transient RGCs. Results are representative of n=1 retina. Shaded areas are s.d.

The second caveat is that the ECT model used here is optimized on natural scenes, not natural movies. As such, it is tractable, but only generates two ‘cell types’: an ON and an OFF type. However, the mammalian retina contains many RGC types that can encode distinct features (e.g. direction selectivity versus chromatic selectivity). Whether there exists any coordination between these RGC mosaics is yet to be determined, but it motivates further exploration of optimizing a model retina on natural movies with chromatic content. This optimization may need to include head and eye movements with global and local motion to mimic animals moving through natural environments.

The anti-alignment of ON and OFF mosaics may support aspects of visual processing beyond efficient coding, such as the formation of orientation tuning in V1[11,40]. A V1 neuron can receive input (indirectly) from two randomly selected nearby ON and OFF RGCs and form an orientation tuned RF if the mosaics are anti-aligned.

Finally, these results raise new questions about the developmental and synaptic mechanisms that produce inter-mosaic coordination. For example, does this coordination require visual experience, or is it genetically programmed? While these questions await further investigation, our results show that distinct cell types can exhibit exquisite and highly specific spatial coordination within a neural circuit, conferring advantages for the efficient encoding of natural stimuli.

Methods

Animals and recording

All procedures involving the use of animals were approved by the Institutional Animal Care and Use Committee. Macaque monkeys (Macaca mulatta) and Long-Evans (LE) rats were used in these experiments. Preparation of primate retinas follows previously described protocols[41,42], and primate data was generously provided by E.J. Chichilnisky. For experiments with rat retinas, animals between the ages of 3 months and 1 year were used. Both main and female primates and rats were included in the study. Animals were housed in 12-hour light-dark cycle with ad lib access to food and water. Rats were dark adapted for 12–16 hours prior to the experiment. Rats were deeply anesthetized by intraperitoneal injection of Ketamine (0.1mg/10mg of body weight) and Xylazine (0.1mg/20mg of body weight), confirmed by extinguished eyeblink and toe-pinch reflex and then euthanized by decapitation. The eyes were immediately enucleated and placed in a petri dish filled with bicarbonate buffered Ames media (Sigma, St. Louis, MO) bubbled with 5% oxygen and 95% carbon-dioxide maintained at 35–36°C and pH 7.4. Following hemi-section and vitrectomy, the retina was detached from the sclera and a small piece ~1.5mm x 3mm was cut from the medial-dorsal region and placed on the multielectrode array (512 electrodes, 60 μm pitch), RGC side down. Euthanasia and retinal dissections were performed in a completely dark room with the assistance of infrared lamps and infrared converting goggles.

Spikes on the MEA were identified by previously described methods[23,43]. Briefly, extracellular voltages were sampled at 20kHz and stored for offline analysis. Spikes were dually sorted using a JAVA based custom automated algorithm and YASS[44,45], and principal component analysis (PCA) was used to reduce the dimensionality of the data and a mixture of Gaussians model was fit to extract clusters of spikes. Errors in spike clustering were corrected by manually adjusting the location and number of Gaussians fit to the data[46]. Each spike cluster was assigned to a putative RGC if the ensemble of spikes exhibited a refractory period of at least 1.5ms with < 10% estimated contamination and at least 100 spikes. Individual RGCs were tracked across multiple stimulus conditions by checking for correspondence in clustered spike shapes and the electrophysiological image[47,48].

Visual stimuli and classification of RGCs into distinct functional types

A gamma corrected OLED display (SVGA+3XL, Emagin Corporation, Santa Clara, CA) refreshing at 60.35 Hz was used to display visual stimuli to the retina. The stimuli were focused on the photoreceptor outer segments using an inverted microscope (Nikon Eclipse TE-2000U). The mean intensity of the visual stimulus at the photoreceptors was ~7000 photoisomerizations/M-cone/sec for rat retinas, and ~1700 photoisomerizations/M-cone/sec for monkey retinas, at photopic light level. At scotopic light level, the intensity values were ~1.7 photoisomerizations/rod/sec respectively for these retinas.

RGCs were classified into distinct types using previously described methods[23,24]. We used binary checkerboard noise, with ~40μm x 40μm checkers for rat retinas and ~55μm x 55μm checkers for monkey retinas. The images refreshed either every one or two frames. We computed the spike-triggered average (STA) stimulus as an estimate of the spatiotemporal RFs of individual RGCs[24]. In conjunction, response properties such as spike rate, autocorrelation and contrast selectivity were used to serially classify RGCs into different groups. The union of RFs obtained from the automated spike sorting algorithm and YASS[44,45] was taken to generate the complete RF mosaic for each group of classified RGCs. Identification of irreducible RGC types was confirmed by observing a mosaic-like arrangement of spatial RFs.

Estimating RF locations and plotting mosaics

The spatial RF of individual RGCs was extracted as follows: Stimulus pixels (checkers) across frames of the STA that exceed four standard deviations above or below the mean intensity across all stimulus pixels were selected. The average time course of pixel intensity values of the selected pixels was computed as an estimate of the temporal RF. The inner product of the temporal RF and each pixel of the STA across time was computed, which collapsed the spatiotemporal STA to a single frame; pixels that evolved in time in a manner similar to the temporal RF had high values and those that did not had small values. This single frame was taken as an estimate of the spatial RF[41]. This procedure assumes the spatiotemporal RF is the outer product of two filters: a spatial and a temporal filter. Singular value decomposition performed on the full spatiotemporal RF yielded similar results. Specifically, the first vector pair (those with the highest associated singular value) yielded an estimate of the temporal and spatial RFs that closely matched the procedure described above, but in general exhibited greater noise[23].

After obtaining the spatial RF, stimulus pixels exceeding 5% of the peak pixel were identified and used to calculate the center of mass of the spatial RF. To plot RF mosaics, the spatial RF was filtered by a Gaussian function with a standard deviation of 0.75 pixels (1.5 pixels for ON brisk sustained cells) and a contour at 61% (~1-s.d. for a Gaussian RF) of the peak pixel was extracted using the contour function in MATLAB (Mathworks, Natick, MA).

Analysis of mosaic coordination

To analyze RF coordination between mosaics of distinct RGC types we computed an inter-mosaic coordination energy (IMCE). We call it an ‘energy’ because it is analogous to postulating a repulsive ‘force’ between heterotypic cell types and integrating these forces to calculate an energy of the system. To compute IMCE, each spatial RF was summarized by a single point at the center of mass (COM) of the RF. Then, the relative spatial distances between the RF COMs across the two mosaics were assessed by computing the inverse of the squared distances between all heterotypic points (Fig. 2d). Specifically, for a pair of points i and j separated by a distance r, their contribution is: To prevent E from diverging as two points come very close to one another, E is set to a constant value if the distance between points is smaller than a threshold value set by r: The value for r was set to 0.20 times the median nearest-neighbor distance between heterotypic points[4] (see also Fig. 2d). Qualitative results did not depend strongly on the choice of this minimum distance. The values of E were then summed over all heterotypic pairs and divided by the total number of pairs to obtain: which provides a metric summarizing pairwise spacing between RF COMs, emphasizing nearby pairs and deemphasizing distant pairs.

The next step of the analysis was to compare the observed E to that for other possible nearby arrangements of the two mosaics, while keeping the intra-mosaic arrangement fixed. To achieve this, all points from one mosaic were rigidly translated while keeping the other mosaic fixed and recomputing the E. This process was repeated for different amounts of shift in different directions (Fig. 2e), yielding a 2D matrix given by, This 2D matrix was z-scored to normalize the values to their standard deviation and shift was normalized to the median homotypic nearest neighbor distance of points from the mosaic kept fixed, to generate a 2D map of IMCE values (Fig. 2 f, g, h): Radially averaging the IMCE map yields a curve that summarizes the spatial relationship between mosaics. If mosaics are aligned, then RFs are close to one another and the IMCE is high at zero shift and tends to fall when one mosaic is shifted with respect to the other mosaic (Fig. 2b, g, j). However, if mosaics are anti-aligned, then RFs are far from one another and the IMCE is relatively low at zero shift and tends to rise as one mosaic is shifted with respect to the other mosaic (Fig. 2a, f, i). To investigate the presence of higher order structures in the 2D IMCE maps, we analyzed the power spectrum and singular value decomposition of 2D IMCE maps from multiple retinas. A concentration of power mainly at low frequencies and a low rank decomposition of the 2D IMCE maps indicated that radial average captures the most dominant feature associated with mosaic coordination. However, absence of higher order structure may reflect limitations in our data (e.g. the number of RGCs sampled over the MEA).

Definition and controls for region of interest (ROI) choice

The lattices formed by the RF COMs are finite and relatively small – limited by the recording area of the MEA and the size of RFs. As a result, the continuity of RF arrangement breaks down near the edge of the measured mosaic. To mitigate potential bias in the estimated IMCE introduced by this discontinuity, we defined a region of interest (ROI) – a geometric area within which the density of RFs is approximately homogeneous for both mosaics and whose edges are at least 1.2 median homotypic nearest-neighbor distance from the boundary (Fig. 3 c, f). For a relatively periodic arrangement of RFs, a shift of 1 median inter-cell distance would be approximately equal to one spatial period, which when applied to one set of points, would move them such that they would occupy positions of their nearest neighbors in the direction of the shift. Shifts limited to <1.2 median homotypic nearest-neighbor distance within an ROI with spatially homogenous density of points, would (on average) preserve the spatial density of points inside the ROI for each shift, thus mitigating bias from discontinuity and variations in spatial density of points at the boundary. For a continuous lattice, exemplified by hexagonal lattices with spatial noise (Fig. 2a, b), the ROI can be a circle or any polygon whose boundary lies at least 1 median inter-cell distance inside the lattice boundary. We tested rectangular, convex hull and non-convex hull ROIs with boundaries that fell within 1.2 median inter-cell distances of the mosaic borders, and shift distances restricted to 1 median inter-cell distance. Results did not qualitatively depend on the choice of ROI shape.

Simulating ON and OFF mosaic coordination using a modified point-process (PIPP) model

We used a modified pairwise interaction point-process (PIPP) model to generate pairs of mosaics with specific types of coordination[21,49]. We first generated a mosaic consisting of N points, by adding points sequentially to a 2-D plane such that each new point lies outside a jittered radius centered around all the other existing points. Thus, the probability of adding a new point “1i” of type 1 at location is: where the homotypic exclusion distance d = d + ζ. The minimum distance d, is set by the areal density of points, such that no two points can lie closer than 2d. The value of the jitter parameter ζ is obtained from an independent draw from normal distribution with mean zero and standard deviation 0.12 ∗ d. This process is repeated until all N points are added to the plane.

To generate a second mosaic that is anti-aligned with the first (Fig. 2a), new points of type 2 are added to the same plane such that each new point “2i” satisfies the homotypic exclusion rule defined by Eq. 7 and the heterotypic exclusion rule: where the heterotypic exclusion distance d' = d + ζ'. The values of d and ζ' are set by the mean nearest-neighbor distance of homotypic points of first mosaic and a draw from normal distribution with mean zero and standard deviation 0.60 ∗ d. To generate a second mosaic that is aligned with the first, (Fig. 2b) the heterotypic exclusion rule is modified to: where, the values of d and ζ' are set by the mean nearest-neighbor distance of homotypic points of first mosaic and a draw from normal distribution with mean zero and standard deviation 0.35 ∗ d. To generate a second mosaic that is independent of the first, new points are added obeying only the homotypic exclusion rule (Eq. 7).

Significance testing for inter-mosaic coordination

To determine if the radial IMCE for the measured ON and OFF mosaic pairs are statistically distinct from those of independent mosaics, we utilized a permutation test by combining mosaics measured from different retinas. Two mosaics from different retinas (and hence different experiments), which we call ‘pseudo-pairs’, ought to exhibit no coordination, i.e., they are statistically independent on average. To estimate the IMCE for these pseudo-pairs of mosaics, we first spatially registered each mosaic pair. This was achieved by computing the center of mass of all the individual RF centers for each mosaic and then aligning the two COMs such that the two mosaics occupy the same region of space. Next, we generated a null distribution of radial IMCE for these pseudo mosaic pairs (Fig. 4a). The radial IMCE curves were further distilled to a single number we term the coordination index (CI) by first estimating the absolute value of the area under the curves. These areas were computed by first fitting the radial IMCE with a damped cosine function: where, a ≥ 0, ϕ, b, c and d are fitting parameters. The damped cosine function fit well the radial IMCE curves of both real and pseudo mosaic pairs, thereby providing an accurate estimate of the area under these curves (Fig. 4a). These areas were then converted into CI by assigning the areas a positive value if the value of E (r = 0) ≤ 0, and a negative value if E (r = 0) > 0. This caused radial IMCEs such as those in Figures 2i and 2j to have relatively large positive and negative CIs, respectively. The CIs for radial IMCE like that in Figure 2k were near zero.

For a given set of real mosaic pairs (e.g. ON and OFF brisk transient RGCs), the average of their CIs was computed. We then compared this average CI value to the distribution of bootstrapped averages over an equivalent number of samples from the distribution of CIs generated from pseudo mosaic pairs (Fig. 4b). Note, the distributions in Figure 4 are mean CI values over these equivalent samples from the ensemble of CI values derived from pseudo pairs. For example, if there were three real mosaic pairs, three samples were drawn from the CIs computed from pseudo mosaic pairs and these were averaged to generate an expected value from sampling three independent (pseudo) mosaic pairs. This process was repeated 1000 times (n=1000 bootstrap samples) to generate a null distribution of average CIs to which the observed CI was compared (e.g. Fig. 4b, c). If the observed CI lay outside the 95% confidence interval of the null distribution, the observed coordination was considered statistically significant.

Efficient coding model

The efficient coding model is an extension of previous work on the optimal coding of natural scenes[6]. Patches of natural scenes[50], each 18×18=324 pixels in size, are provided as inputs to the model. To reduce edge effects, we multiply each input image by a circular mask of a 9-pixel radius, which effectively leaves about 254 pixels in each input image patch. The patches are filtered by an array of linear filters, each corresponding to a model neuron. The total number of model neurons is kept fixed during training. For each model neuron, the output of the linear filter y is passed through a softplus nonlinearity f: to obtain the response r (Fig. 1a) with independently learnable parameters for threshold and slope in the softplus function: The response r represents the firing rate of the modeled RGC. For speed in training, the value of β is kept constant at 2.5. Two i.i.d. Gaussian noises with σ = 0.4 and σ = 3.0, are added to the stimulus and the response respectively. Additionally, we impose metabolic costs for generating spikes by constraining the mean firing rate to be 1 for each neuron. The qualitative results of the optimization did not depend on these parameter choices.

The model is trained using a stochastic gradient method (learning rate of 0.001) to maximize mutual information (MI) between input images and output firing rates. To constrain the mean firing rate, we use the augmented Lagrangian method[51] and reverse the sign of the partial derivatives with respect to the Lagrange multipliers during the gradient descent procedure. Filter weights were initialized as 18×18 white noise images (Fig. 1b) and the initial (log) slope and bias parameters of the nonlinearities were sampled from a uniform distribution [0, 1]. A batch of 100 model neurons was trained for 1,000,000 iterations. The filters were constrained to have unit L2 norm, which is achieved by normalizing the filters after each gradient step. To speed up convergence, we apply two perturbations between 200,000 to 500,000 iterations during training, which we call ‘jittering’ and ‘centering’. The jittering operation is applied every 5,000 iterations, where the filter weights are raised to the power of 1.25; this makes high-amplitude weights of the filters more pronounced while attenuating the low-amplitude weights. The centering operation penalizes the spatial spread of the filters, by adding the mean of the spatial variance of each filter as an additional loss term; this encourages the filter weights to be concentrated around the RF center. Thus, these perturbations can be conceived as a form of simulated annealing applied periodically to allow RFs to orient relative to one another. The eventual shape of the filters, their organization into mosaics, and the anti-alignment between ON and OFF mosaics did not depend on these jittering and centering procedures[6], but they sped convergence of the model.

To quantify the impact of mosaic coordination on the information encoded about natural scenes (see Discussion), we compared the MI that a pair of anti-aligned mosaics conveyed about natural scenes with that of a pair of aligned mosaics. The anti-aligned ON-OFF mosaics were obtained from the efficient coding model. The aligned mosaics were generated by first taking the RF centers of the ON mosaic from the anti-aligned mosaic pair and then duplicating these with a polarity flip to generate the OFF mosaic, such that the ON and OFF units were aligned. The center-surround structure of the RFs was described as a difference of two Gaussians with different standard deviations but with the same mean defined by the RF center. To examine the effect of eliminating the RF surrounds, the surround weights were set in the trained model to zero and the RF centers were unaltered. The mutual information between input images and firing rates of model neurons with these center only RFs was then compared to that from RFs with the original center-surround structure to assess the degree to which the RF surround improves natural scene encoding[6].