High-Fidelity Reproduction of Visual Signals by Electrical Stimulation in the Central Primate Retina

Experimental procedures.

Eyes were obtained from terminally anesthetized macaque monkeys (Macaca mulatta, Macaca fascicularis) of either sex (32 males, 6 females) used by other researchers, in accordance with Institutional Animal Care and Use Committee guidelines. Immediately after enucleation, the eye was hemisected in room light and the anterior portion of the eye and the vitreous humor was removed. The posterior portion of the eye was stored in oxygenated, bicarbonate-buffered Ames' solution (Sigma-Aldrich) at 31–33°C.

In infrared or dim lighting, small (∼2 × 2 mm) segments of retina, along with the attached retinal pigment epithelium (RPE) and choroid, were separated from the sclera. In most preparations, the retina was then separated from the RPE; in one preparation (Fig. 1B), the RPE was left attached. The retina was placed RGC side down on a custom high-density multielectrode array, and a transparent dialysis membrane was used to press the retina against the array from the photoreceptor side. A total of 52 retinal preparations was used: 15 from the raphe region and 37 from the periphery. For raphe preparations, regions 1–4.5 mm (4.5–20°) eccentricity along the temporal horizontal meridian were obtained. For peripheral preparations, segments in both temporal and nasal regions of the retina were obtained, with eccentricities ranging from 5 to 12 mm in eccentricity (22–56° temporal equivalent). Once the retina was mounted on the array, it was superfused with oxygenated, bicarbonate-buffered Ames' solution at 32–35°C.

Visual stimulation and recording.

The retina was typically stimulated with a white noise visual stimulus from a gamma-corrected cathode ray tube monitor refreshing at 120 Hz. The three monitor primaries were modulated independently for spectral variation, or in a coordinated manner for a black-and-white stimulus. The stimulus refresh interval was either 8.37, 16.74, or 33.47 ms and the size of each stimulus pixel at the photoreceptor layer ranged from ∼22–176 μm. For a single preparation (Fig. 1B), the retina was stimulated with a gamma-corrected organic light-emitting diode monitor refreshing at 60 Hz. In this case, the stimulus refresh interval was 16.57 ms and the size of each stimulus pixel at the photoreceptor layer was ∼36 μm. The recordings were obtained using a custom multielectrode stimulation and recording system with 512 electrodes (diameter, 5–15 μm) with either 30 μm pitch, covering an area of 0.43 mm², or 60 μm pitch, covering an area of 1.7 mm² (Litke et al., 2004; Hottowy et al., 2008, 2012). Raw signals from the channels were amplified, filtered (43–5000 Hz), multiplexed, digitized at 20,000 Hz, and stored for off-line analysis, as described previously (Litke et al., 2004). KiloSort (Pachitariu et al., 2016) was then applied to the raw recordings to separate and identify unique cells and their spike times.

Cell type classification.

To classify distinct cell types based on their light response properties, the spike-triggered average (STA) was computed for each cell (Chichilnisky, 2001). The receptive field (RF) of each cell was estimated by determining the pixels within the STA with a signal significantly above the noise level (Gauthier et al., 2009a). The time course of the STA was computed by summing the primaries of pixels within the RF. The spatial component of the STA was computed for each cell by taking the inner product between the time course and the full STA, and the output was fitted with a 2D Gaussian. Clusters of light response properties (RF diameter, first principal component of time course; Fig. 1) and spiking properties (first principal component of the spike train autocorrelation function; see below) were used to identify distinct cell types (Field et al., 2007). Recordings in most peripheral and some central retina preparations typically revealed the five numerically dominant RGC types: ON parasol, OFF parasol, ON midget, OFF midget, and small bistratified. In most preparations, recordings from ON and OFF parasol cells were nearly complete, whereas recordings from ON and OFF midget and small bistratified cells were less complete (Fig. 1).

Linear–nonlinear Poisson cascade model.

