Abstract
High-fidelity electronic implants can in principle restore the function of neural circuits by precisely activating neurons via extracellular stimulation. However, direct characterization of the individual electrical sensitivity of a large population of target neurons, to precisely control their activity, can be difficult or impossible. A potential solution is to leverage biophysical principles to infer sensitivity to electrical stimulation from features of spontaneous electrical activity, which can be recorded relatively easily. Here, this approach is developed and its potential value for vision restoration is tested quantitatively using large-scale multielectrode stimulation and recording from retinal ganglion cells (RGCs) of male and female macaque monkeys ex vivo. Electrodes recording larger spikes from a given cell exhibited lower stimulation thresholds across cell types, retinas, and eccentricities, with systematic and distinct trends for somas and axons. Thresholds for somatic stimulation increased with distance from the axon initial segment. The dependence of spike probability on injected current was inversely related to threshold, and was substantially steeper for axonal than somatic compartments, which could be identified by their recorded electrical signatures. Dendritic stimulation was largely ineffective for eliciting spikes. These trends were quantitatively reproduced with biophysical simulations. Results from human RGCs were broadly similar. The inference of stimulation sensitivity from recorded electrical features was tested in a data-driven simulation of visual reconstruction, revealing that the approach could significantly improve the function of future high-fidelity retinal implants.
SIGNIFICANCE STATEMENT This study demonstrates that individual in situ primate retinal ganglion cells of different types respond to artificially generated, external electrical fields in a systematic manner, in accordance with theoretical predictions, that allows for prediction of electrical stimulus sensitivity from recorded spontaneous activity. It also provides evidence that such an approach could be immensely helpful in the calibration of clinical retinal implants.
- biophysics
- macular degeneration
- multielectrode array
- primate
- retinal electrophysiology
- retinal ganglion cells
Introduction
Restoration of sensory capabilities using neural interfaces of the future will require evoking specific, naturalistic patterns of neural activity in large collections of neurons. This raises a crucial technical challenge: determining which electrodes in an electronic implant activate which cells and to what degree. In practice, with many electrodes and large neural populations, complete measurements of this kind are difficult or even impractical. Thus, any easily measured properties of neurons that could be harnessed to reveal features of their electrical sensitivity could significantly enhance the capabilities of neural interfaces.
A potentially impactful approach would be to leverage knowledge of neuronal biophysics to infer how cells will respond to electrical stimulation, using only measurements of spontaneous electrical activity. Spontaneous activity can be recorded rapidly in parallel on many electrodes, and is not subject to electrical artifacts, and therefore can be analyzed relatively easily. For example, electrodes closer to a particular cell tend to record larger spikes, and correspondingly also require less current to evoke spikes. Thus, spike amplitudes recorded during spontaneous activity could potentially be used to infer electrical stimulation thresholds for every electrode and cell. This kind of inference has been suggested in theoretical and modeling studies, and experimental work in tissue culture (Fohlmeister et al., 1990; Boinagrov et al., 2010; Tsai et al., 2012; Loizos et al., 2014; Radivojevic et al., 2016; Esler et al., 2018a). However, no experimental tests of these ideas have been performed in functioning neural circuits. Thus, it remains unclear whether and how the electrical sensitivity of neurons in intact tissue can be reliably inferred from electrical recordings and used effectively in a neural implant.
Here we propose an approach for electrical sensitivity inference in the context of a future high-resolution retinal implant, and empirically test its potential impact on vision restoration. We used large-scale multielectrode recording and stimulation from hundreds of retinal ganglion cells (RGCs) in the macaque retina ex vivo as a laboratory prototype, and probed the relationship between recorded spikes and sensitivity to electrical stimulation. Electrical sensitivity was systematically related to features of spiking activity, in a consistent manner across hundreds of electrodes and many recordings. This dependence exhibited different properties based on the proximity of each electrode to the soma, axon, and dendrites of the cell, and on the cell type, and was in quantitative agreement with predictions from biophysical models. We show that these trends can be leveraged to infer the responses of individual RGCs to electrical stimulation. Finally, we exploit these inferences to identify optimal electrical stimulation sequences for vision restoration, revealing the potential utility for a future high-fidelity retinal implant.
Materials and Methods
Experimental setup
Custom 512-electrode and 519-electrode systems (Hottowy et al., 2008, 2012; Grosberg et al., 2017) were used to stimulate and record peripheral and central RGCs, respectively, in isolated rhesus macaque monkey retinas (Macaca mulatta) and human retinas. Human eyes were obtained from 3 brain-dead donors (29-year-old Hispanic male, 27-year-old Hispanic male, 47-year-old white female) through Donor Network West. Macaque retinas were obtained from terminally anesthetized animals killed during the course of research performed by other laboratories, in accordance with institutional and national guidelines and regulations. Ocular hemisection was performed in ambient indoor lighting following enucleation, the vitreous was removed, and the posterior portion of the eye was stored in warm, oxygenated, bicarbonate buffered Ames' solution (Sigma-Aldrich) in darkness. Patches of retina ∼3 mm on a side were isolated under infrared light, placed RGC side down on the multielectrode array (MEA), and superfused with Ames solution at 35°C. Electrodes were 8-15 μm in diameter and were electroplated with platinum. The MEAs consisted either of 512 electrodes arranged in a 16 × 32 isosceles triangular lattice with 60 μm spacing, or of 519 electrodes arranged in a 16 × 32 isosceles triangular lattice with 30 μm spacing (Litke et al., 2004). Voltage recordings were bandpass filtered between 43 and 5000 Hz and sampled at 20 kHz. Spikes in the voltage recordings from individual RGCs evoked by light stimulation (which produced no electrical artifacts) were identified and sorted using previously described spike sorting techniques (Litke et al., 2004).
Visual stimulation and cell type classification
To identify the type of each recorded cell, as well as the location and shape of the visual receptive field, the retina was visually stimulated with a dynamic, binary white noise stimulus, and the spike-triggered average (STA) stimulus was computed for each RGC (Chichilnisky, 2001; Chichilnisky and Kalmar, 2002). The STA summarizes the spatial, temporal, and chromatic structure of the light response. The spatial receptive fields and time courses obtained from the STA were used to identify the distinct cell types, as described previously (Chichilnisky and Kalmar, 2002; Field et al., 2007; Rhoades et al., 2019).
Electrical image (EI) computation
An electrical signature for each neuron on the array was obtained by computing the EI (Litke et al., 2004). The EI of each cell (e.g., Fig. 1A) represents the average spatiotemporal pattern of voltage deflections produced on each electrode of the array during a spike from a given cell. A minority (∼3%) of spike clusters identified during spike sorting erroneously merged two different cells, and were excluded from analysis based on visual inspection. The remaining cells were used for subsequent analyses if their EIs featured spike amplitudes greater than twice the SD of the recorded voltage trace (i.e., the electrical recording noise) at a minimum of three recording electrodes.
Identifying cellular compartments and landmarks
As observed previously (Litke et al., 2004; Müller et al., 2015), the biophysical properties of somas, axons, and dendrites cause the spike waveforms recorded on MEA electrodes overlying each subcellular compartment to have distinct shapes (Fig. 1A, top row). Compartments were identified at each EI electrode in two steps. First, waveforms with a larger positive phase than negative phase were labeled as dendrites. 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.06 and 1.32, respectively) in the cumulative distribution of ratios, across tens of thousands of spike waveforms obtained from thousands of recorded parasol and midget cells.
The identification of the cellular compartment at each electrode allowed for the estimation of two morphologic landmarks for each RGC used in several analyses. The location of each RGC soma was estimated by taking the centroid of all identified somatic electrodes, weighted by the negative peak amplitude recorded on each electrode. The orientation of the axon was determined by computing a vector from the soma centroid to the nearest axonal electrode >180 μm away. The axonal initial segment (AIS) location was assumed to be 25 μm along this planar vector from the soma centroid, based on estimates from previous studies while taking into account the bend in the axon as it exits the RGC layer and enters the nerve fiber layer (Maturana et al., 2016; Raghuram et al., 2019). Data analysis using AIS locations of 13 μm (Sekirnjak et al., 2008) and 40 μm (assuming no axon bend) were also tested, and did not affect the results significantly (with a maximum 5.2% reduction in threshold prediction accuracy at 40 μm). The axon midline was estimated by conducting amplitude-weighted fitting of a third-degree polynomial spline to all identified axonal electrodes.
Electrical stimulation and spike sorting
Electrical stimuli consisted of triphasic current waveforms consisting of a cathodal phase flanked by two anodal charge-balancing pulses, with relative amplitudes 2:−3:1. Each phase was 50 μs in duration. The shape of the triphasic pulse was chosen to minimize the recorded electrical stimulation artifact (Hottowy et al., 2012). At each of the 512 or 219 electrodes on the 60 or 30 μm MEA, respectively, triphasic pulses with 39 stimulating-phase amplitudes increasing from 0.1 to 4.1 μA in 10% increments were delivered 25-50 times each. Spikes were elicited with submillisecond temporal precision (Jepson et al., 2013) by directly depolarizing the cell, as shown by synaptic transmission blockade (Sekirnjak et al., 2006). A semiautomated method was used to subtract electrical artifacts from the raw post-stimulation data and assign spikes to cells using waveform templates derived from their EIs, as described previously (Jepson et al., 2013; Mena et al., 2017), supplemented by another, previously described semiautomated clustering spike-sorting method well-suited for analyzing electrodes with small peak waveforms (Madugula et al., 2022). To enable accurate spike assignment, at each electrode, only cells with recorded spike amplitudes of at least 25 µV were analyzed because spikes with lower amplitudes were difficult to distinguish from background electrode noise. For macaque data, stimulation amplitudes above axon bundle activation thresholds, which were estimated at each electrode using a previously described method for the 60 μm pitch MEA (Tandon et al., 2021) or a custom algorithm for the 30 μm pitch MEA (see below), were not analyzed to avoid potential off-array RGC activation. For human data, on the other hand, stimulation amplitudes above axon bundle activation thresholds were included for analysis to compensate for the limited tissue availability. For each stimulus amplitude, the evoked spike probability was computed across repeats. The dependence of spike probability on stimulation current amplitude was modeled by a sigmoidal relationship as follows:
Identification of axon bundle activation thresholds on the 30 μm MEA
Axon bundle activation thresholds on each electrode were determined by an automated method based on a previously described algorithm (Tandon et al., 2021), modified to avoid bias resulting from differences in array geometries, and from smaller axon spike amplitudes for central compared with peripheral RGCs. 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 activated or inactivated, 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 activated electrodes 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 (outer two rings of electrodes) were excluded from analysis because their proximity to the edge complicates the ability to unambiguously distinguish RGC activity from axon bundle activity.
Retinal preparation, cell, and activation curve selection
Of 53 total available peripheral preparations from 50 different macaque retinas, 23 were included for analysis. Additionally, of 11 total available central preparations from 11 different macaque retinas, 10 were included. Central retina preparations were obtained from the raphe region, at 2-4.5 mm eccentricity along the temporal horizontal meridian, and recorded using an MEA with 30 μm electrode spacing (see above). Finally, of 8 total available peripheral preparations from 4 different human retinas, 6 were included. In the chosen preparations, the excised segment of retina covered >80% of the electrodes on the array, and the population RGC firing rate exhibited variance <20% of the mean. Only somatic and axonal cell-electrode pairs recording spikes >25 µV, with activation curve thresholds less than both the bundle threshold on that electrode, which accounted for ∼80% of electrodes that were not able to be analyzed, and the upper current limit of the stimulation range (4 µV) were considered for analysis. Based on these criteria, 246 cell-electrode pairs from 59 ON and 43 OFF peripheral macaque parasol cells, 336 cell-electrode pairs from 185 ON and 70 OFF peripheral macaque midget cells, 1126 cell-electrode pairs from 103 ON and 309 OFF raphe macaque parasol cells, 127 cell-electrode pairs from 60 ON and 13 OFF human parasol cells, and 89 cell-electrode pairs from 28 ON and 23 OFF human midget cells were analyzed.
Derivation and fitting of inference relations
For the relationship between activation threshold and recorded spike amplitude, this inverse functional form was chosen based on theoretical predictions from the literature (Rattay and Wenger, 2010) that electrodes situated <50 μm away from neuron fibers exhibit a linear distance–threshold relationship, while electrodes >50 μm from neuron fibers exhibit a squared distance–threshold relationship. While it is almost certain that most of the MEA electrodes used in this study were <50 μm from the activated RGCs because of the tight juxtaposition of the MEA pressed down on the retinal preparation (Sekirnjak et al., 2008), a squared inverse function of the form
The linear inverse functional form
Inverse functions for each activation curve parameter were fitted by maximizing the log-likelihood of the evoked spike probability data. For simulated data, however, the inverse function was fitted directly to the data using least squares regression. The coefficient of determination (r2) values for inverse curve fits were computed from residuals to a linear regression fit to log-transformed spike amplitude versus activation threshold, or activation threshold versus slope data.
Biophysical simulations
RGC simulations were conducted using the NEURON software package (Hines and Carnevale, 2001; Finn et al., 2020) to create a model cell featuring five ion channel types, along with an extracellular stimulation mechanism, with an integration timestep of 0.001 ms. The properties of the model RGC were based on ion channel densities and compartment sizes adapted from voltage-clamp and morphologic measurements of rat ON-α RGCs (Fohlmeister et al., 1990, 2010; Raghuram et al., 2019). Slight modifications to the model geometry were made for primate peripheral parasol and midget RGCs, guided by previous measurements (Watanabe and Rodieck, 1989; Peterson and Dacey, 1998; Jeng et al., 2011; Grosberg et al., 2017; Kántor et al., 2018). Geometry was further slightly modified to produce the best match to the empirical data, which could be achieved within the bounds on cell diameter reported in the literature, and is listed for each cellular compartment along with fully open ion-channel conductance values in Table 1. The dendritic compartment morphology for both modeled cells was taken from Kántor et al. (2018) (note the dendritic compartment is not shown in Fig. 5A because of its complex geometry).
The extracellular environment was assumed to be an isotropic, uniform medium with a linear resistivity of 1 kOhm-cm (Hadjinicolaou et al., 2016; Kasi et al., 2019). Extracellular stimulation was conducted with a collection of simulated point-electrodes located at a retinal depth of 7 μm relative to the axon, at 60 μm intervals along the length of the soma and axon (see Fig. 5A). The spatial relationship between the electrodes and underlying modeled cell was varied to simulate various effective distances. Specifically, the electrode location was jittered along the axes parallel to the distal axon or sodium-channel band in steps of 5 μm, and along the perpendicular axes in steps of 2 μm. Simulated stimulus waveforms consisted of triphasic current pulses with shapes and amplitudes matching those used in ex vivo experiments (see above). The extracellular potential distribution generated by the stimulus was applied to each simulated cell compartment using the extracellular mechanism from NEURON (Carnevale and Hines, 2006). The mechanism calculates current pathways along and radial (membrane conductance: 1e10 S/cm2) and longitudinal axes (resistivity: 136.6 Ohm-cm) of the RGC membrane and solves for extracellular voltage simultaneously with the membrane voltage, which is then used to drive cell response using membrane dynamics (Fohlmeister et al., 1990, 2010; Raghuram et al., 2019). Poststimulation currents resulting from simulated action potential currents in the model cell were summed across each model element and used to compute simulated recording voltage at each stimulating electrode location. All simulations were conducted at 35°C.
The nondeterministic nature of extracellular stimulation putatively arises from stochastic membrane voltage fluctuations, caused by a combination of rapid ion channel state changes and synaptic inputs. To approximate these fluctuations, random voltage values drawn from a Gaussian distribution were added to the membrane voltage of model segments, at intervals of 0.15 ms. The noise was normalized to reflect the relative ion channel density differences between the distal axonal and AIS (sodium-channel band) segments: noise with an SD of 4 μV was injected at the AIS segments (F. Rieke, personal communication), versus noise with an SD of 1 µV injected at the distal axonal segments. The resulting noise was smoothened by a 1D Gaussian across five neighboring segments at a time along the cell, to avoid unrealistically sharp changes in membrane voltage between neighboring model elements.
Simulated image reconstruction
Simulation of perceived images resulting from application of stimulating current pulses, and selection of optimal electrical stimuli, were based on previous work (Shah et al., 2019) and consisted of several steps.
First, the contribution of RGC activity to visual perception was estimated in each of three retinal preparations, assuming that perception represents an optimal linear reconstruction of the stimulus from the RGC responses (Warland et al., 1997; Brackbill et al., 2020). Optimal linear reconstruction was examined by simulating the response of each cell to 10,000 12 × 6 random black-and-white noise training images, based on the spatiotemporal STA (see above), and then computing reconstruction filters from the simulated responses for each cell as follows. The simulated nonnegative response to each training image was estimated by taking the inner product of the STA and that image, and rectifying the result. Then, the reconstruction filter was obtained by performing least squares regression of the simulated responses against the set of training stimuli.
Next, for each of 15 randomly generated 12 × 6 black-and-white images, 40 electrical stimuli were chosen to (approximately) produce an appropriate neural response in each retinal preparation. Each stimulus consisted of a particular current amplitude on a particular electrode and produced spikes with probabilities from 0 to 1 in one or more RGCs, estimated using either measured, inferred, or fixed activation curves. The fixed curve was assumed using parameters of 1.2 μA for threshold and 12.1 μA−1 for slope, obtained by averaging across all measured responses to electrical stimulation in 53 preparations (see Fig. 6B–G). For each type of activation curve, stimuli were chosen so that the sum of their presumed contributions to perception, based on the linear reconstruction filters weighted by the expected response, minimized the normalized mean squared error (NMSE) between the stimulus and reconstructed image. Normalization was relative to the highest achievable mean-squared error across images and retinas, for ease of interpretation. The number of electrical stimuli (40) was empirically determined as the value beyond which, on average, reconstruction using inferred activation curves no longer reduced the NMSE.
To quantify the effectiveness of using inferred over fixed activation curves for image reconstruction, NMSE values of reconstructed images based on the three curve types were collected for each retinal preparation and target image. Then, NMSEs for inferred and measured activation curves were compared with NMSEs for fixed curves (see Fig. 6B). Across preparations and targets, systematic deviations from the x = y line, which corresponds to using the fixed curve, were captured by slopes obtained from linear regression. Finally, the overall value of performing reconstruction with inferred activation curves was given by |(sm – si)/(sm – 1)|, where sm is the slope obtained with measured activation curves and si is the slope obtained with inference.
Results
An experimental lab prototype of a future retinal implant was used to explore the relationship over space between recorded electrical features and electrical sensitivity of both macaque and human RGCs. Electrical recording and stimulation of RGCs in ex vivo isolated retinas was performed using a large-scale recording system (512 electrodes, 30 or 60 μm pitch; see Materials and Methods). After cell type identification by clustering responses to white noise visual stimulation, 498 ON and OFF parasol and midget cells from 36 macaque retinal preparations, and 48 cells of the same types from three human retinas, were probed with electrical stimulation, eliciting individual, precisely timed, directly evoked spikes (Sekirnjak et al., 2008; Jepson et al., 2013). Next, to obtain a spatial signature of spontaneous spikes produced by each cell, the EI was computed and used to identify cellular compartments (see Materials and Methods). Finally, to determine the extracellular activation characteristics of each cell, responses evoked by electrical stimulation at each electrode were used to identify the activation threshold (current amplitude that evoked spikes with 50% probability), and the map of electrical sensitivity over space was defined as the electrical receptive field (ERF) (Fig. 1).
Features of EIs correspond to the ERFs of ON and OFF parasol cells
Features of the EI, a relatively simple measurement, provided substantial information about the ERF, which captures fully the electrical excitability over space but is far more difficult to obtain. The relationship between the EI and ERF is illustrated in several aspects of the data from a single cell. First, only the electrodes that recorded large spike waveforms (Fig. 1, hollow circles) were able to evoke spiking (Fig. 1, filled circles) within the amplitude range tested (5-200 pC injected charge). However, many such electrodes did not activate the cell at current levels below axon bundle activation threshold (Fig. 1, light blue hollow circles; see Materials and Methods). Second, for the electrodes that did evoke a spike below bundle threshold (Fig. 1; see Materials and Methods), the EI amplitude at each electrode was correlated with the strength of the ERF: electrodes recording higher EI amplitudes (Fig. 1, circle size) also excited the cell more effectively (Fig. 1, circle color). The similarity between EIs and ERFs presumably reflects the fact that both spike amplitude and sensitivity to electrical stimulation are inversely related to the distance between the electrode and the cell (Rattay, 1987; Sekirnjak et al., 2008; Rattay et al., 2012). Third, somatic electrodes closer to the axon (Fig. 1B, electrode at Circle 2) required less current to activate RGCs than somatic electrodes further from the axon (Fig. 1B, electrode at Circle 1), even if the electrodes recorded spikes of similar amplitude (Fig. 1B, radii of Circles 1 and 2). Presumably this is because of a site of spike generation at the AIS. Fourth, electrical stimulation using dendritic electrodes typically did not elicit spiking within the tested stimulation amplitude range (Fig. 1A, top and bottom right). Finally, axonal electrodes exhibited steeper increases of firing probability in response to increasing stimulation current (Fig. 1 bottom, red vs blue sigmoids, slopes), and recorded spike waveforms with smaller amplitudes (Fig. 1 top, red vs blue waveforms, negative peak amplitudes), despite exhibiting activation thresholds in a range similar to those of somatic electrodes (0.2-2 μA). The latter observation suggests that the relationship between EIs and activation curves differs for somatic and axonal electrodes. In what follows, these trends in the relationship between the EI and ERF are explored across many cells and preparations of the macaque and human retina.
Electrical activation threshold is inversely related to recorded spike amplitude
The relationship between EI amplitude and ERF strength observed on axonal electrodes (Fig. 1, marker size and color) was consistent across many electrodes and cells in several retinas. First, recording and stimulation were tested on 61 electrodes overlying the axons of 53 ON and OFF peripheral parasol cells in one retina. Axonal thresholds exhibited a systematic inverse relationship to EI amplitudes (Fig. 2A,B, red markers). To determine whether this relationship holds across preparations, recording and stimulation with 152 axonal electrodes overlying 102 cells were compared across 23 recordings (Fig. 2C, see Materials and Methods). A consistent inverse relationship was observed (Fig. 2C, circle and square markers respectively), similar to the trend observed in a single retina (Fig. 2A,B, red points vs points in Fig. 2C). The relationship between axonal EI amplitudes and activation thresholds was fitted by a first-order inverse function based on theoretical studies of how electrode distance relates to activation threshold at the distances tested (Rattay and Wenger, 2010) as follows:
Similarly to axonal activation thresholds, somatic thresholds exhibited an inverse dependency on EI amplitude, as observed in 55 somatic electrodes overlying 53 ON and OFF parasol cells in two retinal preparations (Fig. 2A,B, blue markers). However, pooling the recording and stimulation features and inspecting the locations of 94 somatic electrodes from the same cells and retinal preparations used for pooled axonal electrode analysis (above) revealed that somatic activation thresholds depended on the distance of the stimulating electrode from the putative location of the AIS (Fig. 2D, marker colors; see Materials and Methods). Specifically, for a given RGC, somatic electrodes near the AIS (Fig. 2D, blue markers below dashed line) required less current to cause activation than electrodes farther away (Fig. 2D, red markers above dashed line), despite recording spikes with similar peak amplitudes. Quantitative characterization of this trend revealed that the fractional difference between measured somatic activation thresholds and estimates based solely on EI amplitude increased roughly linearly as a function of electrode distance from the AIS (r2 = 0.78). Together, these findings suggested the following relationship between spike amplitudes and activation thresholds near the soma:
Most electrodes near dendrites did not reliably elicit spikes (Fig. 1A, bottom, black trace). Across 453 cells, and 34 peripheral and central retinal preparations (see Materials and Methods), on average, only 9% of identified electrodes (12 of 133) recording significant spikes (see Materials and Methods) from each cell were identified as recording from dendrites. None of these electrodes was able to elicit spikes with probability >0.5 over the range of amplitudes tested, compared with 74% for all axonal electrodes and 93% for all somatic electrodes (see Materials and Methods). Thus, the dendritic data were insufficient to explore the EI amplitude-activation threshold relationship and were not further analyzed. At a minimum, these data imply that dendritic electrodes have substantially higher activation thresholds.
Activation curve slope is inversely related to activation threshold
A complete characterization of RGC responsivity to extracellular stimulation of varying strength entails determining not only the 50% activation threshold, but also the spiking probability over a range of stimulus amplitudes (i.e., a full activation curve, Fig. 3A). A simple possibility is that, for each cellular compartment, the activation curve has a stereotyped form across cells, and that this form scales with the injected current according to the electrical impedance between the electrode and the site of activation, which increases with distance. This possibility is corroborated by recent theoretical studies demonstrating that modeled axon fiber activation curve slope scales with increasingly distant electrodes (Rattay and Tanzer, 2022). In this case, normalizing experimentally measured activation curves by their 50% threshold should result in a function that is constant across cells and electrodes. Indeed, normalization revealed strikingly similar activation curves for all somatic electrodes, and for all axonal electrodes, but very different forms for the two compartments (246 electrodes, 102 cells, 23 preparations, Fig. 3C, blue vs red curves). Specifically, axonal activation curves featured a relatively steep “all-or-none” relationship, while somatic activation curves were shallower (mean of resampled Wilcoxon signed rank test: p = 2.3e-6).
The similarity of normalized activation curves within each cellular compartment indicates that the slope is inversely related to the threshold, as would be expected if the main source of variation is the electrical impedance (including the distance) between electrodes and cells. To characterize this relationship further, activation curve slopes were examined as a function of activation threshold, separately for 94 somatic and 152 axonal electrodes in 23 retinas (Fig. 3B,D). This relationship was fitted by an inverse function:
Comparison across cell types, retinal location, and species
The relationships observed between EI amplitude, electrode position, and activation threshold for peripheral macaque parasol cells originate from their broad biophysical characteristics, which could potentially generalize to other RGC types, to RGCs in the central retina, and to human RGCs. To test this hypothesis, electrical recording and stimulation of the somas and axons of peripheral midget cells, parasol cells in the central raphe region, and peripheral midget and parasol cells from human retinas were examined (for counts by type, see Materials and Methods, Fig. 4).
Peripheral macaque midget cells exhibited an inverse relation between recorded EI amplitudes and activation thresholds (Fig. 4A,D), for both somas and axons, which were characterized by a translated version of the relation obtained for peripheral parasol RGCs (Eq. 1; Fig. 4A,D, dashed vs dotted line, r2 = 0.83 for somas, 0.79 for axons). However, unlike with parasol cells, midget RGC somatic electrode thresholds were not obviously influenced by proximity to the AIS (r2 = 0.0001, Fig. 4A, marker colors; see Discussion).
Central parasol cell activation thresholds, measured using an MEA with 30 μm pitch (instead of 60 μm pitch used for all other experiments), also depended on recorded EI amplitudes, with inverse relationships that broadly followed with the fit to peripheral midget cell data (Fig. 4B,E, markers vs dashed curve, r2 = 0.72 for somas, 0.65 for axons). The similarity of electrical recording and activation features between central parasol cells and peripheral midget cells may arise from the similarity in cell size (Watanabe and Rodieck, 1989). Again, as for macaque peripheral midget RGCs, central parasol and midget RGC somatic electrode thresholds were not obviously influenced by proximity to the AIS (r2 = 0.01, Fig. 4B, marker colors).
Finally, to explore the possibility of applying the above findings in a clinical implant, the EI amplitudes and activation thresholds of a more limited collection of human peripheral parasol and midget cells were examined. The inverse relationship between human peripheral midget cell EI amplitudes and activation thresholds for somatic and axonal electrodes roughly coincided with the trend observed for peripheral macaque RGCs (Fig. 4C,F, triangular markers vs dashed curve fit to data in Fig. 4A,D, r2 = 0.67 for somas, 0.61 for axons). Parasol cell EI amplitudes and activation thresholds were broadly consistent with fits to macaque data but were more variable and thus did not provide a strong test (Fig. 4C,F, circular markers vs dotted curve fit to data in Fig. 2C,D, r2 = 0.37 for somas, 0.31 for axons, see Discussion).
The dependence of activation curve slope on activation threshold explored for macaque peripheral parasol cells above could also potentially generalize to different RGC types, eccentricities, and to the human retina. The macaque parasol cell inverse slope-threshold relationships for somatic and axonal electrodes (Fig. 3B,D, Eq. 3) provided a reasonable prediction of the midget cell, central retina, and human retina data (r2 = 0.72, 0.64, and 0.84, respectively, for somatic electrodes, 0.84, 0.87, and 0.86, respectively, for axonal electrodes; Fig. 4G–I, blue and red dashed lines vs corresponding markers).
Biophysical simulations approximately corroborate experimentally observed extracellular activation trends
Simulated extracellular stimulation of mammalian RGCs (Fohlmeister et al., 2010; Fohlmeister and Miller, 1997) reproduced the main empirically observed trends in electrical activation properties. Simulations were conducted using the NEURON software package (Carnevale and Hines, 2006) (see Materials and Methods) to implement a previously published electrical circuit model of a rat RGC, with soma and axon diameters modified to reflect primate parasol and midget cell geometry and to match the empirical data (Watanabe and Rodieck, 1989; Peterson and Dacey, 1998; Jeng et al., 2011; Grosberg et al., 2017; Kántor et al., 2018) (see Materials and Methods, Fig. 5A). The model cells were simulated in a uniform medium with isotropic resistivity. To simulate nondeterministic activation, spontaneous membrane voltage fluctuations (mimicking noise channel noise and synaptic inputs) were modeled as Gaussian noise (SD = 1 µV for axonal compartments, 4 µV for somatic compartments) in each model segment (see Materials and Methods). Stimulating electrodes were modeled as a row of point current sources with 60 μm spacing along the cell and axon, at a retinal depth 7 μm from that of the axon (Fig. 5A, black, blue, and red row of bars spanning the top). Voltages at each simulated electrode were obtained by summing the voltage fluctuations resulting from action potentials in the model segments (Gold et al., 2006).
Analysis of these simulations approximately matched several key findings in the recorded data. First, simulated electrodes located over the soma (Fig. 5A, blue rectangle) recorded spikes with 3-5 times the amplitude of electrodes overlying the axon (Fig. 5A, blue rectangle), but activated the model cell with similar current thresholds (Fig. 5C,D, blue vs red markers). Second, stimulation and recording over a range of electrode-cell distances (see Materials and Methods) resulted in inverse relationships between spike amplitude and activation threshold for both somatic and axonal electrodes, exhibiting significant overlap with the corresponding experimental data for both parasol and midget cells, albeit with slight differences in trends (Fig. 5C,D, red and blue circular and triangular markers vs dashed and solid lines). Third, electrodes recording somatic spike waveforms closer to the AIS had lower activation thresholds than electrodes farther from the AIS (Fig. 5B, circle colors). Fourth, simulated stimulation of dendrites did not produce activation within the tested current range (Fig. 5B, leftmost circle). Finally, stimulation with increasing current produced steeper activation curves at axonal electrodes than at somatic electrodes. The simulated somatic activation curve slope and axonal curve slope were within the range observed in the empirical data (Fig. 5E, blue and red curves). Together, these results suggest that the experimentally observed relationships between the electrical features of RGCs and their activation properties are mediated by well-understood biophysical mechanisms.
Activation curves can be accurately inferred from EIs
The systematic relationship between the EI and ERF probed above proved useful for inferring the parameters of ON and OFF parasol and midget cell activation curves using only recorded spikes. To test this, the previously fitted inverse relations were used to estimate activation thresholds, and these estimates were compared with measured thresholds, for the somas and axons of 102 parasol cells (Fig. 6A, example ON Parasol cell), across 246 electrodes from 23 different preparations (Fig. 6B,C,E,F). The accuracy of inferred thresholds was compared with that of a single, naive estimate of the activation threshold obtained by averaging measured activation curve parameters over all axonal and somatic electrodes across 53 retinal preparations (Fig. 6B,C,E,F, horizontal dashed lines, average threshold 1.2 μA, average slope 12.1 1/μA). Inferred ERFs were qualitatively similar to measured ERFs for individual cells (Fig. 6A left vs right, circle colors). Across axonal electrodes, 59% of activation thresholds were estimated within 25% of the measured value, compared with 33% obtained by averaging (Fig. 6E). Across somatic electrodes, 50% of inferred thresholds were within 25% of the measured value, compared with 28% obtained by averaging (Fig. 6B). Activation slopes inferred from these estimated activation thresholds were within 25% of the measured slope value for 58% of axonal electrodes, compared with 14% using the average slope (Fig. 6F), and for 67% of somatic electrodes, compared with 22% using the average slope (Fig. 6C). The r2 values for inferred versus measured somatic activation threshold, somatic activation slope, axonal activation threshold, and axonal activation slope were 0.64, 0.63, 0.59, and 0.69, respectively.
These inferred activation curves provided accurate estimates of cellular activation probabilities over the range of applied stimulus currents. This was determined by pooling data across 246 electrodes and comparing measured curves to the corresponding inferred curves (using the fitted inverse relations for threshold and slope inference, respectively), or to a single, naively estimated sigmoidal activation curve derived from average parameters (as above). The difference between each measured and estimated sigmoidal activation curve was estimated by computing the absolute difference of the integrated area between the curves. The separation between fitted and inferred activation curves was significantly smaller than the separation between fitted and averaged curves for somas and axons (Fig. 6D,G, Wilcoxon signed rank test between green and purple histograms: p = 0.002 for somas, p = 0.02 for axons). The mean separation between inferred and measured somatic activation curves was 24% smaller than between the fixed and measured curve; for axonal activation curves, the corresponding figure was 32% (Fig. 6D,G, vertical dashed lines).
Inferred activation curves can be used to inform stimulation choices for vision restoration
The inferred electrical responsivity of RGCs, including both the activation curve slope and threshold (Fig. 6), was useful for optimizing electrical stimulation over space using an algorithm designed for use in a high-resolution epiretinal implant (Shah et al., 2019). This was demonstrated by estimating the cumulative potential visual perception (estimated by image reconstruction from RGC spikes) (Warland et al., 1997; Brackbill et al., 2020), resulting from a collection of current pulses delivered by electrodes across the array (see Materials and Methods), separately in three different retinal preparations. In each preparation, activation of each ON or OFF parasol cell was assumed to contribute a particular image component to visual perception, corresponding to the optimal linear reconstruction filter obtained using responses to white-noise visual stimuli (see Materials and Methods). Spiking probabilities were computed using either measured responses to electrical stimulation (ground truth), inferred responses based purely on electrical recordings, or a fixed activation curve representing the average across 246 electrodes and 102 cells in 53 preparations. Each image was translated into a spatial pattern of stimulation, using 40 individual electrodes and current amplitudes selected to minimize the NMSE (see Materials and Methods) between the reconstructed and true image.
Qualitative comparison between reconstructed images revealed that choosing electrodes based on inferred activation curves resulted in substantially more accurate reconstruction than could be obtained by using the fixed curve, and approached the best achievable reconstruction obtained using fully measured activation curves (Fig. 7A, example reconstructions of two target images in a single preparation). Comparing the accuracies of 45 reconstructed targets across three retinal preparations revealed that choosing stimuli using inferred activation curves provided a 52% improvement in NMSE over using fixed curves (see Materials and Methods, Fig. 7B, dashed red, green, and black dashed line slopes), compared with a 94% improvement obtained with the measured curves. In other words, the inference method developed here produced approximately half of the gains that would be obtained with exhaustive electrical stimulation and recording in a given retina. These results suggest that inference of activation curves from EIs has the potential to produce more accurate visual perception with a future, closed-loop epiretinal implant.
Discussion
The present data reveal a systematic relationship between the electrical features of neurons and their responsivity to electrical stimulation over space, most notably an inverse relationship between spike amplitude and stimulation threshold, and an inverse relation between threshold and activation curve slope, for each cellular compartment. These trends were observed across cell types in human and macaque retina. They were also replicated in simulations with a biophysical model, indicating that they arise from well-understood properties of neurons. Finally, we showed that the observed relationships can be harnessed to accurately infer the extracellular activation properties of cells over space, and data-driven simulations reveal that this inference should be useful in guiding electrical stimulation choices to optimize prosthetic vision.
Several key features of electrical activation properties presented here are supported by previous studies. First, the inverse relationship between spike amplitude and activation threshold agrees with theoretical predictions, which suggest that the distance–threshold relationship is linear at electrode distances <∼50 μm (Rattay, 1987; Rattay and Wenger, 2010). The linear relationship in turn supports the observed inverse relationship between recorded spike amplitude and threshold. In addition, reports from previous experimental and theoretical studies indicate that the AIS, which is positioned between the soma and distal axon of RGCs, is highly excitable (Boiko et al., 2003; Fried et al., 2009; Rattay and Wenger, 2010; Radivojevic et al., 2016; Esler et al., 2018b; Werginz et al., 2020), consistent with the increase in activation threshold with distance from the AIS for electrodes near the soma. However, applying the inference relations presented here to other devices may require the collection of at least some electrical stimulation data because stimulation thresholds (but not recorded spike amplitudes) may vary based on electrode geometry (Lempka et al., 2011; Cao et al., 2015; Viswam et al., 2019). Regardless, simulations reveal that the present observations could be used to limit the number of required measurements (Shah et al., 2019).
In a notable previous study, inference of activation thresholds from the recorded features of electrically evoked spikes was performed in 14 cultured cortical neurons (Radivojevic et al., 2016). Inference was performed using multivariate regression on eight features of recorded spikes, revealing that the two critical features for inference were spike amplitude and time of spike onset. These two features approximately correspond to those used for threshold inference in the present work: spike onset times often identify the cellular compartment because of the stereotyped time course of action potential propagation.
Several limitations of the present work bear mentioning. First, the results were obtained from healthy primate retinas, but degenerate retinas would be more relevant for application to epiretinal implants. Experimental and computational studies of electrical stimulation using rat models of retinal degeneration indicate that recorded spike waveforms and stimulation thresholds do not significantly change with degeneration (Sekirnjak et al., 2009; Loizos et al., 2018). Earlier studies in mice suggesting that thresholds increase during degeneration were performed using long pulses and large electrodes (Suzuki et al., 2004; O'Hearn et al., 2006; Jensen and Rizzo, 2008) that likely produce indirect network-mediated activation of RGCs (Suzuki et al., 2004). It remains to be seen whether the relationship between recorded spikes and direct activation thresholds of the kind measured here is the same in diseased and healthy primate retina.
Second, the relationship between threshold and recorded EI amplitude (as well as distance from the AIS for somatic electrodes) exhibited significant variability, with an average r2 = 0.62 across the data from macaque (peripheral and central, parasol and midget cell, axonal and somatic). The observed variability could arise from a number of factors, not the least of which is the varying eccentricity of the retinal preparations considered in this study. However, despite the variation in threshold, data-driven simulations indicate that the observed trends can support clinically useful inference in a future implant (Figs. 6 and 7). The variability in data from the human retina was even greater; and although the trends were broadly consistent with trends observed in macaque data, the data did not provide a decisive test of the species similarity.
Third, for analysis of macaque RGC data, only activation thresholds below axon bundle threshold were included in this study, and estimating axon bundle thresholds currently requires electrical stimulation and collection of electrical response data (Grosberg et al., 2017; Tandon et al., 2021). In theory, it is perhaps possible that axon bundles could be removed during device implantation on the retina by incising the tissue on the distal edge, inducing Wallerian degeneration, but this approach has not been tested and may damage the surrounding tissue. Despite the complication of axon bundles, it was still possible to choose many stimulating electrodes and current levels that selectively activated target RGCs for image reconstruction (Grosberg et al., 2017) (Fig. 7).
Fourth, a single current pulse matching the underlying activation threshold of a given target cell delivered by a given electrode may simultaneously activate other RGCs, potentially leading to unwanted visual percepts. This can be addressed by careful selection of stimulating electrodes and currents that activate RGCs with separated response curves and, more importantly, by using inferred activation curve slopes to deliver subthreshold or suprathreshold current pulses to reduce the probability of unwanted activation. Indeed, 92% of stimulating currents chosen by the image reconstruction algorithm were either 20% more or 20% less than the inferred or measured target cell activation threshold (Fig. 7).
Finally, somatic activation thresholds for macaque peripheral midget cells and central parasol cells did not exhibit the inverse dependency on proximity to AIS observed for peripheral and central parasol cells. This could be because of the difficulty of accurately estimating the location of the AIS in small cells with an MEA featuring 60 μm electrode pitch, since the trend is clear in central parasol cells, which have EIs that are approximately the same size as peripheral midget cells. Conversely, it may be that there is a true difference in AIS sensitivity of parasol and midget cells, although this observation is not suggested by any previous findings.
Despite these limitations, inference of RGC activation properties purely from spontaneous activity using biophysically grounded relationships has the potential to improve the function of a future retinal implant, by wholly or partially eliminating the need to perform painstaking calibration of a large-scale retinal interface. This is significant because electrical stimulus calibration generally requires sequential testing of each individual electrode, which is time-consuming, and more importantly is greatly complicated by electrical artifacts, which make it difficult to identify spikes evoked by electrical stimulation. These requirements could pose substantial technical challenges in the clinical setting and highlight the value of the inference approach described. A hybrid approach could entail performing full electrical calibration only for cells with small recorded spikes, to supplement the more reliable inferred electrical activation curves of cells with large recorded spikes. Alternatively, initial estimates of electrical thresholds obtained from recording alone, or using inferred sensitivity as priors in Bayesian estimation (Shah et al., 2019), could be used to more quickly identify thresholds using limited electrical stimulation and recording. Inference also allows for rapid recalibration of electrical stimulation over time, potentially compensating for small movements of an implanted device, which could prove critical for long-term stability.
Footnotes
This work was supported by National Institutes of Health National Eye Institute F30-EY030776-03 to S.S.M.; Polish National Science Center Grant DEC-2013/10/M/NZ4/00268 to P.H.; Pew Charitable Trust Scholarship in the Biomedical Sciences to A.S.; a donation from John Chen to A.M.L.; and Stanford Medicine Discovery Innovation Award, Research to Prevent Blindness Stein Innovation Award, Wu Tsai Neurosciences Institute Big Ideas, and National Institutes of Health National Eye Institute R01-EY021271, R01-EY029247, and P30-EY019005 to E.J.C. Human eyes were provided by Donor Network West. We thank Donor Network West for cooperation; all of the organ and tissue donors and their families for giving the gift of life and the gift of knowledge, by their generous donations; J. Carmena, T. Moore, W. Newsome, M. Taffe, S. Morairty, J. Horton, and the California National Primate Research Center for providing macaque retinas; K. Jensen, M. Allen, and K. Tinajero, K. Berry, K. Williams, B. Morsey, J. Frohlich, and M. Kitano for helping obtain access to macaque retinas and metadata; D. Palanker, F. Rieke, S. Mitra, P. Li, T. Carnevale, M. Hines, F. Rattay, P. Tandon, E. Wu, R. Vilkhu, V. Fan, M. Zaidi, C. Rhoades, N. Brackbill, G. Goetz, S. Wienbar, and the Stanford Artificial Retina team for helpful discussions; K. Kish and J. Weiland for helpful discussions and graciously sharing the base code for the biophysical modeling simulations used in this work; H. Nguyen for contributions to electrical spike sorting in the central retina; R. Samarakoon and S. Kachiguine for technical assistance; and L. Jepson, P. Li, M. Greshner, G. Field, J. Gautier, A. Heitman, and G. Goetz for participating in data collection.
The authors declare no competing financial interests.
- Correspondence should be addressed to Sasidhar S. Madugula at sasidhar{at}stanford.edu