Patient-specific computational models of retinal prostheses

Retinal prostheses stimulate inner retinal neurons to create visual perception for blind patients. Implanted arrays have many small electrodes, which act as pixels. Not all electrodes induce perception at the same stimulus amplitude, requiring clinicians to manually establish a visual perception threshold for each one. Phosphenes created by single-electrode stimuli can also vary in shape, size, and brightness. Computational models provide a tool to predict inter-electrode variability and automate device programming. In this study, we created statistical and patient-specific field-cable models to investigate inter-electrode variability across seven epiretinal prosthesis users. Our statistical analysis revealed that retinal thickness beneath the electrode correlated with perceptual threshold, with a significant fixed effect across participants. Electrode-retina distance and electrode impedance also correlated with perceptual threshold for some participants, but these effects varied by individual. We developed a novel method to construct patient-specific field-cable models from optical coherence tomography images. Predictions with these models significantly correlated with perceptual threshold for 80% of participants. Additionally, we demonstrated that patient-specific field-cable models could predict retinal activity and phosphene size. These computational models could be beneficial for determining optimal stimulation settings in silico, circumventing the trial-and-error testing of a large parameter space in clinic.


Introduction
Electronic visual prostheses activate neurons in the visual pathway to create light perception for blind patients 1 . Electric current alters the transmembrane potential of nearby neurons, opening voltagesensitive ion channels to induce action potentials. Arti cially induced spots of light are called phosphenes, and visual prostheses aim to create interpretable scenes composed of phosphenes. Retinal prostheses are implanted at the back of the eye and activate inner retinal neurons for patients with degenerated photoreceptors 2 . Clinical testing has demonstrated that people with retinal prostheses can detect large objects, and in some cases distinguish letters, although at a rate much slower than natural reading 3 . Retinal prostheses still function ten years post-surgery, which supports the safety of long-term implantation and stimulation 4 .
Despite advancements in the eld of arti cial vision, retinal prostheses continue to face substantial limitations. Even the best-restored visual acuity is below the threshold of legal blindness 5 . The future success of these devices depends on their ability to activate target neurons with improved spatiotemporal resolution, while avoiding off-target stimulation 6 . To achieve this, researchers are developing microelectrode arrays with an increasing number of electrical contacts. Early commercial retinal prostheses, like the Argus II, have fewer than one hundred electrodes 7 . However, future generation implants could have hundreds or thousands 2 . Not all electrodes will induce a percept at the same current amplitude, requiring the establishment of a visual perception threshold for each individual electrode.
Manually programming these devices in clinic using a trial-and-error process will place an excessive burden on both clinicians and patients. Furthermore, phosphenes created by single-electrode stimuli can vary in shape, size, and brightness 8 .
Computational models provide a tool to predict inter-electrode variability and automate device programming. Prior statistical models have found a signi cant correlation between electrode-retina distance and perceptual threshold [9][10][11] . Other data-driven models have been used to predict the visual perception resulting from a set of retinal stimulation parameters 12,13 . An alternative to the data-driven approach, eld-cable models aim to explicitly model the retinal circuitry and its response to stimulation [14][15][16][17][18][19][20] . This two-part technique uses a three-dimensional bioelectric eld model to solve for the spatial distribution of electric potential throughout the tissue, before calculating the effects on neural activity with cellular cable models 21 .
In this study, we investigated inter-electrode variability across seven epiretinal prosthesis users with both statistical and patient-speci c eld-cable models. We evaluated the utility of these computational approaches for predicting perceptual threshold. Additionally, we demonstrated that patient-speci c eldcable models constructed from imaging data are useful for predicting retinal activity and phosphene size.