To summarize the light response properties of the major RGC types, the linear–nonlinear Poisson (LNP) cascade model was fitted to responses to white noise (Chichilnisky, 2001; Chichilnisky and Kalmar, 2002). At the retinal eccentricities recorded, the RF diameters of raphe RGCs are approximately half the diameter of those in the periphery (Figs. 1, 2). Therefore, to compare model parameters in the raphe and periphery, stimuli with a 22 μm pixel width were used in the raphe and stimuli with a 44 or 55 μm pixel width were used in the periphery. Both stimuli had a 16.74 ms refresh interval. Preparations recorded at different temperatures (32, 33, and 35°C) were all used, and there was a modest effect of temperature (i.e., lower temperatures resulted in slightly slower response time courses). However, the effects reported here were present at all temperatures (albeit with different magnitudes), and thus are not attributable to temperature differences.

First, the RF diameter of each cell was determined by fitting a 2D Gaussian to the spatial component of the STA (see above) and calculating the geometric mean of each SD parameter multiplied by 2. To summarize the response time course of each cell, the STA time course of each cell was normalized (L2) and fitted with a difference of cascades of low-pass filters of the form f(t)=p1(t/τ1)ne−n(t/τ1−1)−p2(t/τ2)ne−n(t/τ2−1), where f(t) is the time course, t is the time index, and n, p1, τ1, p2, and τ2 are free parameters (Chichilnisky and Kalmar, 2002). Then, the time of zero crossing of the time course was estimated by numerically solving for the root of the fitted function using the bisection method. Last, the contrast–response relationship for each cell was estimated from the STA, stimulus, and spike train (Chichilnisky, 2001; Chichilnisky and Kalmar, 2002) and fitted with a smooth function of the form f(x)=αL(βx+γ), where x is the generator signal, L is the logistic cumulative distribution function, and α,β, and γ are free parameters. The nonlinearity index for each cell was defined as the logarithm of the ratio of f′(x+σ) and f′(x−σ), where f′(x) is the first derivative of the contrast–response relation at x (computed numerically using the central difference method) and xσ is the 1 SD generator signal value estimated from the distribution of all generator signal values for that cell during the recording.

Analysis of spike train autocorrelation functions.

The autocorrelation function of the spike trains from white noise visual stimulation were computed for each parasol and midget cell. Briefly, spike times were binned with 1 ms precision and a histogram of pairwise spike time differences (range, 1–500 ms) was calculated and then the spike counts within each bin were then converted to a spike rate. To determine differences in firing between raphe and peripheral parasol cells, the autocorrelation was normalized to unit sum and the sum of the first 25 elements (1–25 ms) for ON parasol cells and the first 100 elements (1–100 ms) for OFF parasol cells was calculated. These quantities provide estimates of the probability that, after the occurrence of a spike, a second spike occurs within 25 and 100 ms, respectively.

Analysis of spike train cross-correlation functions.

The cross-correlation functions of spike trains obtained during white noise visual stimulation were computed between all homotypic pairs of ON parasol and OFF parasol cells. Briefly, spike times were binned with 1 ms precision and a histogram of pairwise spike time differences ranging from −1000 to 1000 ms (0 ms inclusive) was calculated, and then the spike counts within each bin were converted to a spike rate. To determine differences in correlated firing kinetics, nearest-neighbor homotypic pairs of ON parasol and OFF parasol cells were identified and an estimate of the baseline spike rate was determined by calculating the median value of the spike rate from −1000 to −995 ms time offset and from 995 to 1000 ms time offset. Next, the cross-correlation was smoothed with a Gaussian kernel (σ = 5 ms), the noise estimate was subtracted, and the resulting signal was upsampled 10×. To determine the width of the cross-correlation, the full-width at half-maximum (FWHM) was computed.

Analysis of receptive field overlap.

To determine the RF overlap between nearest-neighbor homotypic pairs of ON and OFF parasol cells, an estimate of the RF of each RGC was generated from its 2D Gaussian fit and the correlation coefficient was calculated between the RF of each cell and the RF of its nearest homotypic neighbor.

Analysis of electrical images.

