Modeling extracellular stimulation of retinal ganglion cells: theoretical and practical aspects

Objective. Retinal prostheses use electric current to activate inner retinal neurons, providing artificial vision for blind people. Epiretinal stimulation primarily targets retinal ganglion cells (RGCs), which can be modeled with cable equations. Computational models provide a tool to investigate the mechanisms of retinal activation, and improve stimulation paradigms. However, documentation of RGC model structure and parameters is limited, and model implementation can influence model predictions. Approach. We created a functional guide for building a mammalian RGC multi-compartment cable model and applying extracellular stimuli. Next, we investigated how the neuron’s three-dimensional shape will influence model predictions. Finally, we tested several strategies to maximize computational efficiency. Main results. We conducted sensitivity analyses to examine how dendrite representation, axon trajectory, and axon diameter influence membrane dynamics and corresponding activation thresholds. We optimized the spatial and temporal discretization of our multi-compartment cable model. We also implemented several simplified threshold prediction theories based on activating function, but these did not match the prediction accuracy achieved by the cable equations. Significance. Through this work, we provide practical guidance for modeling the extracellular stimulation of RGCs to produce reliable and meaningful predictions. Robust computational models lay the groundwork for improving the performance of retinal prostheses.


Introduction
Retinal prostheses aim to provide artificial vision for blind people, using an implanted neurostimulation device [1]. They have been used to treat profound vision loss caused by retinal degenerative diseases, such as retinitis pigmentosa and age-related macular degeneration [1]. The devices are implanted in or near the eye to target inner retinal neurons (ganglion and bipolar cells) that survive, even during late-stage disease [1]. Retinal prostheses have been shown to improve users orientation and mobility, and allow them to locate high-contrast objects [1]. However, the visual acuity possible with artificial vision remains limited, with a best-reported value of 20/400 [2]. To improve outcomes, it is critical to understand the fundamental mechanisms underlying the electrical excitation of retinal neurons. Computational models provide a tool to investigate how various anatomical and biophysical factors contribute to retinal activation, and to design improved stimulation paradigms and devices. In this work, we used a two-part technique to model the electrical stimulation of retinal tissue. We used finite element analysis to model the current flow from the stimulating electrodes. We then coupled the calculated electric fields with multi-compartment cable models, to capture the membrane dynamics of target neurons.
The most widely used family of models describing retinal ganglion cell (RGC) biophysics was initially developed by Fohlmeister, Miller and colleagues at the University of Minnesota [3][4][5][6]. This seminal model was based on whole cell patch clamp recordings from neurons in the tiger salamander retina [3]. They derived cable equations describing ion channel gating kinetics by fitting a model to experimental data, using the same mathematical structure as Hodgkin and Huxley [7]. In 2010, Fohlmeister et al adapted the original model to describe mammalian RGC membrane dynamics using patch clamp data from rat and cat retinal neurons [6]. Importantly, they included four ion channels in the RGC membrane (sodium, calcium, calcium-activated potassium, and delayed rectifier potassium). They used temperature as an independent variable to determine the ion channel distribution across various cellular regions and calculate experimental Q10 values.
The Fohlmeister cable equations have been the basis for hundreds of modeling projects. Despite their widespread use, documentation of RGC model structure and parameters is sparse and inconsistent. Model implementation can influence performance and simplifying assumptions may limit the significance of the model predictions. In this paper, we provide detailed instructions for creating a mammalian RGC cable model, and share all code necessary to replicate the model. Next, we investigated several practical aspects of modeling the extracellular stimulation of RGCs. A critical question is the extent to which the neuron's three-dimensional (3D) shape will influence model predictions. We conducted sensitivity analyses to study how dendrite representation, axon trajectory, and axon diameter influence membrane dynamics and corresponding activation thresholds. Another practical consideration is how to maximize computational efficiency. When modeling a large number of neurons (e.g. a densely populated retina [8]), runtime remains a limiting factor in spite of modern computing resources. First, we optimized the spatial and temporal discretization of our cable model. Second, we applied several simplified threshold prediction theories to compare their prediction accuracy with our RGC cable model. The overall aim of this work was to provide shareable tools for modeling extracellular stimulation of RGCs, to investigate how morphometric factors influence model predictions, and to test several strategies for decreasing computation time.

