Biophysical model of axonal stimulation in epiretinal visual prostheses

Visual prostheses aim to restore vision to people blinded from degenerative photoreceptor diseases by electrically stimulating surviving neurons in the retina. However, a major challenge with epiretinal prostheses is that they may accidentally activate passing axon fibers, causing severe perceptual distortions. To investigate the effect of axonal stimulation on the retinal response, we developed a computational model of a small population of morphologically and biophysically detailed retinal ganglion cells, and simulated their response to epiretinal electrical stimulation. We found that activation thresholds of ganglion cell somas and axons varied systematically with both stimulus pulse duration and electrode-retina distance. These findings have important implications for the improvement of stimulus encoding methods for epiretinal prostheses.


I. INTRODUCTION
Microelectronic retinal prostheses aim to restore vision to people blinded from retinal degenerative diseases by electrically stimulating surviving neurons in the retina. The evoked neuronal responses are transmitted to the brain and interpreted by the patient as visual percepts ('phosphenes'). Two devices are already approved for commercial use: Argus II (epiretinal, Second Sight Medical Products Inc. [1]) and Alpha-IMS (subretinal, Retina Implant AG [2]). In combination with stem cell therapy and optogenetics, various sight restoration options should be available within a decade [3].
However, a major challenge with current retinal prostheses is the inability to achieve focal tissue activation. Because epiretinal prostheses sit on top of the optic fiber layer, these devices may accidentally stimulate passing axon fibers, which could antidromically activate cell bodies located peripheral to the point of stimulation [4], [5]. This can cause nontrivial perceptual distortions [6], [7] that may severely limit the quality of the generated visual experience [8].
Although previous studies have modeled the effect of epiretinal electrical stimulation on the ganglion cell response, most of them either focused on a single neuron (e.g., [4], [9], [10]) or did not consider axons as potential sites of spike initiation (e.g., [9], [11]). In addition, axonal stimulation has been shown to vary as a function of various stimulus parameters [5], the theoretical foundation of which remains poorly understood.
To examine this issue, we developed a computational model of a small population of morphologically and biophysically detailed retinal ganglion cells, and simulated their response to epiretinal electrical stimulation.

A. Neuron model
The spatial setup of the model is shown in Fig. 1. We populated a ganglion cell layer with a simulated mouse ganglion cell based on 3D morphological data from neuromorpho. org (Cell 030102, ID NMO 06380 [12], [13]). The cell's soma and dendritic tree occupied a 3D box of size 79 µm × 229 µm × 18.7 µm. Starting at the soma, the axon first traveled upward to a distance of z ≈ 30 µm, before it bent and traveled along the surface of the simulated optic fiber layer. Ensuring the axon extended well beyond the test area of the extracellular stimulating electrode, we linearly extended the axon by ≈ 900 µm, starting from the end point of the experimentally traced axon. We then placed a simulated disk electrode (200 µm diameter) at varying distances (100 − 800 µm) above the optic fiber layer.
We considered five nonlinear ion channels: Na + , Ca 2+ , non-inactivating K + (delayed rectifier), inactivating K + (A type), and Ca 2+ activated K + , as well as a leak current [14]. Following Kirchoff's law, the membrane potential of each compartment was given as: where C m = 1 µF cm −2 , E N a = 35 mV, E Ca = 132 mV, E K = −75 mV, and E leak = −65 mV. The rate constants for m, h, c, n, a, and h A all solved the first-order kinetic equation: where both α and β were functions of voltage (see Table I for channel gating kinetics).
The gating of I K,Ca was modeled as: where the Calcium dissociation constant, [Ca 2+ diss ], was set to 1 × 10 −6 M, and j was set to 2 [10]. Internal Calcium concentration, [Ca 2+ ], was allowed to vary in response to I Ca . We assumed that the inward flowing Ca 2+ ions are distributed uniformly throughout the cell and that the free Ca 2+ above a residual level, [Ca 2+ ] res = 1 × 10 −7 M, was actively removed from the cell with a time constant τ Ca = 50 ms: where F was the Faraday constant, 3/r was the ratio of the surface to volume of the spherical cell soma, and the factor of 2 on F was the valency. Following [4], [14], the five ion channels were distributed with varying densities across dendritic, somatic, and axonal compartments, simulated by varying the value of g for each channel (see Table II).
All simulations were run either on an Intel(R) Core(TM) i7-5930K CPU operating at 3.5 GHz using the Brian simulator [15], or an NVIDIA Titan Xp graphics card using the Ion Channel State var α β Na + Sodium m   [4], [14]). brian2cuda extension. Typically, 10 µs time steps were used.
To ensure fine numerical solutions, the ganglion cell was modeled with almost 3000 compartments.