Statistical models for perceptual threshold
We created linear regression models to determine the in uence of four explanatory variables on visual perception threshold for epiretinal disc electrodes. The explanatory variables were electrode impedance, electrode-retina distance, retinal thickness beneath the electrode, and brotic tissue thickness on the electrode surface (if present). We collected data from seven participants implanted with the Argus II retinal prosthesis (Second Sight Medical Products Inc., Sylmar, CA). Table 1 summarizes the results of our regression analyses, including coe cients of determination and p-values. For most participants (5/7), there was a signi cant negative correlation between retinal thickness and perceptual threshold (i.e., Fig. 1a). For some participants (3/7), there was a signi cant positive correlation between electrode-retina distance and perceptual threshold (i.e., Fig. 1b). Two participants demonstrated a signi cant negative correlation between electrode impedance and perceptual threshold. Counter to our expectation, we found no signi cant correlation between brotic tissue thickness and perceptual threshold.
Based on the results of our multiple regression analysis, the measured explanatory variables accounted for 55-85% of the variability in perceptual threshold for all but two participants. For CL01 and CL02, they only explained 5-11% of the variability in perceptual threshold. It is worth noting the restricted range of electrode-retina distances and perceptual thresholds for these two participants (Fig. 1c). In other words, the epiretinal electrodes were close to the retina and perceptual thresholds were generally low.
We t linear mixed models with random intercept and slope to determine if the in uence of explanatory variables on threshold was non-zero on average, and if the effect varied across participants (Fig. 1d).
Only retinal thickness had a signi cant xed effect on perceptual threshold (p = 0.02). The effects on perceptual threshold of retinal thickness, electrode-retina distance and electrode impedance varied by individual (p < 0.0001). Patient-speci c eld-cable models to predict perceptual threshold From the same group of subjects, we built patient-speci c eld-cable models for ve eyes implanted with the Argus II retinal prosthesis 19 . We excluded two subjects because no signi cant explanatory variables were identi ed during regression analysis. We constructed three-dimensional models from imaging data and used a two-part technique to model the electrical stimulation of retinal tissue. First, we calculated the electric elds generated by each stimulating electrode using nite element analysis 19,[21][22][23] . Second, we functionalized the bulk tissue models with multi-compartment cable models of retinal ganglion cells (RGCs) to predict neural activity in response to the electric elds 19 . Using these models, we calculated the action potential threshold for RGCs beneath each stimulating electrode.
We compared the action potential thresholds predicted by patient-speci c eld-cable models with participant visual perception thresholds using linear regression (Fig. 2). Our model predictions were signi cantly correlated with perceptual data for 4/5 participants.
In general, thresholds predicted in silico were lower than human reported thresholds. This was especially true when the electrode-retina distance was less than 100 µm. Changes to model conductivity values within the biological range did not eliminate this discrepancy 24,25 . Patient-speci c eld-cable models to predict phosphene shape and size In addition to threshold prediction, the patient-speci c eld-cable models can estimate the activation patterns of retinal tissue (Fig. 3). These patterns depend on the position of the electrode with respect to the retina and the angle of nearby retinal ganglion cell axons 12 .
We compared the patient-speci c eld-cable model predictions of retinal activation to a prior study analyzing the effect of increasing stimulation amplitude on phosphene size 8 . The prior study found that increasing pulse amplitude caused an increase in phosphene size at an average rate of 1.17 deg 2 /X threshold, with the slope varying slightly across electrodes (Fig. 4a). The study used data from nine electrodes in a single Argus I patient 8 . Our models predicted a similar pattern of amplitude modulating phosphene size with slope varying between 1.05 and 2.72 deg 2 /X threshold across participants and electrodes (Fig. 4b).