RGC membrane dynamics
We defined the biophysical properties of our RGC cable model following previous work [6]. We set the temperature of the simulation environment to 37.1 • C. The cytoplasmic (axial) resistivity was 136.6 Ω cm and the membrane capacitance was 1 µF cm −2 across all compartments [6]. The overall differential equation governing membrane voltage was: Table 1. Rate constants for voltage-gated ion channels [6].
For the voltage-gated ion channels (Na + , K + , Ca 2+ ), the channel state variables (m, h, n, c) are governed by equations of the form: The voltage-dependent rate equations for each state variable are shown in table 1. From Fohlmeister et al (2010), the equations were adjusted by the Q10 values corresponding to 37.1 • C (see [6], table 2).
Unlike the other ion channels, the calciumactivated potassium channel (K Ca + ) is ligand-gated according to the following equation: Calcium ion concentration is driven by a pump mechanism, with the equation, where F is Faraday's constant (96 489 • C) and r is the depth (0.1 µm) at which the calcium ion concentration [Ca 2+ ] i is measured: Each ion has a unique Nernst potential, defined in table 2. These potentials are constant, except for calcium.
The maximum ion channel conductance values, listed in table 3, represent the ion channel density in each region. The RGC axon is subdivided into four distinct regions [6,9]. The sodium channel band (SOCB) has an ion channel density 3-5× higher than the neighboring regions, making it prone to excitation [9,10]. Since Fohlmeister et al defined a range of appropriate conductance values, we adopted exact parameters following Raghuram et al [11].
Other physical characteristics of each region (e.g. length, diameter) are provided in table 4. We implemented our multi-compartment cable model and governing equations using NEURON v 7.7 in a Python simulation environment [13]. All code used to build the RGC cable model is online and freely available at: https://github.com/Kathleen-Kish/ Retinal_Ganglion_Cell.

Applying extracellular voltage
To predict the RGC response to extracellular stimulation, we used the finite element method (FEM) to determine the extracellular potential produced by the stimulating electrode at each compartment of the cable model (V e ) [14]. The geometry of our simplified 3D model is shown in figure 1. We conducted finite element analysis in COMSOL Multiphysics Version 5.6 (Stockholm, Sweden) using the AC/DC electric currents (ec) module. This physics interface is used to compute the electric field, current, and potential in conductive media. For these simulations, we represented the active electrode as a 1 A current terminal (surface boundary condition), and assigned bulk tissue conductivities to each domain. We designated the bottom boundary of the sclera as an electrical ground (0 V). We obtained the tissue conductivities and layer thicknesses from prior work [15][16][17]. Using COMSOL, we generated a physics-controlled finite element mesh with extra fine element size. We used a quasi-static solver to calculate the electric potential (ϕ) distribution throughout the mesh. This solver employed the conjugate gradient method to solve Laplace's equation, shown below: ∇ · (σ∇ϕ) = 0 After solving the potential fields generated by the stimulation, we interpolated the spatially dependent FEM solution to find the extracellular potential at the center of each neuron compartment. These V e values were exported from COMSOL as a .txt file and imported to Python. To couple the cable model with the applied extracellular potentials, we used the 'extracellular' mechanism in NEURON (xraxial = 1 × 10 9 S cm −1 , xg = 1 × 10 10 S cm −2 , and xc = 0 µF cm −2 ). Since biological tissue conductivities are predominantly linear at retinal stimulation frequencies (5-100 Hz [1]), we scaled the extracellular potentials calculated for a 1 A current by the timedependent stimulus pulse parameters, and integrated over time to calculate the membrane voltage response [18]. We chose electrode size and stimulus pulse parameters (biphasic, cathodic-first, 0.45 ms/phase) to match the Argus II device [19]. A version of our COMSOL model including the 3D geometry is posted on GitHub, and a full version including the mesh and solution is available upon request.

