Insights from Computational Modelling: Selective Stimulation of Retinal Ganglion Cells

Improvements to the efficacy of retinal neuroprostheses can be achieved by developing more sophisticated neural stimulation strategies to enable selective or differential activation of specific retinal ganglion cells (RGCs). Recent retinal studies have demonstrated the ability to differentially recruit ON and OFF RGCs – the two major information pathways of the retina – using high-frequency electrical stimulation (HFS). However, there remain many unknowns, since this is a relatively unexplored field. For example, can we achieve ON/OFF selectivity over a wide range of stimulus frequencies and amplitudes? Furthermore, existing demonstrations of HFS efficacy in retinal prostheses have been based on epiretinal placement of electrodes. Other clinically popular techniques include subretinal or suprachoroidal placement, where electrodes are located at the photoreceptor layer or in the suprachoroidal space, respectively, and these locations are quite distant from the RGC layer. Would HFS-based differential activation work from these locations? In this chapter, we conducted in silico investigations to explore the generalizability of HFS to differentially active ON and OFF RGCs. Computational models are particularly well suited for these investigations. The electric field can be accurately described by mathematical formulations, and simulated neurons can be “probed” at resolutions well beyond those achievable by today’s state-of-the-art experimental techniques.

One possible primary reason for the inability of existing retinal neuroprostheses to provide better visual perception may be the indiscriminate activation of different neuronal types across large regions of the retina, providing conflicting information to higher visual centers. To address this problem, we require an improved understanding of how different functional retinal ganglion cell (RGC) types respond to artificial electrical stimulation. In particular, if different RGC types can be selectively or differentially activated in a desired temporospatial sequence, the elicited signals may be interpreted more accurately by the brain, giving rise to visual percepts of greater meaning and utility.
Previous reports have indicated that 1-6 kHz high-frequency electrical stimulation (HFS) may elicit differential excitation of different RGCs in a manner similar to RGC responses to light stimuli in a healthy retina [5][6][7][8]. These studies suggest a great promise in eliciting RGC responses that parallel RGC encoding: one RGC type exhibited an increase in spiking activity during electrical stimulation, while another exhibited decreased spiking activity, given the same stimulation parameters. To test whether a larger range of HFS parameters can improve or even maximize the differential excitation of ON and OFF RGC pathways, we began with in silico investigations using biophysically and morphologically detailed computational models of ON and OFF RGCs using the NEURON computational environment, to evaluate the performance of a range of electrical stimulation amplitudes (10-70 μA) and frequencies (1-10 kHz) on RGC responses.
In addition, in order to investigate the effect of ON and OFF RGC dendritic morphologies on HFS-induced responses, we developed a neural morphology generator, capable of generating RGCs with tunable morphological properties, including the dendritic field radius, total dendritic length, and stratification level. Neuronal morphology has been reported to play a vital role in shaping response properties as well as the integration of neuronal inputs in many cell types throughout the central nervous system (CNS) [9,10]. It is therefore likely that similar morphological dependence is also present in RGCs.
Finally, we used a population-based computational model of ON and OFF layers to explore the performance of electrical stimulation with clinically relevant electrode sizes and locations, as well as stimulation duration.