Electrical images (EIs) were computed for each cell by averaging the voltage signal on the array during the time of each spike (Litke et al., 2004; Petrusca et al., 2007). Features of spike waveforms obtained from the EI were then used to classify them into distinct cellular compartments (Litke et al., 2004; Müller et al., 2015). In general, somatic and dendritic waveforms exhibit biphasic structure, whereas axonal waveforms exhibit triphasic structure, reflecting biophysical differences between compartments. Based on these known differences, a heuristic-based automated method (Madugula et al., 2022b) was applied to label each electrode as either “dendritic,” “somatic,” “axonal,” or “mixed” (a superposition of either somatic and dendritic or somatic and axonal waveforms; Madugula et al., 2022b). First, waveforms with a larger positive than negative phase were labeled as dendritic. Second, the waveforms at remaining electrodes were identified as somatic, mixed, or axonal, based on the ratio of the first positive peak to the second positive peak. Thresholds for each respective compartment were determined by two empirically observed inflection points (0.05 and 1.6, respectively) in the pooled distribution of ratios obtained from many cells and electrodes in peripheral recordings.

To determine spike amplitudes within each compartment, for each cell, electrodes recording somatic and axonal signals were determined and the maximum absolute value of the peak negative voltage on each electrode was calculated. Then, the maximum voltage over somatic and axonal electrodes was computed. Each preparation in this analysis was acquired using the electrode array with 30 μm interelectrode spacing to ensure consistent estimates of spike amplitudes (see Fig. 4B).

Axon spike conduction velocity was estimated from the electrical images of ON and OFF parasol and midget cells, as described previously (Zaidi et al., 2022; Madugula et al., 2022a). Briefly, axonal electrodes for each cell were identified (see above) and pairwise velocity estimates were computed by dividing the distance between each electrode by the time difference in the negative peak amplitude, for each pair of axonal electrodes. Then, a velocity estimate was computed by averaging the estimate from each pair of electrodes weighted by the product of the peak amplitudes.

Analysis of electrical features for distinguishing cell types.

Axon spike conduction velocities and the spike train autocorrelation function were used together to distinguish ON and OFF parasol and ON and OFF midget cells, to demonstrate the effectiveness of cell type classification based on intrinsic electrical features alone (Richard et al., 2016; Zaidi et al., 2022). Each preparation in this analysis (see Fig. 4C–G) was acquired using the electrode array with 30 μm interelectrode spacing to ensure consistent estimates of spike conduction velocities. To quantify the degree of separability between parasol and midget cell spike conduction velocities, a 1D support vector machine was fitted to the spike conduction velocity–cell type label pairs (L2-penalty with a linear kernel). The learned optimally separating hyperplane (a point in the 1D space) was used to quantify the classification accuracy on the training data. The goal of this analysis was merely to determine the degree of separability in the data, so cross-validation with held-out data was not performed.

Parasol and midget cells were highly distinguishable using spike conduction velocity (see Fig. 4D,F). To determine whether ON and OFF types within the parasol and midget cell classes could be separated by features of the autocorrelation function, the autocorrelation function was calculated (see above), and the resulting spike count vector was normalized by its L2 norm. Then, principal components analysis (PCA) was performed on a data matrix composed of either ON and OFF parasol or ON and OFF midget cell autocorrelation functions, and projections onto the first two principal components were analyzed (see Fig. 4E,G, insets). To quantify the degree of separability between ON and OFF (parasol or midget) spike train autocorrelation functions, a 2D support vector machine was fitted to the projections onto the first two principal components along with the corresponding cell type labels (L2 penalty with a linear kernel). The learned optimally separating hyperplane (a 1D line in the 2D space) was used to quantify the classification accuracy on the training data. As above, cross-validation with held-out data was not tested.

Responses to electrical stimulation.