Discussion
Retinal prosthesis users report a variety of visual perception thresholds and phosphene shapes 10,26 . In this work, we developed statistical and patient-speci c eld-cable models to predict inter-patient and interelectrode variability. Such models, used individually or in combination, have the potential to expedite the process of " tting" a prosthesis for each patient, taking into account their particular anatomy and implant position.
Statistical models are a data-driven approach for predicting perceptual threshold. De Balthasar et al. previously investigated the effect of electrode-retina distance, retinal thickness, and impedance on perceptual threshold for six Argus I users 9 . The authors found a signi cant correlation between electroderetina distance and threshold for 1/6 participants, between retinal thickness and threshold for 3/6 participants, and between impedance and threshold for 5/6 participants 9 . Other prior studies have con rmed the correlation between electrode-retina distance and perceptual threshold, but have measured distance using the entire electrode array 10 or groups of four electrodes 11 . In this study, we examined the correlation of electrode-retina distance, retinal thickness, brotic tissue thickness, and electrode impedance with perceptual threshold for individual epiretinal disc electrodes using regression analysis.
One novel nding was that perceptual threshold increases systematically as retinal thickness decreases, with a xed effect across participants. This suggests that retinal thickness could be a proxy for retinal health, and the number of viable neurons in some cases. Therefore, we recommend pre-operative OCT to guide array placement toward healthy retinal regions. Our study does not take into account the effects of cystoid macular edema, which could dynamically change the thickness measurements and could potentially explain some of the variability seen between patients. In addition, our regression analysis corroborates prior studies that identi ed an effect of electrode-retina distance on perceptual threshold, but we found that the strength and slope of correlation varies signi cantly across individuals. Finally, we found that brotic tissue growth on the surface of the microelectrode array did not affect perceptual threshold. This supports the ndings of Rizzo et al., who found that 50% of Argus II users developed a brosis-like hyper re ective tissue at the array interface, but that those patients did not experience any deterioration in visual performance. 27 In general, data-driven models take an experimental dataset and t a lter to capture the relationship between input and output variables. This approach is effective because it can be quick and does not require prior knowledge of the system. However, data-driven models are essentially a black box, providing limited understanding of the underlying physiological mechanisms.
On the other hand, pairing three-dimensional bioelectric eld models with multi-compartment cable models aims to predict thresholds and phosphene patterns using rst principles of electrophysiology. These predictions emerge directly from imaging data in absence of any perceptual data. Our patientspeci c eld-cable models predicted activation thresholds that signi cantly correlated with perceptual thresholds for 4/5 participants. However, our results show that the slope of the correlation is not consistent across participants (Fig. 2). Individual device users have unique perceptual criteria for reporting phosphenes, and we do not know how many neurons must re to produce a phosphene. There are other cases in the eld of neuromodulation where each patient-speci c eld-cable model was t with its own optimal conductivity value to improve performance 28 . This indicates that some perceptual data collection may be required to calibrate the models, even with the patient-speci c eld-cable approach that could conceivably produce threshold predictions based solely on imaging data. Furthermore, our patientspeci c eld-cable models did not predict threshold better than the statistical models; the prediction accuracy was similar. On the other hand, they provided direct estimates of retinal activity that relates to the shape and size of phosphenes.
The impact of this work includes expanding our novel proof-of-concept study conducted in a single participant to a cohort of ve participants 19 while also re ning the original approach for greater e ciency. Our original work used both ultrasound and optical coherence tomography imaging to create a patient-speci c eld-cable model. We later found that ultrasound images, which we used to locate electrical ground, did not affect model predictions. Another improvement to our prior study was implementing a mammalian RGC model with dendrites that we optimized for run-time based on our sensitivity analysis, as opposed to a simpli ed, amphibian RGC model 29 Future work could apply this modelling framework to other retinal prostheses. In addition, we could apply optimization paradigms to these models to program patient-speci c stimulation settings. For example, we could evaluate multi-electrode paradigms that attempt to focalize phosphenes 6 . This virtual device programming session would help avoid manually testing a large parameter space in the clinic.
This study had several limitations. Current optical coherence tomography images cannot provide cellularscale imaging of the retina, and therefore cannot directly measure the number of viable RGCs 30 . Although we identi ed a signi cant, xed effect of retinal thickness on perceptual thresholds, we had to estimate the relationship between retinal thickness and neuron density. Removing RGCs from the model where the retina was less than 100 µm produced a good match with experimental data. However, the development of high-resolution retinal imaging that could directly measure cell count would improve our modelling framework. Secondly, our model did not include other inner retinal neurons (i.e., bipolar, amacrine, and horizontal cells). We assumed that epiretinal prostheses primarily cause direct RGC activation 1 . However, if this modelling framework were applied to intraretinal or subretinal prostheses, we should include a degenerate retinal network model 31,32 . A limitation of the patient-speci c eld-cable modelling approach is the time needed to build and run these models. To use them at a large-scale, we would recommend further automating our methodology.