Computational Model of ON and OFF RGC Clusters
ON and OFF RGC clusters were implemented using the computational software NEURON [11]. Techniques used to model individual RGCs have been described in detail previously [12][13][14].
Firstly, the morphological structures of different RGCs were simulated by a customized neural morphology generator [13,15] (see Fig. 1A). The RGC soma was initially defined as a point at the origin. With the soma as the center, a number of random carrier points, which serve as the basis of dendritic growth, were distributed within a circular planar region having a user-defined radius. An algorithm, based on the minimum spanning tree algorithm [16], generated dendritic branches by connecting unconnected carrier points to nodal points of the tree. At each step, a sweep through all nodes starting from the soma was undertaken to find the carrier point closest to the tree. A cost function was used to calculate the weighted distance e d between a carrier point and a node in the tree, as follows: where d e is the Euclidean distance between a carrier point and a node in the dendritic tree; d p is the length of the path along the corresponding branch from the soma to the carrier point, which is the sum of d e and the length of the branch from the soma to the corresponding node; and b f is a balance factor, which weighs d e and d p against each other in the cost function. The carrier point with the shortest e d was chosen as the candidate point to be connected to the corresponding node. After creating the dendritic tree, the soma was then extended into a 15-μm segment. A 50-μm-long axonal hillock and a 1000-μm-long axon were added subsequently. The vertical distance between the axon and the soma was set to 10 μm, and the first 50-μm segment of the axon was defined as the axon initial segment (AIS). RGC dendritic morphological parameters [17], including dendritic field radius (μm), total dendritic length (μm), and stratification level (μm), were adjusted based on published data for ON and OFF RGC morphologies in guinea pig retina (see Table 1).
Secondly, the ionic model used in this study can be represented by the equivalent cable equation: where V m represents membrane potential, x is the axial cable distance, σ is the intracellular conductivity (mS•cm À1 ), A is the local cell surface to volume ratio (cm À1 ), and membrane capacitance (C m ) per unit membrane area was set to 1 μF•cm À2 . The intracellular axial resistivity (1/σ) was set to 110 Ω•cm. The simulation temperature was 37 C. J ion (mA•cm À2 ) represents the total cell membrane ionic current, consisting of seven time-dependent currents and one leakage current: I Na , I K , I KA , I Ca , and I KCa were defined in a previous RGC model [18], while the hyperpolarization-activated non-selective cationic current (I h ) and the low-threshold voltage-activated calcium current (I CaT ) are newer additions, with both latter currents known to contribute to RGC excitation [19][20][21]. All gating variables, except those for I CaT , satisfied the following first-order ordinary differential equation (ODE): where x is a gating variable, with α x , β x its opening and closing rates, respectively, which are typically functions of membrane potential V m. For I CaT , second-order dynamics were used: where the inactivation process for I CaT was modelled using two transition steps: h T and d T [22]. In this chapter, ON and OFF RGC models shared the same kinetic-defining parameters for all ionic currents. Ionic channel distributions across different cellular regions were set to be cell-specific to reproduce the stimulus dependency of recently published in vitro whole-cell patch-clamping data [14,23]. The estimated kinetic- R radius of dendritic field area, L total length of dendrites, S vertical distance between the soma and the dendritic tree layer. All parameters were estimated based on published data of ON and OFF RGC morphologies [17] defining parameters of each rate and maximum membrane conductance values (mS/cm 2 ) per region in each cell are listed in Tables 2 and 3.

ON and OFF Layer Simulation
To explore the generalizability of HFS-induced differential activation with more clinically relevant electrode size, location, and stimulation duration, we conducted in silico investigations using population-based neuronal models. ON and OFF RGCs were uniformly distributed on a square grid with 40-μm lateral and 40-μm longitudinal distances between neighboring cells. In total, 21 Â 21 cells were simulated in each layer (see Fig. 1B).

Extracellular Electrical Stimulation and Electrode Settings
To simulate extracellular stimulation, we used a mono-polar circular electrode disk with ground located at infinity and approximate the extracellular domain to be homogeneous. The extracellular potential at each spatial point was adapted based on an analytic formula [24][25][26]: where r and z are the radial and axial distance, respectively, from the center of the disk for z 6 ¼ 0 and R is the radius of the disk. The disk potential V o can be determined by the electric stimulation current I and extracellular resistivity ρ e (500 Ω•cm) [25]: For simulation of single ON and OFF RGC stimulation, we defined a 3D Cartesian (x, y, z) coordinate system, with the soma as the origin, so that the upper surface of the RGC dendritic field was aligned in the x-y plane and the RGC axon was aligned with the y-axis. A hexapolar electrode array (each disk electrode of 15-μm radius, with a center-to-center distance of 60 μm) was positioned at the location (0, À40, À50) μm, where (0, 0, 0) μm was the local 3D coordinates of the soma [23]. Cathodic-first, charge-balanced, symmetric, constant-current biphasic stimuli, each with a pulse width of 50 μs per phase, were used without an interphase interval. The extracellular stimulus amplitude ranged from 10 to 70 μA in 5-μA Table 2 Ionic current formulations and kinetic parameters Membrane current For population-based simulations, ON and OFF layers were stimulated using a subretinally placed large diameter electrode (Φ ¼ 200 μm), positioned 200 μm above the RGC soma layer. The stimulus frequency was set to 10 kHz. The extracellular stimulus amplitude ranged from 10 to 350 μA in 10-μA steps. All pulse trains were reduced to 40 ms in duration. All elicited spikes were observed and counted at the soma. Differential activation was determined from the difference in averaged total spike numbers between one cell cluster and the other.