Cell morphometry
Naturally, RGCs exhibit substantial variations in both somatodendritic and axonal morphometry [20]. To study the effects of cell morphometry modifications on model behavior, we conducted several sensitivity analyses.

Dendritic arbor
First, we studied the importance of including a dendritic arbor as part of an RGC cable model when modeling extracellular stimulation. An RGC model published by Schiefer and Grill does not include any dendrites [21]. On the other hand, many recent publications include a full branched dendritic morphology traced from ex vivo images [6,[22][23][24][25]. Werginz suggests that using a simplified equivalent cylinder to represent the dendrites produces nearly identical results as a more complex model [26]. To reconcile the conflicting approaches in previous studies, we designed an experiment comparing an RGC model with a full branched dendritic arbor, an RGC model with an equivalent cylinder used to represent the dendrites, and an RGC model with no dendrites (figure 2). To identify a suitable dendrite tracing, we searched the Neuromorpho database with keywords 'human' and 'ganglion' [27]. Two publications were identified, one containing five intrinsically photosensitive RGCs [28] and the other containing forty-seven RGCs from the mid and peripheral retina [12]. We selected a single cell tracing (ID: NMO_110421) with a mid-size dendritic field (168 × 183 µm). This is consistent in size with a parasol cell found 2-4 mm from the fovea [29] and would be reasonably targeted by electrodes implanted in the mid-peripheral region [30]. To create the equivalent cylinder model, we used a short vertical dendritic section (length = 10 µm, diameter = 4 µm) attached to a longer horizontal dendritic section (length = 1620 µm, diameter = 2 µm), following the methods described in prior work [26]. The equivalent cylinder approximates the somatodendritic membrane as a lumped impedance by matching the surface area of the original model (11 560 µm 2 ) [31].
We compared the response of the models to extracellular stimulation from a disc electrode [32]. We placed the electrode at a height of 50 µm above the retina, which is within the range of electroderetina distances for current clinical devices [30,33]. We moved the electrode in a two-dimensional grid (1 × 1 mm) above the soma, with a 50 µm step-size. We calculated the action potential threshold at each electrode location using a bisection algorithm (with convergence of 0.1 µA).

Axon trajectory
Next, we investigated the RGC axon pathway as it ascends from the soma and enters the nerve fiber layer. Anatomical studies show that axon trajectories vary among cells [34]. Prior models have made assumptions about the axon; for example, that it follows a 90 • circular arc [21,35] or ascends linearly [9,36].
For this analysis, we systematically explored the influence of RGC axon trajectory on activation threshold. We calculated the path of an ellipse with one vertex at the soma and the other at the nerve fiber layer [37]. We adjusted the distance between ellipse vertices to create variable curvature and steepness (figure 3). We tested RGCs with multiple soma depths (35, 55, 75 µm) to represent natural variations [38]. We compared the response of the models to extracellular stimulation from a disc electrode. We placed the electrode at a height of 50 µm above the retina, and moved it horizontally along the length of the axon. We calculated action potential threshold at each electrode location using a bisection algorithm (with convergence of 0.1 µA).

Axon diameter
As described in table 4, RGC axons are divided into several regions that vary in both physical dimensions and ion channel densities. Early studies identified that RGC axons have a narrow segment with an average length of 75 µm, surrounded by a larger diameter region on either side [39]. More recently, Fried et al used immunochemical staining to identify an SOCB adjacent to the narrow segment, densely populated with voltage-gated Na + channels, that is on average 40 µm long [10]. Experimental thresholds are lowest when the stimulating electrode is placed over this band [10].
Axon diameter in all four regions has been inconsistent in prior computational models. When measured anatomically, RGC axon diameter has been generally shown to scale with cell size [6,11]. We conducted a sensitivity analysis to investigate the effect of axon diameter on membrane dynamics and action potential threshold for extracellular RGC stimulation. We placed the disc electrode in three locations: directly above the soma, above the SOCB (150 µm offset in the +x direction), and above the distal axon (500 µm offset in the +x direction). We adjusted the axon hillock diameter between 2-4 µm and the narrow region diameter 0.6-1 µm [11]. We tapered the SOCB diameter to connect these two regions smoothly [9]. In all simulations, we recorded membrane voltage at all compartments to classify the action potential propagation behavior. We calculated action potential threshold at each electrode location using a bisection algorithm (with convergence of 0.1 µA).

Model run time
In general, more complex models require higher computation time. To maximize efficiency, it is beneficial to find the simplest model that produces reliable predictions. We investigated two strategies for improving the run time of our RGC cable model. First, we conducted a sensitivity analysis to identify the minimum temporal and spatial resolutions necessary to provide consistent predictions. Second, we examined the ability of simplified techniques, such as the activating function (AF), to predict activation threshold.

Spatial and temporal discretization
Multi-compartment cable models predict a spatiotemporally continuous solution for the membrane voltage by solving at a finite number of points in space (sections) and time. NEURON discretizes time and space to solve the relevant partial differential equations using a backward Euler integration method [40]. As a result, the integration time step (∆t) and the spatial interval between sections (∆x) both contribute to the solution accuracy. In general, the runtime of a cable model is directly proportional to the product of ∆t and ∆x [40].
To establish a 'ground truth' solution, we set ∆x = 1 µm and ∆t = 1 µs and calculated the action potential threshold using a bisection algorithm. From there, we incrementally increased the value for ∆x until activation threshold changed by more than 0.1 µA, documenting the change in computation time. We adjusted ∆x independently for each cell region and tested multiple electrode positions to ensure robust results. Then, we similarly increased the value for ∆t incrementally, documenting the change in activation threshold and computation time.

AF for threshold prediction
The second spatial derivative of extracellular potential along an unmyelinated axon, called the AF (AF = ∂ 2 V e /∂x 2 ), drives a neuron's response to an applied stimulus [41,42]. The maximum value of the AF represents the point of peak depolarization on the cell membrane, thus, the most probable action potential initiation site ( figure 4). The AF has also been used as a predictor for threshold [43][44][45][46]. This approach represents a significant reduction in computational demands because it does not require solution of the RGC membrane voltage.
While the AF approach is appropriate for a uniform axon, its validity remains unclear for neurons with non-uniform ion channel density and physical dimensions, such as RGCs. Werginz et al found that while the AF predicts the RGC membrane's initial response to stimulation, uneven axial current flow during the stimulus pulse can limit its relevance [36]. Esler et al suggested using a weighted AF (AF w ) to linearly approximate the cellular integration of transmembrane currents [46]. In this framework, the AF value at each neural compartment (ƒ n ) is multiplied by a weight value (w n ) representing its influence on the SOCB. The authors found that the sum ( ∑ ƒ n ·w n ) across a certain number (n) of compartments could accurately predict RGC firing behavior [46].
We tested the ability of both the AF and weighted AF (AF w ) to predict threshold for our RGC cable model. We used two methods to calculate the AF. First, we used Python's gradient operator to calculate the second spatial derivative of the applied potential field along the trajectory of the RGC. To produce a smooth function with this technique, we discretized the finite element tetrahedral mesh with a cubic (third-order) shape function in COMSOL. Secondly, we implemented a passive RGC model in NEURON and used the 'i_membrane' variable as a proxy for the AF. This variable measures the net transmembrane current density for a given compartment. The transmembrane current (mA cm −2 ) in response to the applied extracellular potential field after one time step (t = 5 µs) is proportional to the AF, while also incorporating differences in axial resistance (R a ) caused by non-uniform section diameter. Finally, to calculate the weighted AF (AF w ), we followed the methods described by Esler et al [46]. We calculated individual compartment weights (w n ) by finding the relative SOCB depolarization resulting from injecting 1 nA current into each compartment, using a passive model. We tested the sum ( ∑ ƒ n ·w n ), varying n = 1 to n = 50 compartments and ordering compartments from highest to lowest weight [46].
To determine the accuracy of these simplified threshold predictors, we calculated their value for 441 uniquely positioned RGCs (figure 5). We used full cable model solutions to determine the threshold electrode current amplitude and corresponding extracellular potential along the RGC. Then, we defined the AF threshold as the AF peak when applying the threshold V e vector. For a reliable predictor, AF threshold should fit an exponential curve based on electrode-axon distance [43,44]. We quantified prediction accuracy by calculating the coefficient of determination (R 2 ) for the exponential curve. We calculated electrode-RGC distance from the center of the electrode to the center of the soma.

Dendritic arbor representation influences activation threshold map
To characterize the effects of dendritic arbor on activation, we built three versions of an RGC cable model: (1) with a full-branched dendritic morphology traced directly from a human cell, (2) with an equivalent cylinder used to represent the dendrites, and (3) a simplified model with no dendrites (see figure 2). We compared the response of each cell model to extracellular stimulation. Figure 6 shows the distribution of activation thresholds (µA) in response to a biphasic stimulus pulse, when a disc electrode was moved in a 1 × 1 mm grid above each cell model. Compared to the full branched morphology (figure 6(a)), the equivalent cylinder representation (figure 6(b)) predicted the same absolute minimum threshold (21.7 ± 0.1 µA) located above the SOCB. However, the elongated dendritic geometry created a low threshold region spatially overlapping the equivalent cylinder. Removing the dendrites altogether (figure 6(c)) increased the threshold magnitude (absolute minimum: 26.4 ± 0.1 µA), especially when the electrode was near the soma. This model lacks the active ion channels on the dendritic membrane, which can integrate the voltage produced by the electrode and increase the overall excitability of the neuron. Based on these results, we used the fullbranched dendritic morphology for our remaining analyses.

Axon curvature modulates activation threshold profile
We altered our RGC cable model by creating a series of elliptical axon trajectories with variable curvature and steepness, including several soma depths (see figure 3). Again, we compared their response to biphasic extracellular stimulation with a disc electrode. Figure 7 shows the profile of activation thresholds (µA) as we shifted the stimulating electrode horizontally along the length of each axon. The minimum threshold always occurred when the electrode was above the SOCB. In general, steeper  trajectories with sharper curvature had lower activation thresholds, particularly when the electrode was near the soma. The maximum threshold decrease (between gradual and steep trajectories) was 37%, 53%, and 60% for somata at depths of 35, 55, and 75 µm, respectively.
Steeper trajectories may have lower activation thresholds for two reasons: smaller SOCB-electrode distance and greater SOCB angle. Both factors influence the electric field gradient across the SOCB, which drives activation. To separate these effects, we looked at two scenarios. First, we identified axons with an SOCB angle of 6 • , and varying SOCB-electrode distance ( figure 7(d)). Increasing the SOCB-electrode distance systematically increased threshold. Then, we identified axons with an SOCBelectrode distance of 140 µm, and varying SOCB angle ( figure 7(e)). Increasing SOCB angle systematically decreased threshold. This supports the claim that threshold decrease for steep trajectories is due to a combination of SOCB-electrode distance and SOCB angle.

Axon diameter affects action potential propagation
The final morphometry modification to our model was to change axon diameter (specifically in the axon hillock and narrow region) and evaluate action potential propagation. First, we adjusted the axon hillock diameter between 2 and 4 µm. Figure 8 shows the membrane voltage dynamics of our model as we decreased axon hillock diameter. In this analysis, we applied the threshold stimulus pulse with a disc electrode directly above the soma. When the axon hillock diameter was 4 µm, the action potential propagated initiated in the SOCB and propagated in both directions. When the axon diameter was 3 µm, there was a prolonged latency before the action potential invaded the soma. This provided time for the SOCB membrane to repolarize and experience a secondary 'echo spike' , which propagated down the axon. When the axon hillock diameter was 2 µm, the action potential could not depolarize the somatodendritic membrane. The same behavior occurred when the stimulating electrode was located above the SOCB (150 µm offset) or distal axon (500 µm offset). As shown in figure 8, there was a minor decrease in activation threshold with axon hillock diameter when the electrode was above the soma. No change in threshold was observed when the electrode was above the SOCB (threshold: 22.5 ± 0.1 µA) or distal axon (threshold: 38.1 ± 0.1 µA).
Adjusting the narrow region diameter between 0.6 and 1.0 µm caused analogous trends in the model. Action potentials propagated bi-directionally only when the narrow region diameter was greater than 0.8 µm. Narrow region diameter had no influence on activation threshold. In some situations, increasing Na + and K + conductance in the narrow region could counteract a smaller diameter. For example, settinḡ g Na = 250 mS cm −2 andḡ K = 125 mS cm −2 allowed for bidirectional propagation, even with a 0.6 µm narrow region.

Optimized spatial and temporal discretization
At first, we built our multi-compartment cable model with a compartment length of 1 µm and integrated with a time step of 1 µs. We incrementally decreased the spatial and temporal resolution, measuring the change in activation threshold and computation time. Table 5 shows the spatial resolution of our final model by region. The 'ground truth' model contained 4260 sections, while the reduced model had 960 sections. Computation was 7.4 times faster for the spatially optimized model, with no differences in predicted threshold (±0.1 µA). Table 6 summarizes the effects of increasing the integration time step on simulation time required to calculate threshold for one neuron. We chose a value of ∆t = 5 µs for our final model. Overall, performing this sensitivity analysis allowed us to solve the RGC cable equations 12.75 times faster, with a negligible influence on activation threshold.

Simplified threshold predictors fall short of full cable model solutions
We assessed both an AF approach and a weighted AF (AF w ) approach to reduce the computational demands required to predict threshold. We compared these simplified approaches to the threshold predictions generated with the full RGC cable model. First, we calculated the second spatial derivative by applying Python's gradient operator twice to the extracellular potentials along the neuron. We calculated the AF threshold and spatial location of the AF peak for 441 uniquely positioned RGCs. Figure 9(a) shows AF threshold plotted against electrode-RGC distance. The exponential fit had an equation of y = 0.67 × 10 −9.23x + 0.25 with R 2 = 0.04. The low coefficient of determination indicates poor threshold prediction accuracy. Figure 9(a) also indicates the spatial location of the AF peak using color. The AF peak was most often located at the axon hillock, due to its changing orientation with respect to the soma (e.g. soma aligned with z-axis, axon hillock begins shifting in x-direction). However, the peak occurred in the axon when it passed directly beneath the disc electrode. The distribution of peak locations contradicts the cable model, in which the action potential always initiated in the SOCB.
Using transmembrane current (mA cm −2 ) after one time step (t = 5 µs) as a substitute for the AF only produced modest improvements in prediction accuracy. We calculated the threshold and spatial location of the peak transmembrane current for 441 uniquely positioned RGCs. Figure 9(b) shows 'i_membrane' threshold plotted against electrode-RGC distance. The exponential fit had an equation of y = 4.71×10 −14.4x + 0.52 with R 2 = 0.19. By accounting for non-uniform section diameter, the instantaneous transmembrane current incorporates differences in axial resistance (R a ) and provides a better estimate of the AF peak. The peak was more likely to occur in the large diameter soma and less likely to occur in the narrow axon. However, the AF peak was still only located in the SOCB 30% of the time with this method, which does not agree with the cable model.
The AF w approach involved scaling the AF by specified compartment weights and finding the sum across n compartments ( ∑ ƒ n ·w n ), ordering compartments from highest to lowest weight [40]. We first derived compartment weights (w n ) by finding the SOCB depolarization resulting from injecting 1 nA current into each compartment, using a passive model, and normalizing [46]. Figure 9(d) shows the resulting weights, coded by color. We calculated AF w across n = 1-50 compartments for 441 uniquely positioned RGCs at threshold. Figure 9(c) shows weighted AF threshold ( ∑ ƒ n ·w n ) plotted against electrode-RGC distance when n = 25. Notably, this approach generated both positive and negative values that were largely dependent upon the AF polarity at the SOCB. The polarity varied based on the neuron's location with respect to the stimulating electrode (figure 9(e)). As a result, we did not see a consistent value emerge at threshold, no matter how many compartments were included in the sum.

Discussion
This work had three main outcomes. We provided a functional guide for modeling extracellular RGC stimulation. We described how morphometric factors influence model predictions, adding to the prior work. Finally, we determined temporal and spatial resolutions that optimize run time versus accuracy.
RGC dendrites contain active ion channels that contribute to action potential generation [4,48]. Prior models have made simplifications including eliminating dendrites or using an equivalent cylinder representation [21,26]. Our sensitivity analysis showed that including a full branched dendritic morphology was important to produce reliable activation threshold maps. We limited our analysis to a single parasol cell in the mid-peripheral region [12,29]. Prior work demonstrated that increasing dendritic field diameter from 150 µm to 500 µm had no significant influence on threshold [47]. Therefore, while including a full branched dendritic morphology is important for future models, the dendritic field diameter of parasol cells is unlikely to influence predictions. On the other hand, developing models of midget cells, which are prevalant in the foveal region and have much smaller dendritic arbors, should be investigated [49]. Unfortunately, no human midget cell tracings were available on the Neuromorpho database at the time of this publication.
RGC axons ascend from a soma in the inner retina to the nerve fiber layer, and the path they take varies naturally among cells [34]. Prior models have made various assumptions about the curvature of this path [9,21,35,36]. Our sensitivity analysis revealed that RGC axon trajectory influences activation thresholds, specifically due to the orientation and distance of the SOCB in relation to the stimulating electrode. Regardless, the lowest thresholds always occurred when the electrode was above the SOCB, consistent with experimental results [10]. Furthermore, the average threshold profile is similar across all soma depths. These results demonstrate that future models should incorporate a range of axon trajectories or use the average trajectory and clearly state their assumptions.
The diameter of RGC axons (including their various subregions) has not been consistent in prior models. In our cable model, axon diameter had a minimal effect on activation threshold, but did influence action potential propagation. Prior models of intracellular current injection have similar findings. A model of neocortical pyramidal neurons established that axon hillock diameter can influence the efficacy with which a spike invades the soma [50]. Sheasby and Fohlmeister found that the diameter of the narrow region can influence whether an action potential will propagate uni-directionally or bi-directionally [5]. Importantly, cells with a high somatodendritic surface area and low axonal surface area can experience increased latency for an axonal spike to enter the soma or even failure of a spike to enter the soma altogether [5]. Our sensitivity analysis showed that for action potentials to propagate bi-directionally with extracellular stimulation, we must similarly avoid a large impedance (i.e. surface area) mismatch between the axon and soma. Based on experimental evidence, we believe that action potential propagation in both the orthodromic and antidromic direction is most realistic. Specifically, calcium imaging data shows somatic activation in response to axonally initiated spikes [51]. Therefore, the surface area of the axon (which is directly proportional to diameter) must be large enough to depolarize the somatodendritic surface area, given certain ion channel densities. For our particular cell tracing, an axon hillock diameter of at least 3 µm and narrow region diameter of at least 0.8 µm were necessary. Additionally, there was a small range of axon diameters that caused an axonal 'echo spike' (see figure 8(b)). It is difficult to garner experimental evidence about whether echo spikes occur in nature or are simply an artifact of the cable model. Therefore, it may be up to a modeler's discretion whether to permit them in a simulation.
We also investigated two strategies for improving the computational efficiency of our RGC cable model. First, we optimized the temporal and spatial resolutions. The optimized section lengths are summarized in table 5. We found an ideal integration timestep of 5 µs. Computation time is proportional to the product of ∆t and ∆x, allowing our optimized cable model to run 12.75 times faster without impacting threshold predictions. Given that the original full-resolution cable model took 112 min to run, this reduction can be significant depending on the number of neurons included in the simulation.
Secondly, we found that simplified threshold prediction techniques could not reliably replicate the predictions generated with our full RGC cable model, in agreement with Werginz et al who observed that spikes are generally generated in the SOCB [36]. The varying ion channel densities and diameters between cell regions violate the assumptions of the AF framework that require homogeneous cell properties [41]. Figure 10 exemplifies how the AF, which is proportional to the instantaneous transmembrane current, does not clearly predict cell firing behavior. In figure 10(a), we see an example where the SOCB is hyperpolarized at the pulse onset, but the transmembrane current changes direction mid-pulse. By the end of the cathodic phase, current is peaking in the SOCB, instigating spike initiation. This behavior is caused by the high SOCB conductance (ion channel density) compared to neighboring regions. Figure 10(b) shows a comparison where the AF peak is located in the SOCB to begin with, and subsequent depolarization and spike initiation occur more rapidly. In general, the AF predicts the initial membrane response, and largely depends on where the cell is located in relation to the stimulating electrode (see figure 9(e)). However, for RGCs, this relationship could not reliably indicate what would happen by the end of the pulse, due to disproportionate axial current flow into the SOCB. We could not overcome these unpredictibilities, even by using a weighted AF as suggested by Esler et al [46]. Perhaps this disparity is because they drew conclusions based on a limited Figure 10. Example of RGC firing behavior that cannot easily be explained by the AF. (a) The left plot shows membrane voltage over time, with the red arrow indicating the initial hyperpolarization of the SOCB. Despite this hyperpolarization, the action potential still initiates in the SOCB. The middle plot shows the transmembrane current for each compartment at the time of the pulse onset (proportional to AF). The SOCB compartments are highlighted red, and SOCB current is negative. The right plot shows spatial transmembrane current at the end of the pulse. The SOCB compartments are highlighted red, and SOCB current has changed polarity. (b) Example where the SOCB is depolarized right away. The first plot shows membrane voltage over time. The second plot shows spatial transmembrane current at the pulse onset (proportional to AF), with the SOCB compartments highlighted in red. Immediately, SOCB current is positive. Overall, these plots highlight the challenges of using simplified threshold predictors derived from second spatial derivative of the extracellular potential field for non-uniform axons. (c) Location of the neurons in (a) and (b) with respect to the stimulating electrode. number of RGC locations, reducing the variability of AF shape.
There were several limitations of this work. First, we did not analyze the effects of changing regional maximum ion channel conductances. In all simulations, we used the constant values provided in table 3. These values are consistent with experimentally derived ion channel densities across four mammalian RGCs [6]. However, future work could include a more systematic analysis of ion channel conductance, spanning across the physiological range. Additionally, we limited our model to epiretinal stimulation, in which the stimulating electrode is located directly above the nerve fiber layer. In doing so, we considered only direct RGC activation, disregarding indirect activation via other cells (bipolar, amacrine) in the retinal network. For subretinal electrodes, it may be important to include a network model. Finally, RGC structural changes (e.g. dendritic field reduction, neurite sprouting) that may be induced by retinal degeneration were not included in our model [52,53]. The presence and severity of these structural changes likely depends on the stage of disease progression and retinotopic location [52,53].
For the purpose of this work, we calculated the extracellular potentials generated by the stimulating electrode using finite element analysis in COMSOL. In certain situations, a researcher may not have access to this software or may want to simplify their simulations. It is possible to calculate extracellular potential generated by a disc electrode using a simplified equation, as described by Wiley and Webster [54].
However, using this approach would disregard the effects of non-uniform conductivity for various tissue types. Alternatively, an open-source finite element analysis software (e.g. FEBio, Netgen) could be used.

Conclusion
Through this work, we aim to provide practical guidance for modeling the extracellular stimulation of RGCs to produce reliable and meaningful predictions. Additionally, we intend to increase the accessibility of these methods by sharing our code. Reliable computational models lay the groundwork for improving the performance of retinal prostheses. They allow for the design of novel stimulation paradigms and hardware, which could ultimately improve the quality of life for millions suffering from retinal degenerative diseases.

Data availability statement
All data that support the findings of this study are included within the article (and any supplementary files).