To identify the responses to electrical stimulation of each recorded RGC, a current pulse was delivered to each electrode on the array, in a random sequence across electrodes, while recording RGC responses on all electrodes simultaneously (Jepson et al., 2013; Grosberg et al., 2017). The pulse was triphasic, positive first, and charge balanced (50 μs/phase; relative ratios, 2:–3:1). The current amplitude range tested (second phase) was 0.1–4.1 μA (∼40 amplitudes, logarithmic scale), with each amplitude repeated 25 times. In all electrical stimulation experiments, the electrode array with 30 μm interelectrode spacing was used.

An automated spike-sorting method was used to identify RGC responses to electrical stimulation. First, for each stimulating electrode, RGCs were identified as candidates for electrically evoked activity if their EIs (obtained during visual stimulation) exhibited peak amplitudes >0.5 times the SD of the recorded electrical noise on that electrode. Each RGC was then assigned a set of significant electrodes from its EI recording signals with power greater than two times the power of electrode noise on each electrode to be used to determine responses. Cells that did not have at least one electrode with signal power greater than two times the power of the electrical noise threshold were excluded from analysis because their signals could not be distinguished from noise. Then, for each current amplitude, across all RGCs and their respective significant electrodes, voltage traces recorded in the 3 ms period following stimulus application were grouped using unsupervised clustering. For each pair of clusters, signals were subtracted, aligned, and averaged to obtain differences; typically, these difference signals reflect the firing of one or more cells in one cluster but not the other. These residuals were iteratively compared with each of the recorded spike waveforms obtained from the EIs of all RGCs and their respective significant electrodes to identify the set of RGCs contributing to the responses in each trial for each amplitude.

This algorithm produces response probabilities as a function of current level, for each RGC and stimulating electrode, which is in general monotonically increasing (Hottowy et al., 2012; Jepson et al., 2013). The relationship between current amplitude and response probability for each RGC and stimulating electrode was then fitted by a sigmoidal curve of the form p(a)=1e−k(a−b), where p(a) is the probability of activation, a is the current amplitude, and k and b are free parameters, using maximum likelihood estimation. The stimulation threshold (i.e., the current amplitude required to elicit a response probability of 0.5) was extracted from the fit (Jepson et al., 2013).

Analysis of activation thresholds.

To determine activation the threshold within each compartment, for each cell, electrodes recording somatic and axonal signals were determined (see above) and the activation threshold for each electrode was extracted from the sigmoidal fit. Then, the minimum activation threshold was computed over somatic and axonal electrodes, for each cell.

Identification of axon bundle activation thresholds.

Axon bundle activation thresholds on each electrode were determined by an automated method based on a previously described algorithm (Tandon et al., 2021). The algorithm was modified to avoid bias resulting from differences in array geometries (the present work used a smaller, hexagonal array, whereas the algorithm described in the study by Tandon et al., 2021, was developed using a larger, rectangular array) and axon spike amplitude differences between central and peripheral RGCs (see Fig. 4B). For each preparation, a threshold voltage was first determined to identify electrodes that recorded significant axonal signals in response to electrical stimulation, as follows. For each RGC recorded during white noise visual stimulation, the electrodes recording axonal signals were identified as described above, and the average axonal spike amplitude was determined. The median axonal spike amplitude across all recorded RGCs was computed and was taken to be the threshold voltage. Next, to determine the axon bundle activation threshold, for each stimulus current applied, electrodes were first identified as either recording a signal or not, depending on whether the recorded signal was above the threshold voltage. Activity on the array was identified as an axon bundle activation event when the collection of electrodes recording signals formed a contiguous path reaching at least two nonadjacent edges of the electrode array. The bundle activation threshold was then defined as the minimum current level at which an axon bundle activation event was evoked. Electrodes near the border of the array (outermost two rings of electrodes) were excluded from analysis because their proximity to the edge precludes the ability to unambiguously distinguish RGC activity from axon bundle activity.

Analysis of selectivity.