Differential Activation of Individual ON and OFF RGCs
Using a Large HFS Parameter Space  Low-threshold voltage-activated calcium current, I L . leakage current stimulation frequencies (from 1 to 10 kHz). The elicited spikes were observed and counted at the soma. Our model predicted that the elicited RGC spike counts were highly dependent on stimulating frequency and amplitude. The colors in Fig. 2A1-2 denote the number of evoked spikes from ON or OFF cells for a given stimulation frequency and pulse amplitude. The total spike count of the ON cell (panel A1) reached a plateau as the stimulus current surpassed a certain threshold at frequencies up to 6.5 kHz. However, with frequencies higher than 7 kHz, the total spike number increased initially with stimulus amplitude, followed by a sudden decline with further amplitude increases, creating a non-monotonic surface in the frequency-amplitude topological space. In contrast, the OFF cell (panel A2) exhibited a non-monotonic profile at stimulation frequencies higher than 2 kHz. Nevertheless, both activation maps indicated a decreasing stimulation threshold trend with respect to HFS pulse trains with increasing stimulation frequency, where the threshold was defined as the stimulation amplitude capable of eliciting 10% of the maximal spike number of each non-monotonic spike-stimulus profile.
As shown in Fig. 2B, with increasing stimulation frequency, both ON and OFF RGCs exhibited an increased slope of the rising phase in spikes/μA (the epoch in which spike counts increased with increasing stimulation current) and, concomitantly, an earlier onset of the falling phase (in which the total spike numbers saturated or declined). Interestingly, the stimulus-dependent response of the ON cell became relatively stable only at stimulation frequencies higher than 9 kHz, while the OFF cell response tended to be unchanged already at frequencies higher than 5 kHz, thus indicating the ON/OFF-cell-specific frequency dependency.
The differential activation map shown in Fig. 2C provides an alternative visualization of differential activation of RGC types at each stimulation frequency and amplitude. Each grid point was defined as the difference of total spike number (spikes/200 ms) of ON and OFF cells. Our model suggested that differential activation of the ON RGCs was maximized at high stimulation amplitudes (>45 μA) and frequencies (between 3 and 10 kHz). In contrast, HFS pulse trains across all tested frequencies induced robust differential activation of OFF RGCs with different stimulation amplitudes ranging from 10 to 50 μA. Moreover, in Fig. 1C, the threshold at which differential activation began for both cell types gradually reduced as the stimulus frequency was increased from 1 to 6 kHz. The stimulation current range for preferentially activating ON RGCs increased when the stimulus frequency increased from 2 to 5 kHz, then gradually decreased when the stimulation frequency increased from 7 to 10 kHz. In contrast, the stimulation current range for preferentially activating OFF RGCs was mostly stable across all frequencies.
Based on the information provided by Fig. 2, potential optimal stimulation parameter combinations can be chosen to selectively excite ON and OFF cells. For example, with 7-kHz stimuli of 60 μA, the ON RGC was strongly activated, while simultaneously blocking the OFF RGC spontaneous spikes. In contrast, with 1-kHz stimuli of 30 μA, the OFF RGC was strongly activated, while the ON RGC remained silent. . Juxtaposition of the ON and OFF spike count against stimulating amplitude, at frequencies ranging from 1 to 10 kHz. (c). Differential RGC excitability map, defined as the difference of total spike count between ON and OFF RGCs, indicating stimulation parameters which can preferentially activate one cell type, while minimally activating the other type

