Foveal eccentricity can influence activation threshold in subretinal electrical stimulation

Purpose: Retinal prosthetic devices have been trialled in patients and encouraging results have been reported. Various factors including retinal disease state and effective electrode placement during surgery, would affect the outcome of retinal prosthesis implementation. In this study, we explore the possible effect of implant foveal eccentricity on retinal activation threshold using computational modelling. Method: Five healthy and five early stage dry age-related macular degeneration human volunteers were scanned using optical coherence tomography. The 3D retinal structure of each individual was extracted from the images and smoothly meshed. A set of five hexapolar electrode configurations was placed in the sub-foveal, 1 mm and 2 mm eccentricity positions from the centre of the retina within the subretinal space. The activation threshold for individual electrodes at all three sites of stimulation were simulated, using bidomain finite element modelling, and compared for all subjects. Results: We found that sub-foveal pit electrode placement required the lowest activation threshold, followed by the 1 mm eccentricity electrode placement. The 2 mm eccentricity electrode placement required the highest activation threshold. Conclusion: Our simulations indicate that clinically minute changes in foveal eccentricity of the subretinal implant can increase activation threshold and therefore, with a given stimulus current, retinal prosthetic performance drops with foveal eccentricity. Translational Relevance: Retinal prosthetic devices have been introduced as a potential therapy for patients suffering from outer retinal diseases. Human subjects with retinal implants demonstrated a varying range of responses to these devices. In this study, we introduce an anatomically-based patient-specific model of the human retina, using it to explore the effect of foveal eccentricity on activation threshold.