To summarize the selectivity of electrical stimulation, a selectivity metric was computed for each ON and OFF parasol cell. First, at each stimulating electrode, the response probability at each current level below axon bundle activation threshold was determined for the target cell (p_target) and every other nontarget parasol cell (p_nontarget). Then, the maximum value over current amplitudes of the quantity (p_target) (1 – p_nontarget) was determined. This was repeated for each nontarget cell, and the minimum value of (p_target) (1 – p_nontarget) over nontarget cells was determined at this stimulating electrode, to provide a worst case estimate of stimulation selectivity for this cell, at this stimulating electrode. This was repeated for each stimulating electrode, and the maximum value of (p_target) (1 – p_nontarget) over stimulating electrodes was defined as the selectivity index for the cell.

Inference of visual perception.

To estimate the quality of restored vision that could be achieved with the measured responses to electrical stimulation, a recently developed stimulation algorithm (Shah et al., 2019b, 2022) was applied in a data-driven simulation. This method assumes that the brain reconstructs visual stimuli from RGC spike trains linearly and optimally, and that responses from individual RGCs that occur within a single “integration time window” of downstream visual processing (e.g., tens of milliseconds) are summed linearly to produce perception.

In this framework, an optimal linear reconstruction filter is associated with each RGC that quantifies its contribution to the image reconstruction each time it fires a spike. To compute the optimal reconstruction filter for each cell, parameters from the LNP model were used to generate responses to white noise images displayed for 100 ms. Specifically, the spike rate was simulated by taking an inner product between the STA and the stimulus and rectifying the result. The STAs were scaled such that the average parasol cell spike rate matched that observed experimentally (approximately three spikes on average per 100 ms flash in the peripheral retina). The stimuli are denoted by the matrix S∈ℝimages x pixels, and the responses in a matrix R∈ℝimages x cells. A reconstruction matrix, W∈ℝcells x pixels, was computed using linear least-squares regression, via the normal equations: W=(RTR)−1RTS. In this framework, the brain reconstructs an estimate of the true stimulus S from the RGC responses R by the linear operation Ŝ=RW (Warland et al., 1997). A single set of reconstruction filters was learned for a set of white noise images with varying pixel sizes (352, 220, 176, 110, 88, and 55 μm), using 10,000 training images each (60,000 images total) to test performance at a range of spatial frequencies. To produce linear reconstruction filters for naturalistic stimuli (Brackbill et al., 2020), this procedure was repeated with 60,000 grayscaled naturalistic images from the ImageNet database (Fei-Fei et al., 2009).

To determine the reconstruction that could be achieved with the responses to electrical stimulation, RGC response probabilities for each cell–electrode pair, at current levels below the threshold for axon bundle activation, were assembled in a matrix, D∈ℝcells x (electrodes x amplitudes), where the columns are distinct electrical stimuli and the rows are the responses of each recorded cell. The electrical stimulus delivered at each time step was chosen such that the change in the image reconstruction due to the evoked responses of RGCs would maximally reduce the expected mean-squared error (MSE) between the reconstructed image and the target image (Shah et al., 2019b, 2022). The algorithm sequentially chooses stimuli greedily at each time step to minimize the expected MSE until no further reduction in error occurs. No two stimuli resulting in a cellular response probability of >0.04 were delivered within 5 ms of each other to avoid the refractory period, and each stimulation pulse was assumed to be 150 μs in duration, to be consistent with the experimental data and hardware configuration. This method was applied to a total of 90 white noise images (15 for each pixel size) and 100 naturalistic stimuli from the ImageNet database (Fei-Fei et al., 2009), all separate from the images used for training the linear reconstruction filters.

To approximate the optimal reconstruction, or the reconstruction that could be achieved with perfect control over the activity of each RGC, a convex solver was used to solve for the optimal non-negative responses under squared loss between the reconstruction and the true stimulus. This response vector was then discretized by rounding each floating point response to the nearest integer value. Then, for each stimulus, the resulting response vector was applied to the reconstruction matrix to produce the approximate optimal reconstruction.

Convolutional autoencoder.