B. Neural population model
The single cell model was replicated to produce a population of 11 × 6 retinal ganglion cells distributed over a 1000 µm × 500 µm area. To simulate a range of excitation profiles, we uniformly sampled values for g Na , g K , and g A from the naturally occuring range of values reported in [4], [14] (see ranges in Table II).

C. External stimulus application
External stimulation was given either by an idealized point source or a circular disk electrode.
The electric field V (r) of a point source was given by: where σ = 110 Ω cm was the resistivity of the retinal extracellular solution (typically Ames medium), I was the amplitude of the constant current pulse, and r was the distance from the stimulating electrode to the point at which the voltage was being computed. Nonuniformities in the electric field arising from the presence of the model cell were not considered. For stimulation by a circular disk electrode in a semiinfinite homogeneous isotropic half-insulating volume conductor, the extracellular potential was given by [16]: for z = 0, where r and z were the radial and axial distances from the center of the disk, V 0 was the disk potential, σ was again the medium conductivity, and a = 100 µm. The extracellular potential was calculated at the center of each   [4], [9]). compartment and uniformly applied to the surface of the entire compartment. The constant voltage model of the disk electrode may be converted to a constant current model (since the extracellular space is modeled as purely resistive) with the addition of a constant multiplicative factor [9], [10].

A. Excitation thresholds along the axon of a ganglion cell
In physiological experiments with retinal ganglion cells and a small stimulating electrode, the lowest activation thresholds occurred when the electrode was positioned over the sodium channel band (located on the proximal axon, 40 − 80 µm from the soma) [4]. To assess the activation threshold of our model ganglion cell, we applied a 0.2 ms monophasic cathodic pulse originating from a point source (see Eq. 5) located at z = 50 µm above the axon. The result is shown in Fig. 2. Similar to rabbit ganglion cells [4], the model cell exhibited its lowest activation threshold when the electrode was centered over the sodium channel band, with thresholds rising across the thin segment and eventually plateauing in the distal axon.

B. Dynamic range to achieve focal stimulation
The effectiveness with which a stimulus activates only neurons near the electrode can be estimated by studying the relative thresholds for activation of the somatic region (including the axon hillock, sodium channel band, and thin segment) versus activation of the (distal) axon. For example, if thresholds are higher in the distal axon than in the soma, it may be possible to avoid axonal stimulation by utilizing a relatively low stimulus amplitude. We thus measured the range of currents that activated ganglion cell somas without additionally activating axons of passage ('dynamic range') by systematically varying the vertical offset (z) of an epiretinal disk electrode from the population of simulated ganglion cells. We gradually increased the current amplitude of a 0.2 ms monophasic cathodic pulse until the simulated cells started spiking, and identified the site of spike initiation (i.e., soma versus axon). The resulting excitation thresholds are shown in Fig. 3. Electrode-retina distances below 200 µm tended to activate somas and axons with similar thresholds, thus suggesting that axonal stimulation cannot be avoided at these distances. However, as the vertical offset of the electrode increased, both thresholds and dynamic range increased. The best dynamic range was achieved at an electrode-retina distance of 800 µm.

C. Response of population model to epiretinal stimulation
A recent electrophysiological study found that stimulus pulse durations two orders of magnitude longer than those typically used in existing epiretinal prostheses can avoid activation of ganglion cell axons [5]. The authors of the study showed that this was due to pulse durations on the order of 10 − 100 ms activating bipolar cells, which in turn activated  Maps are oriented such that the optic disc lies to the right of the image.
ganglion cells, thus confining retinal responses to the site of the electrode. Interestingly, we found a similar effect when we simulated the population response to epiretinal stimulation of 0.1 − 50 ms pulse duration (Fig. 4). Here, we first determined the threshold current for ganglion cells located directly below the electrode. Each dot in Fig. 4 then indicated the number of spikes emitted by a ganglion cell in response to stimulation at 2× threshold. As pulse duration increased, the spatial extent of ganglion cell activation was noticeably reduced.

IV. DISCUSSION
We found that activation thresholds of ganglion cell somas and axons varied systematically with both stimulus pulse duration and electrode-retina distance. A key prediction of the model is that focal stimulation (i.e., to activate ganglion cell somas without also activating their axons) might be achieved at relatively large electrode-retina distances. However, such distances are also known to give rise to higher activation thresholds and larger phosphene sizes [17]. In addition, the use of longer pulse durations led to a reduction in the spatial extent of axonal activation. Previous studies [5] have attributed this effect to increased activation of bipolar cells, which in turn activate ganglion cell somas. However, since the present model did not include bipolar cells, the observed effect can be attributed to the spatiotemporal integration properties of the ganglion cell model.
These results suggest that focal activation of the retina might be achieved by a clever combination of stimulus parameters and electrode configuration. In the future, we will extend the model to feature diverse classes of ganglion and bipolar cells. Furthering our understanding of the retinal response to various stimulation patterns and electrode configurations will remain crucial to the continued improvement in the design of epiretinal visual prostheses.