Introduction
Retinal prosthetic devices have been introduced as a potential therapy for patients suffering from outer retinal diseases, such as retinitis pigmentosa (RP) and age-related macular degeneration (AMD). The underlying strategies, while different in many aspects, share the common idea of replacing degenerated photoreceptor cells with an array of electrodes stimulating the intact retinal ganglion cells (RGCs) [1]. It is often believed, although unproven, that human subjects with retinal implants were able to exhibit light perception [2,3], detect the orientation of moving light spots [4,5], recognize objects in a high contrast environment [6,7] and recognize large font alphabets [8,9]. Despite these enormous successes, inconsistent performance in patient populations have also been reported [4,10]. In addition to gross retinal physical properties such as overall retinal thickness and the level of retinal scarring, other factors such as the accuracy of electrode placement during surgery [11], as well as foveal eccentricity of the implant [12] are reported to play a significant role in sub-retinal implant performance. Development of a patient  The detected borders were digitized to provide a 3D data-cloud of each retinal surface (for illustration purpose the distance between the retinal surfaces has been enlarged). (C) The 3D data cloud of the retinal surface with realistic dimensions. The OCT imaging protocol and properties as well as our image processing protocol are well detailed in our previous work [22,23]. specific computational model of the human retina, able to provide estimates of activation threshold in individual electrodes, can assist clinicians to implement more effective electrode placements for individual patients. We envision that these patient specific models could also enhance the performance of existing FDA and CE approved medical devices (such as Argus II).
Previous attempts at providing a reliable computational modelling platform for the human retina have focussed on determining the activation threshold of a single RGC [13,14], or providing insight to the spread and sensitivity of retinal cells in response to a single electrode stimulation [15], or the optic nerve using a volume conductor model [16]. On the other hand, continuum mod]elling of the human retina [17,18] provides a macroscopic approximation of the retinal network, allowing researchers to explore the response of the retina to a stimulus electrode array. In such models however, the ultrastructure of the retina is typically simplified to a homogenized rectangular slab [17,[19][20][21], neglecting anatomical variations among individuals. Hence, these models could not replicate the effect of foveal eccentricity on retinal prosthesis performance. In short, in our approach, the RGC layer is modelled as an electrically active volume where, Ve and Vi denote the extracellular and intracellular potentials, respectively, Cm is the membrane capacitance per unit area and Jion is the ionic current per unit area. The intracellular potential is resistively tied to a resting potential Vr through a conductance gr, representing the intracellular connection between the soma and more distal intracellular cellular regions. In this study, we introduce an anatomically-based patient-specific model of the human retina, using it to explore the effect of foveal eccentricity on activation threshold.
The B-scans of each subject were automatically segmented by detecting the five major retinal borders using the OCT machine software (figure 2(A)). The detected borders were sampled at a spatial resolution of 46 μm, reconstructing each retinal surface with 6903 data points ( figure 2(B)). The data cloud was mapped from image coordinates to global Cartesian coordinates to give a 3D discrete retinal representation (figure 2(C)).
The acquired data cloud was used to generate a 3D finite element computer mesh. To do so, each retinal surface was divided into 16×16 regions in the x-y plane with spacing varying from 184 μm at the fovea to 736 μm at the periphery ( figure 3(A)). The partitioning points were linearly connected together, reconstructing each retinal surface using 256 nodes and 225 rectangular elements ( figure 3(B)). Subsequently, the 2D retinal surface meshes were interconnected to form the 3D volumetric retinal geometry ( figure 4).   The mesh, although representative of the general retinal anatomy, had sharp edges and discontinuous surface tangents in adjacent elements (see dark bands in retinal surfaces in figure 4(B). These discontinuities could result in unreliable modelling results. To achieve a smooth uniform mesh, suitable for finite element analysis, an iterative fitting process was applied to the initial mesh, based on an algorithm developed by Bradley et al [24]. Using this fitting process, discontinuities in the initial mesh were removed and the root square error between the data cloud and the associated retinal surface was minimized to 9 μm, as shown in figure 5. This method was applied to the retinas of all subjects to generate anatomical-based meshes. However, the para-foveal OCT images of Subject 5DS were of very poor quality, preventing the image segmentation algorithm from detecting the retinal borders ( figure 6). Therefore, only the sub-foveal pit region was reconstructed, and the other retinal regions of this subject were eliminated from this study.

Retinal electrical formulations
The electrical potential within the deeper retinal layers and the vitreous domain was modelled using a Poisson equation (equation (1)) [20]. The active response of the RGC layer was modelled using bidomain equations (equations (2)-(3)), which is well-documented in our previous work [18,[27][28][29], The zoomed-in subfigure shows dark areas around the boundaries of each element indicating discontinuity in surface normals. The ∼3 mm legend 4 includes the electrically passive layer, which is essential for our modelling approach. The nodes are specified by green points and the retinal surface data clouds are specified by the colour code of figure 2.  Table 2. Initial values of model variables [18].

Variable
Initial value In these formulations, σ e is the extracellular conductivity tensor, σ i is the RGC intracellular conductivity tensor, V e and V m are the extracellular and membrane potentials respectively, and I m is the membrane current determined from equation (4): Here, C m is the cell membrane capacitance per unit area (0.01 μF.mm −2 17 ), and A m denotes the RGC surface to volume ratio. This parameter is dependent on cell density and retinal thickness and varies with eccentricity from the central human fovea. Figure 7 shows the RGC density in different locations of the healthy human retina. In the central fovea ( figure 7(A)), RGC density steeply increases from 2500 to 200 00 cell.mm −2 , reaching an average maximum cell density of 350 00 cell.mm −2 at 1 mm eccentricity from the central fovea [30]. Thereafter, RGC density smoothly declines reaching 100 00 cell.mm −2 at 2 mm eccentricity from the central fovea. The average GCL+IPL thickness of our population study was 49 μm in the central fovea, 117 μm at around 1 mm, and 117 μm at around 2 mm eccentricity. The RGC axon initial segment was assumed to be an electrically compact cylinder parallel to retinal structure with average diameter of 0.46 μm [31] and 90 μm length [18]. Figure 7(B) shows the calculated surface to volume ratio for each foveal eccentricity based on the above parameters [32].
In equation (4), J ion is the ionic current density, modelled based on Abramian et al [18]:

.
Here, V Na , V K and V L are the equilibrium potentials for each ion, g Na , g K , g K.A , g K.Ca , g L are the membrane conductances of each ionic current, and m, n, h, a and h A are gating variables determined by first-order kinetics (equation (6)): where α and β are the forward and reverse gating rates, and are dependent on membrane potential (

Electrode placements and retinal stimulation modelling
To minimize electric crosstalk between concurrent stimulating electrodes, a hexapolar electrode configuration was employed [33]. To minimize electric crosstalk between concurrent stimulating electrodes, a hexapolar electrode configuration was employed consisting of a central, active, stimulating electrode, surrounded by six electrodes, each set to ground voltage (i.e. 0 V) [25,26]. An array of five hexapolar electrodes was implemented in the subretinal space at the centre of the sub-foveal pit, as well as at 1 mm and 2 mm eccentricity positions (figure 8), stimulating the retina using a biphasic cathodic-anodic stimulus, which is explained elsewhere [34], with pulse duration 0.5 ms. The activation threshold for each individual hexagon was measured in all subjects for all three electrode placement sites, i.e. sub-foveal pit, 1 mm and 2 mm eccentricities. For this purpose, the stimulating current was increased from 9 mA.mm −2 in intervals of 1.5 mA.mm −2 . The minimum required current that could evoke an area of activation (V m >60 mV) at the surface of the active region (RGC) within a time window of 0.5 ms from stimulus onset was considered as the activation threshold. The pulse duration (0.5 ms) was chosen to only elicit direct RGC activation [35], and the constraint of evoking an area of activation was chosen to replicate the shape of light perception in experimental studies [2]. In order to incorporate the shunting effect of neighbouring return electrodes on the stimulus current, the threshold of activation for all hexagons was measured simultaneously (i.e. all hexagons were active).

Results
Each hexagon was labelled according to its location within the array from 1 to 5, as well as its placement site. Figures 9 and 10 illustrate the activation threshold of individual electrodes along the entire retina for healthy and diseased subjects respectively. The electrodes were sorted based on their threshold of activation.
Comparing activation thresholds among all cases showed that the sub-foveal pit placement required the lowest activation threshold (mean 12.75± 3.63 mA.mm −2 ), followed by the 1 mm eccentricity electrode placement (mean 14.76±3.40 mA.mm −2 ). The 2 mm eccentricity electrode placement required the highest activation threshold (mean 32.80± 6.87 mA.mm −2 ). There were some exceptions however in this trend for the healthy and, to a lesser extent, for the diseased subjects. Subjects 3 H, 4 H, 2DS and 3DS required a lower activation threshold for the 1 mm eccentricity electrode placement compared to the sub-foveal pit placement. In these cases, the mean electrode distance to target cells was higher at the subfoveal placement compared to the 1 mm eccentricity (table 4), which could account for the decreasing activation threshold at 1 mm in those cases.   Figure 11 illustrates activation thresholds versus electrode distance to target cells for the three electrode placement sites. There was a positive correlation between activation threshold and electrode distance to target in all three sites, albeit with different slopes. The sub-foveal pit electrode placement in the diseased cases showed the lowest correlation between activation threshold and electrode distance to target (R 2 =0.013, p=0.57). These values show that although the linear regression fit doesn't account for much of the variation of the data (R 2 <0.1), it is nonetheless significant (p>0.05). In our subjects, this electrode placement site was the only region experiencing retinal deformation after dry AMD onset. For the sub-foveal pit to 2 mm eccentricity positions, there was increased correlation between activation threshold and electrode distance to target cells, and a better fit was achieved using linear regression (R 2 >0.1 and p<0.05). Moving from sub-foveal pit to 1 mm and 2 mm eccentricity electrode positions, the covariance between activation threshold and electrode distance increased from 0.42 to 0.74 and 0.78 in healthy retinae, and from 0.11 to 0.58 and 0.72 in diseased retinae, respectively. Since the heterogeneity of retinal structure decreased from the sub-foveal pit to the 2 mm eccentricity position, we posit that the stimulus current could traverse a more direct route from electrodes to target cells, and thus the extracellular potential gradient at the surface of the RGC layer would be more dependent on the electrode distance.

Discussion
In our modelling, the sub-foveal pit and 1 mm eccentricity electrode placements required lower activation thresholds compared to the 2 mm eccentricity electrode site, predicting lower retinal prosthesis performance for the larger eccentricity within a particular range of stimulus current. This result agreed with a previous clinical study in humans using subretinal implantation [12], indicating that foveal implantation provides higher performance compared to non-foveal implantation, attributed to the role of the central retina in high-accuracy sharp vision. In our simulations, shorter electrode distance to target cells and larger surface to volume ratios (in the sub-foveal pit and 1 mm eccentricity regions), decreased the activation thresholds, thus predicting higher performance in these regions.
Measuring the activation threshold of each hexagon for all three electrode sites revealed a positive correlation between electrode distance to target cells and activation threshold. This result was in agreement with a clinical study using epiretinal implantation [10]. In contrast, our simulations were undertaken with subretinal stimulation, with outer retinal layer thicknesses chiefly accounting for the electrode distance and therefore, there was a positive correlation between activation threshold and retinal thickness. Furthermore, there were sites in the active RGC region with similar distance to the injecting electrodes, but requiring different levels of current for activation (see figure 11). This has been repeatedly observed in experimental studies, and may be attributed to either the different RGC types present [36], the result of direct or indirect RGC activation [35], or due to the pulse polarity [37]. In our study, this difference in activation thresholds could be explained by the heterogeneity in retinal structure, which ultimately affects the electric current path from injecting electrodes to target cells.
Unlike the experimental study of Rizzo et al [38], our simulations indicated a mean lower activation threshold for diseased retinae (18.6±9.56 mA.mm −2 ) compared to healthy cases (20.92±10.66 mA.mm −2 ). This difference, although not significant (p>0.05), can be explained in several ways. The mean age of diseased retinae in our study was older (67 years) than healthy retinae (25 years), which has potentially lowered retinal thickness for the diseased retinae [39]. Also the onset of disease might have resulted in decreased retinal thickness Figure 9. Threshold of activation for individual hexagons at three placement sites for healthy subjects. The colour bars follow the label of the associated hexagons (red sub-foveal, blue 1 mm and yellow 2 mm). [40]. Therefore, in modelling subretinal direct RGC activation, diseased retinae with shorter electrode distance to target cells (mean 211±27.09 μm) compared to healthy retinae (mean 220±20.03 μm) required lower activation thresholds. Another contributing factor could be the different types of retinal pathology examined. In our study, due to limited incidence of RP, we only had access to AMD subjects, whilst Rizzo et al [38] used RP patients. AMD and RP target different retinal regions (i.e. the central retina in AMD and the peripheral/para-foveal regions in RP), and alter retinal circuitry and structure in different ways. Furthermore, our model parameters and formulations, including RGC activation and retinal layer conductivities, were based on animal studies. These Figure 10. Threshold of activation for individual hexagons at three placement sites for diseased subjects. Data for subject 5DS only included the sub-foveal electrode placement due to low OCT image quality at 1 and 2 mm eccentricities. The colour bars follow the label of the associated hexagon (red sub-foveal, blue 1 mm and yellow 2 mm). parameters can vary between human and animal retinae, which could result in lower activation thresholds in our simulations. Experimental studies using subretinal implants in wild-type and diseased mice showed noticeably higher electrical activation thresholds in diseased retinae compared to healthy cases [41]. In these experiments however, different types of RGC as well as both direct and indirect RGC activation were considered. One of the limitations of our simulations is that only one RGC type and only direct activation were simulated. Modeling of the RGC stimulation based only on retinal thickness is simplistic, since it does not take into account the changes in RPE conductivity, which can become highly heterogenous, especially in advanced AMD. Here, due to lack of available clinical data, retinal layer conductivities were kept fixed in diseased and healthy cases. Another limitation of our study has been the use of AMD OCT imaging data as our 'diseased' dataset. We acknowledge that the most relevant patient group for visual prostheses are inherited retinal degenerations. We also acknowledge that our control and test groups were not age-matched. These limitations were due to clinical data availability at the time of this study. We hope to use this foundation in the future combined with more physiologically appropriate model of retinal pathology such as retinitis pigmentosa. We expect that incorporating changes in retinal layer conductivity as well as including indirect RGC activation will alter the activation threshold for diseased retinae in our modelling.

Conclusion
The aim of this study was to develop patient-specific models of human retinal geometry to investigate the possible effect of foveal eccentricity on retinal implant performance. This was achieved by creating 3D finite element meshes of human retinae from both healthy and diseased subjects, using OCT images and employing Poisson and bidomain equations to represent retinal activation. The key contribution of this approach was to enhance our understanding of 3D current distribution through retinal layers under multi-electrode subretinal stimulation. Our simulations enabled, for the first time, investigation of the effect of retinal curvature on the electrode distance to target cells, and allowed replication of activation threshold variations within and across subjects. Our model could also reproduce the effect of electrode placement (i.e. sub-foveal pit, 1 mm and 2 mm eccentricity positions), on device performance for each subject. These predictions make the model a guide for future surgical retinal prosthesis implantation, allowing for improved individual performance of retinal prostheses.

Author contribution
Farzaneh Shalbaf was the doctoral student on this project [42]. She created the computational model and solved it as part of her doctoral studies. She also created the first draft of this manuscript. Nigel H Lovell was the co-supervisor of this project. He advised on the direction of the manuscript and edited its later versions.
Socrates Dokos was the co-supervisor of this project, and this project is based on his earlier work. He edited the later version of this manuscript.
Mark Trew was the co-supervisor of this project, and this project and oversaw the development of the model, while ensuring its accuracy. He also edited the earlier versions of the manuscript.
Ehsan Vaghefi was the main supervisor of this project. He ensured the efficacy of the computational model created here. He finalized the manuscript and submitted it.

Additional information
There are no competing financial interests to be reported by any of the authors of this manuscript.