To enhance reconstructions of naturalistic images, a previously described convolutional autoencoder (CAE; Parthasarathy et al., 2017) was applied to linear reconstructions. The CAE was developed as follows. First, linear reconstruction filters were learned for each cell using least-squares regression (see above) on a set of 60,000 training images from ImageNet (Fei-Fei et al., 2009). Then, linearly reconstructed images were generated for each image in the training set using the LNP model generated responses and reconstruction filters. This linear reconstruction was then passed through an eight-layer convolutional neural network with four downsampling layers and four upsampling layers. The Adam optimizer (Kingma and Ba, 2014) with MSE loss and a learning rate of 0.0004 was used for model training. A validation set of 1000 images from ImageNet unseen by both the linear model and the CAE was used to evaluate model accuracy. Twenty training epochs with a batch size of 32 were used, after which validation error did not substantially decrease further.

Analysis of image reconstructions.

To quantify the quality of the white noise image reconstructions, two quantities were calculated. First, the set of all relevant pixels was determined for each preparation, indicating the regions covered by the RFs of all ON and OFF parasol cells. Then, the normalized MSE (NMSE) was determined by calculating the squared error between the reconstructed and original images over this set of relevant pixels and normalizing by the squared L2-norm of the original image within this set of pixels.

Next, the fraction of reconstructed pixels with polarities opposite that of the original stimulus was calculated. A subset of the aforementioned set of relevant pixels was determined, indicating regions covered by the RFs of ON and OFF parasol cells chosen to be stimulated by the algorithm. Then, in each reconstructed image, the polarity of each pixel in the reconstruction was compared with that in the original stimulus and the fraction of pixels with opposing polarities was calculated. For the naturalistic image reconstructions, in addition to NMSE, SSIM, a perceptual similarity metric (Wang et al., 2004; Wu et al., 2022), was used to quantify image quality in a way that takes image structure into account. Note that in some of the peripheral preparations (see Fig. 7B), SSIM was slightly larger for empirical naturalistic image reconstructions compared with optimal reconstructions after (but not before) applying the CAE. This is possible because the objective function for optimization in training the CAE was MSE, rather than SSIM. Separately, unlike in the raphe, the CAE had little effect on image reconstructions in some peripheral preparations because their respective linear reconstructions were coarse, because of the relative low density of RGCs.

Simulated midget cell receptive fields and reconstruction filters.

To determine the extent of midget cell activation and the noise corruption of reconstructions that would occur when targeting only parasol cells, a raphe preparation with partial recordings of midget cells was used to simulate the electrical stimulation and visual responses of complete populations of ON and OFF midget cells (see Fig. 8A–C). First, RF centers were generated from a hexagonal lattice. Ideally, the spacing of the lattice would be informed by the typical neighbor distance distribution of the midget cells. However, the midget cell populations from each raphe preparation were too sparse to reliably estimate the typical neighbor distance (Fig. 1A; see also Fig. 7A). Since the typical neighbor distance is proportional to the typical RF diameter (Gauthier et al., 2009b), the ON (or OFF) midget cell typical neighbor distance was therefore estimated by computing the median ON (or OFF) parasol typical neighbor distance and scaling this quantity by the ratio of median ON (or OFF) midget cell receptive field diameter to the median ON (or OFF) parasol cell receptive field diameter. This value was used to set the spacing of the lattice. The lattice locations were then jittered by adding Gaussian noise to each x and y coordinate independently, with mean 0 and SD equal to 10% of the original lattice spacing. Next, using the partial populations of ON or OFF midget cells, kernel density estimates were fitted to distributions of the larger of the two (major) SD parameters from the Gaussian fit to the RFs, the ratio of the major SD to minor SD parameters, and the tilt of the receptive fields, across all ON or OFF midget cells. To avoid extreme values, the kernel density estimate was fitted to the subset of the data within the 30th and 70th percentile of each distribution. Receptive fields were then generated by randomly sampling SD parameters and tilts from the fitted kernel density estimates, resulting in simulated receptive field mosaics (see Fig. 8B).