Simulating Population-Based RGC Activity Under
Clinically Relevant Conditions Figure 3 shows the ability of the model to predict differential activation in a population of 21 Â 21 pairs of ON and OFF cells using a subretinally placed large diameter electrode (200 μm), positioned at 200 μm above the RGC layers. Figure 3A demonstrates the total spike number (spikes/40 ms) elicited among the whole simulated ON and OFF populations. Our model suggested that 10-kHz HFS pulse trains were still able to induce differential activation of ON and OFF RGC populations with different stimulating amplitude parameters. Differential excitability shown in Fig. 3B was determined by the difference in the total spike numbers between the ON population and OFF population. Our model suggested that the activation of the ON RGC population was maximized at higher stimulation amplitudes (>200 μA). In contrast, excitation of OFF population is maximized at amplitudes ranging from 20 to 80 μA. Excitation maps shown in Fig. 3C provide the differential stimulus dependency of ON and OFF populations during HFS. For example, when stimulated at a small stimulation amplitude (60 μA), a large OFF population (C2) was strongly activated, while only local ON RGCs (C1) located close to the electrode were weakly excited. However, when stimulation amplitude was gradually increased above a certain level (>120 μA), local OFF RGCs below the electrode started to be inhibited while a large population of ON RGCs were still being increasingly excited.

Discussion and Conclusion
The simulation results illustrated in this chapter suggest the possibility of translating recent laboratory advances in differential neural activation to large-scale, clinically relevant conditions by (1) relaxing the constraint requiring stimulation electrodes to be near the RGC cell bodies, which are impossible to locate under clinical conditions, and (2) translating the HFS-based differential activation from epiretinal to subretinal stimulation. Computational RGC models provide the ability to investigate neural modulation by changing key stimulation parameters. One advantage of the ⁄ ä When stimulation amplitude was gradually increased to certain level (>120 μA), local OFF RGCs below the electrode started being inhibited while ON RGCs were still increasingly excited computational approach is that the model-generated response space map can be made arbitrarily large and fine-grained for thorough exploration of stimulus parameters. This is difficult, if not impossible, to achieve through biological experiments due to the invasiveness of intracellular recordings. Moreover, simulation can be used to guide in vitro experimental design. For example, it is worthwhile to investigate population-based RGC responses to HFS predicted by our model using high-density multielectrode arrays [27] or using a calcium imaging technique [28].
Discrimination between ON and OFF RGCs with electrical stimulation is an initial step toward improving artificial vision. Until recently, retinal stimulation has not been able to provide differential activation of ON and OFF RGCs. Such co-activation is highly unnatural, providing conflicting information to higher visual centers, and potentially degrading the efficacy of retinal implants. This chapter, built on previous in vitro [5,8,23,29] and in silico [6,7,26,30] studies, demonstrated that preferential or differential activation of individual and population-based RGC types could be achieved. Here, we further showed that the effect was possible over a wide range of HFS parameters. In particular, the ON RGC could be targeted at relatively higher stimulation amplitudes and frequencies, while the OFF RGC could be targeted with lower stimulation amplitudes across all tested frequencies. The precise mechanism underlying differential RGC activation remains largely unknown. Further modeling and in vitro studies are still required to better understand the factors that shape the response of a retinal neuron to biphasic HFS. In particular, efforts should be devoted to assessing the contribution of intrinsic RGC properties including cell-specific ionic channel distributions in shaping RGC spiking profiles. To the best of our knowledge, no studies have identified the distribution of different ionic channel subtypes between the RGC types. Experimental studies in spinal sensory neurons reported that different sodium channel subtypes may respond differentially to high stimulus frequencies [31]. It is therefore likely that a similar frequency dependence is also present in RGCs. Further experiments based on variable ionic channel identifiers [32][33][34] will help us to better understand the reason for the unique stimulus dependency of each RGC type.
The HFS-based stimulation strategy described here may be useful for closely mimicking the natural encoding of RGC visual patterns. Specifically, the ON ganglion cells showed an increase in spike counts (spikes/200 ms) as the stimulus current was increased, while the OFF RGC responses were inhibited by the increased stimulus. In addition, our results suggest that differential activation of the ON RGC may be maximized within stimulation frequencies of 5-10 kHz, as shown in Fig. 2. However, it should be noted that higher frequencies can degrade stimulation efficacy [29,35]. Therefore, a balance between current amplitude and HFS frequency may be necessary for a practical stimulation strategy.
In summary, the modelling approach can predict where the optimal stimulation parameter space is likely to be without detailed experimental investigations, providing insights into stimulation strategies that may contribute further to the development of retinal prostheses.