Data collection
For this study, we obtained data from seven subjects implanted with the Argus II retinal prosthesis (Second Sight Medical Products Inc., Sylmar, CA). We recruited two subjects from the W.K. Kellogg Eye Center (University of Michigan, Ann Arbor, MI) and collected data for the purpose of this study (Clinical Trial ID: NCT03635645). We obtained informed consent following approval from the University of Michigan's Institutional Review Board. The study adhered to the tenets of the Declaration of Helsinki and national regulations for medical device clinical trials. Five subjects were patients at the Cole Eye Institute (Cleveland Clinic, Cleveland, OH). Previously collected de-identi ed data was extracted from medical records, as speci ed by a data transfer and use agreement.
We obtained optical coherence tomography (OCT) scans of the implanted eye for each study participant. For UM01 and UM02, scans were taken using a Heidelberg Engineering Spectralis system. Images spanned 30° x 25° of the visual eld, using 62 B-scans. Each B-scan was 768 pixels (8.8 mm) by 496 pixels (1.9 mm) and the scan-to-scan spacing was 122 µm. For the remaining participants, scans were taken using a Zeiss Cirrus system. Images spanned 20° x 20° of the visual eld, using 128 B-scans. Each B-scan was 400 pixels (6 mm) by 120 pixels (2 mm) and the scan-to-scan spacing was 46.875 µm.
First, we used the OCT scans to make measurements of electrode-retina distance, retinal thickness, and brotic tissue thickness for all visible electrodes. The Argus II implant has sixty electrodes arranged in a 6x10 grid, with the standard labels A1-F10 (see Fig. 8a). We made measurements using Mimics Research Version 24.0 (Materialise NV, Leuven, Belgium). We identi ed individual electrodes using the coronal view (fundus view), and made corresponding measurements using the axial view (B-scans). The opaque platinum electrodes occlude the scanning light, creating a dark shadow on the underlying tissue. This re ection artifact con rmed the exact electrode locations on the B-scans, and we made measurements at the center of the largest shadow obtained by a B-scan. Since the OCTs could not precisely center every electrode in line with a scan pattern, our measurements occurred at the center of the chord that was closest to the full electrode diameter. We measured the retinal thickness as the distance from the inner to outer boundary, without subdividing retinal layers. Figure 5 shows example measurements for two electrodes.
We also obtained electrode impedance and perceptual threshold data for each participant. We used the Argus II Clinician Fitting System (CFS) to measure the impedance of each electrode. We used the Hybrid Threshold program on the CFS to determine perceptual threshold for individual electrodes, as described in prior work 19 . Perceptual threshold is the current amplitude at which the participant sees a percept 50% of the time, based on a Weibull function sigmoid curve t to "yes"/"no" responses. Electric stimuli were biphasic, cathode-rst current pulses with 0.45 msec pulse width and 20 Hz frequency. We used a train of ve identical pulses for each trial (250 msec duration). Perceptual threshold collection occurred over a 1-2 day period for each participant.