Next, STAs were generated using the Gaussian parameters from the simulated receptive fields. As with the parasol cells above, the STAs were scaled such that simulated responses to white noise resulted in, on average, three spikes within the simulated ON and OFF midget cell population, consistent with empirical results from parasol and midget cells in the periphery during flashed white noise visual stimulation (data not shown). Using least squares, a set of white noise linear reconstruction filters was learned using 60,000 training images (10,000 from each pixel size), and a set of naturalistic image reconstruction filters was learned using 60,000 training images from the ImageNet database (Fei-Fei et al., 2009).

Simulated midget cell responses to electrical stimulation.

To estimate responses to electrical stimulation, recordings during single-electrode stimulation of the partial raphe midget cell populations were analyzed similarly to the parasol cells (see above) by determining response probability as a function of current amplitude and fitting a sigmoid to the resulting relationship, for each cell-electrode pair. For each cell–electrode pair, the bivariate relationship between response probability and soma distance from electrode and current amplitude was determined. This was repeated for each cell–electrode pair, and the data were pooled over all cell–electrode pairs. From the pooled data, the cell–electrode distances and current amplitudes were binned and the average response probability within each bin was calculated, providing a lookup table for a response probability, given a distance from the stimulating electrode and current amplitude (see Fig. 8C).

Estimation of midget cell noise.

The electrical stimuli provided to the parasol cells for the white noise reconstruction simulation (see Fig. 6A,B) were then used to generate expected responses for each simulated midget cell. Specifically, for each electrical stimulus and cell, the distance from the electrode and the current amplitude were used to determine the response probability (see Fig. 8C). Over the entire stimulation sequence, the expected response for each midget cell would therefore be the sum of the probabilities from each stimulus. The resulting expected response vector was then passed as input to the linear decoder, yielding the midget cell reconstruction. Finally, this resulting reconstruction was summed with the parasol reconstruction to obtain an estimate of the expected perception when midget cells are ignored during parasol cell stimulation (see Fig. 8D).

Early stopping criterion.

To reduce the magnitude of midget cell activation and consequent noise added to white noise reconstructions, an early stopping procedure was developed based on an estimate of the midget cell noise. In its original form, the stimulation algorithm used here (Shah et al., 2019b, 2022) continues to provide electrical stimuli greedily until no further reduction in expected MSE between the target stimulus and the expected parasol cell reconstruction can be achieved. In the modified procedure, the MSE between the target image and the reconstruction with the midget noise added is calculated at each electrical stimulus, and the “early stopping point” is defined as the stimulation number at which the global minimum in MSE is achieved. The exact optimal early stopping point is unknown, and its estimation is directly affected by the estimate of the midget cell population, which is based on the true population, but is simulated and stochastic. To minimize potential bias in the early stopping point, 10 distinct midget cell populations were generated as above using 10 different seeds for the receptive field and electrical stimulation response generation. Then, the early stopping point for each midget cell population was determined for each of the 90 white noise stimuli, and the median stopping point across the midget cell populations was calculated for each stimulus. Then, a distinct midget cell population was used for evaluation (see Fig. 8B–D,F) using the early stopping points determined from the separate “training” midget cell populations.

Analysis of noise contributed by midget cells.

To quantify the extent of image corruption by the activation of midget cells, fractional error was calculated as error_{parasol+midget}/error_parasol, where error_{parasol+midget} is the error of the reconstruction with midget cell noise added to the parasol cell signals and error_parasol is the error of the reconstruction from parasol cells only. For white noise images, error was either the root MSE (RMSE) or the fraction of incorrect pixels. For naturalistic images, error was either RMSE or 1 - SSIM.

Experimental design and statistical analyses.

In nearly all statistical tests, a two-tailed Mann–Whitney U test was performed using the statistics package in SciPy (Virtanen et al., 2020). In the statistical test comparing the fraction of RGCs that could be activated with a probability of 0.5 without activating axon bundles (see Fig. 5B), a custom resampling procedure was implemented. All p values ≥ 1 × 10⁻¹⁰ are reported as exact values; all others are reported as p < 1 × 10⁻¹⁰.