Statistical models for perceptual threshold
We used regression analysis to create statistical models for perceptual threshold. An electrode was included in the analysis if it was visible on at least one B-scan and the perceptual threshold was less than 677 µA. We provide a summary of excluded electrodes in Supplementary Table S1. Explanatory variables were electrode-retina distance (mm), retinal thickness (mm), brotic tissue thickness (mm), and electrode impedance (kΩ). We t individual linear regression models to nd the in uence of each explanatory variable on perceptual threshold for each participant. We t a multiple regression model to nd the cumulative in uence of all explanatory variables on perceptual threshold for each participant. In each case, we calculated the coe cient of determination (R 2 ) and p-value. We used a signi cance level (α) of 0.05. We performed statistical analysis using R version 4.3.0 (R Core Team, 2023).
Additionally, we used linear mixed models to determine the in uence of explanatory variables on perceptual threshold across participants. We t a random intercept and slope model for each explanatory variable using the 'lme4' package to explore how variation in perceptual threshold depended on xed effects across individuals and on individual-level random effects. Signi cance of xed effect was determined using Satterthwaite's degrees of freedom. Signi cance of random slope was determined by likelihood ratio test.
Patient-speci c eld-cable models for perceptual threshold Our methodology for creating patient-speci c eld-cable models was adapted from our prior proof-ofconcept study, conducted in a single participant 19 . We created models of the implanted eye for participants with at least one signi cant explanatory variable identi ed during regression analysis. As a result, we excluded participants CL01 and CL02 from this portion of the study, and we built ve patientspeci c eld-cable models.
We segmented OCT scans in Mimics Research Version 24.0 (Materialise NV, Leuven, Belgium) using the multiple slice edit tool. Prior to segmentation, we corrected the jitter in scans obtained with the Zeiss Cirrus system using the image processing toolbox in MATLAB R2021a (MathWorks Inc., Natick, MA, USA).
We converted the segmented images into a nite element mesh (FEM) using the 3-Matic Research Version 16.0 (Materialise NV, Leuven, Belgium) (Fig. 6b). Pre-processing included wrapping each domain with a gap closing distance of 25 µm. Sixty circular electrode surfaces were created on the MEA surface by using cylinders to create intersection curves. Cylinders were 200 µm in diameter and spaced at 525 µm pitch and were aligned using the electrode shadow artifacts. The resulting electrode surfaces created by the intersection curve operation are shown in Fig. 6c.
We created a processing pipeline to re ne each patient-speci c FEM. To limit edge effects, we extruded the edges of the retina and choroid to create a rectangular model with 25 x 17 mm dimensions. We created a vitreous domain that extended 18 mm above the retinal surface. We performed an adaptive remesh of all surfaces with a maximum triangle edge length of 0.05 mm for the MEA, and a maximum triangle edge length of 0.1 mm for the retina, choroid, and vitreous. We built a grid-based non-manifold assembly (grid resolution: 0.01 mm) to align the triangle nodes of adjacent surfaces. Finally, we calculated a tetrahedral volume mesh for the entire model. The number of tetrahedral volume elements for each patient-speci c FEM was between 3.3 to 4.0 million.
We conducted nite element analysis in COMSOL Multiphysics Version 5.6 (Stockholm, Sweden) using the AC/DC electric currents (ec) module, following the methods described in our previous work 19 . We represented the active electrode as a surface current terminal (1A) and assigned a oating potential boundary condition to inactive electrodes 23 . We designated the outer boundary of the choroid as the electrical ground (0V). We assigned a tissue conductivity value to each domain (σ choroid = 0.503 S/m 33 , σ retina = 0.100 S/m 24,34 , σ vitreous = 1.5 S/m 33 , σ platinum = 9.43x10 6 S/m 23 , σ brosis = 0.15 S/m 35 ). We used a contact impedance condition to model the thin, resistive retinal pigment epithelium membrane at the boundary between the retina and choroid (thickness = 10 µm, σ = 0.001 S/m 24,36 ). We modelled the MEA substrate as a perfect insulator (σ = 0 S/m). We used a quasi-static solver to calculate the electric potential distribution throughout the nite element mesh by solving Laplace's equation (Fig. 7).
To predict the neural retina's response to stimulation, we populated the retinal domain of the FEM with multi-compartment cable models of retinal ganglion cells. For each electrode in our FEM, we uniformly distributed 250 somas within a 700 µm radius beneath the electrode using Lloyd's algorithm 37 . We calculated axon trajectories using the equations de ned by Jansonius et al. 38 (Fig. 8a). We positioned each soma 55 µm below the retinal surface mesh and positioned each axon 15 µm beneath the retinal surface mesh (Fig. 8b). Figure 8c shows the coordinates of an example target RGC population for one electrode.
The cable equations for RGC membrane dynamics are described in detail in our prior publication, and our biophysical model is open-source on GitHub 29 . Brie y, we determined the extracellular potential at the center of each neuron compartment from the FEM solution. We scaled these extracellular potentials by the stimulus waveform, and applied them to the neurons, integrating over time to calculate the membrane voltage response. To match clinical experiments, electric stimuli were biphasic, cathode-rst current pulses with 0.45 msec pulse width. We found the action potential threshold for each individual neuron using a bisection algorithm with a convergence of 0.1 µA.
We found the absolute minimum action potential threshold for any neuron in each electrode's target RGC population. We used the same electrode inclusion criteria as described for the statistical model. By using a uniform distribution of neurons, our patient-speci c eld-cable models did not capture the in uence of retinal thickness on perceptual threshold. To remedy this, we removed neurons in regions where the retina had degenerated to a thickness less than 100 µm 39 . We re-computed the absolute minimum action potential threshold in each electrode's target RGC population. We compared the model-predicted thresholds to each participant's perceptual threshold using linear regression.

Phosphene predictions
We compared the patient-speci c eld-cable model predictions of retinal activation to a prior study analyzing the effect of increasing stimulation amplitude on phosphene size. 8 Nanduri et al. collected phosphene drawings in response to single-electrode stimulation from nine epiretinal disc electrodes. 8 In concordance with our study, electric stimuli were biphasic, charge-balanced, 0.45 msec/phase, cathode-rst pulse trains that were 500 msec in duration. 8 Stimulus amplitude was modulated between 1.2x and 6x threshold, while holding frequency constant at 20 Hz. 8 The study found that increasing pulse amplitude caused an increase in phosphene size at an average rate of 1.17 deg 2 /X threshold (range of slopes: 0.73-1.92). 8 To replicate this study in our patient-speci c eld-cable models, we chose a single electrode from each model. We calculated the retinal area activated as stimulus amplitude was modulated between 1.2x and 6x threshold. We calculated retinal area (mm 2 ) using the "shapely.MultiPoint" tool in Python to t a convex hull around activated cell bodies (Supplementary Figure S1). Retinal area was converted to degrees of visual angle using a conversion factor of 288 µm per degree. 40 Declarations Data Availability The datasets generated during and/or analyzed during the current study are not publicly available out of concern for patient privacy, as data from medical records is considered sensitive, but de-identi ed datasets are available from the corresponding author on reasonable request.

Code Availability
The computer code used to model the extracellular stimulation of retinal ganglion cells using cable equations is available at: https://github.com/Kathleen-Kish/Retinal_Ganglion_Cell  Linear regression analysis comparing visual perception thresholds with action potential thresholds in the patient-speci c eld-cable models. Thresholds predicted by the patient-speci c eld-cable models for four electrodes (dashed black line).
Retinal ganglion cell somas are shown as black dots, and colored contours show the distribution of action potential thresholds beneath each electrode. The absolute minimum action potential threshold for each electrode is noted.  Optical coherence tomography scans with measurements. The top image shows a B-scan for participant CL05, with measurements for electrode-retina distance (0.31 mm) and overall retinal thickness (0.19 mm) for electrode C6. The bottom image shows a B-scan for participant UM01, with measurements for brotic tissue thickness (0.08 mm) and retinal thickness (0.33 mm) for electrode E7. In this case, the microelectrode array is opposed to the retina and the electrode-retina distance is equal to the brotic tissue thickness.   Example of electric potential distribution calculated in the patient-speci c FEM for electrode D5, participant CL03.

Figure 8
Neuron placement paradigm (a) The x-y coordinates of a target retinal ganglion cell population beneath electrode A3 for participant UM01. Somas were distributed in a 700-µm circular region beneath the electrode, and axon x-y coordinates were calculated using equations de ned by Jansonius et al. 38

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.