Computational modeling of endovascular peripheral nerve stimulation using a stent-mounted electrode array

Objective. Endovascular neuromodulation has attracted substantial interest in recent years as a minimally invasive approach to treat neurological disorders. In this study, we investigated with a computational model the feasibility of stimulating peripheral nerves with an endovascular stent-mounted electrode array. Approach. Anatomically realistic FEM models were constructed for the pudendal and vagal neurovascular bundles. The electromagnetic fields generated from electrical stimuli were computed using Sim4Life NEURON models to predict dynamic axonal responses. Main results. The models predict that the stimulation thresholds of the endovascular stent-electrode array configurations tested are comparable to that of ring electrodes and are dependent on the inter-electrode distance and orientation of the device. Arranging multiple electrodes along the longitudinal axis of the nerve lowers surface charge density without sacrificing axon recruitment, whereas arranging electrodes along the circumference of the blood vessel reduces the risk of misalignment but lowers axon recruitment. Significance. Overall, this study predicts that the endovascular stent-electrode array is a feasible stimulation option for peripheral nerves, and the electrode array can be flexibly optimized to achieve the lowest stimulation threshold.


Introduction
In recent years, there has been significant progress in the development of peripheral nerve interfaces (PNIs). Interfacing with peripheral nerves using various electrical neurostimulation modalities has enabled restoration of function lost to injury, such as restoring sensation and motion of an amputated or paralyzed limb [1,2], modulating bladder overactivity secondary to spinal cord injury [3,4], and supporting respiratory function [5]. Despite successful therapeutic attempts, there are still obstacles to be overcome for designing implantable PNIs intended for chronic use.
The design of a PNI faces a trade-off between invasiveness and stimulation selectivity or recording resolution. Extraneural PNIs are positioned outside nerves and are suitable for chronic applications that require long-term interface stability. Conventional extraneural designs include nerve cuffs, book electrodes, helical cuffs, and flat interface neural electrodes [6][7][8][9]. Chronic implantation of such electrodes may inflict damage to the nerve through constriction, traction, or relative motion, leading to morphological changes such as thinning of myelin, perineural thickening, proliferation of endoneural connective tissue, and reduced density of myelinated axons [10,11]. Chronic inflammation and fibrotic encapsulation of PNI devices are also commonly associated with extraneural PNIs [12][13][14]. Fibrous tissue develops between the contact surfaces of the nerve epineurium and PNI. Consequently, the electrodes essentially migrate away from the target nerve, circumventing the spatial advantage of the PNI designed to be in direct contact with the nerve trunk. Fibrotic encapsulation also increases tissue resistivity around the implanted electrodes, which may reduce the stimulation efficiency or the quality of recorded signals [15,16]. The fibrous tissue between a chronically implanted cuff and nerve surface can reach a thickness of 500 µm [10], which may significantly alter the electric field generated by the electrodes. Endovascular PNIs have been recently explored as a less invasive alternative approach to deliver stimuli to or record from target nerves without the need of open surgery or physical contact with the neural target. Preclinical and clinical attempts of endovascular peripheral neuromodulation include the induction of bradycardia through stimulation of the vagus nerve [17,18], treatment of sleep apnea through phrenic nerve stimulation [19], and induction of muscle contractions through femoral nerve stimulation [20]. The most recent advancements in endovascular PNI design include the Stentrode TM (Synchron, Brooklyn, NY), a self-expanding stent with an incorporated electrode array [15,[21][22][23]. The radially expanding stent provides fixation of the electrode array on the blood vessel wall, preventing electrode migration and facilitates vessel wall incorporation of the device. It has been previously explored as a minimally invasive approach for intracranial electrical recording and brain stimulation [24]. Acute and chronic pre-clinical studies have shown the efficacy and safety of a stentmounted electrode array implanted in a cortical vein in a sheep model for stimulation and recording [15,21,22], and the first-in-human trial of a motor neuroprosthesis using the stent-electrode array is underway [23]. However, to date, attempts have not been made to expand the application of the endovascular stentelectrode array to the peripheral nervous system. Similar to other PNIs, applications of the stentelectrode array are expected to be constrained by a number of factors, including the relative position of the target blood vessel, neural anatomy and electrode configuration. Computational modeling is a powerful tool for evaluating the theoretical effects of such factors on neuromodulation efficacy. A previous modeling framework developed by Teplitzky et al [25] investigated the feasibility of endovascular deep brain stimulation with ring electrodes and guidewire electrodes. Our previous work has also evaluated the effect of the relative positions of the nerve and blood vessel on endovascular neurostimulation efficacy in an idealized pelvic nerve model [26]. However, little is known about how the electrode geometry of the stent-electrode PNI quantitively affects neurostimulation efficacy.
In this study, we developed a computational modeling framework to investigate the effect of interelectrode distance, electrode orientation, shape and configuration on axon fiber recruitment of two known peripheral neuromodulation targets: (1) the pudendal nerve for restoring bladder control in people with spinal cord injury [3,4] and (2) the cervical vagus nerve for treating diseases such as epilepsy, depression, and arthritis [27]. Based on available anatomical data, two approaches to reconstruct realistic tissue volumes were demonstrated. A finite element model was used to compute the electromagnetic field induced in the tissue volume conductors by an electrical stimulus. The simulated electromagnetic field was coupled with multi-compartment cable models of axons to predict the axonal dynamic responses and to estimate stimulation thresholds. Understanding the impact of electrode array configuration on fiber recruitment will be an important theoretical basis for the optimization of endovascular PNI design.

Methods
For the pudendal nerve and vagus nerve, a finite element method (FEM)-NEURON hybrid model was constructed using Sim4Life (Version 6.2, ZMT Zürich MedTech AG, Switzerland). A tissue volume conductor model based on the FEM was solved for the electromagnetic field applied by electrical stimuli delivered by the endovascular electrodes, which was interpolated along the axons and passed on to multicompartmental neuron models to simulate axonal responses.

FEM tissue conductor model
A cross-section of the key anatomical structures at the level of the proximal end of the pudendal canal (Alcock's canal) was segmented (figure 1(a)) [28]. The compound pudendal nerve within the pudendal canal was chosen as the stimulation target due to its proximity to the blood vessels within the canal, which may reduce the current threshold required to achieve adequate axon recruitment. Anatomically realistic reconstruction of the fascicular structure was performed based on histological sections of a human female pudendal nerve at the selected level [29]. An image of the histological section showing the fascicular map was imported into Sim4Life and segmented manually ( figure 1(b)). Since the local tissue anatomy near the pudendal neurovascular bundle remains relatively constant throughout the length of the canal, the segmented geometry was extruded longitudinally to create 3D tissue volumes.
The FEM tissue conductor model consisted of several key anatomical features observed near the pudendal nerve. The model included a thin layer of connective tissue surrounding the pudendal neurovascular bundle, vessel walls of the internal pudendal artery and vein and the blood flowing within, the ischiorectal fossa filled with fatty tissue, and the obturator internus muscle ( figure 1(c)). The disk electrodes were modeled based on the original Stentrode TM (Synchron Inc., Brooklyn, NY) design [21] with a thickness of 50 µm and a diameter of 750 µm, which conformed to the inner surface of the pudendal artery vessel wall. The struts of the Stentrode TM design were not included in the model as they are electrically insulated from the electrode array and, therefore, assumed to have negligible impact on  [29]). The fascicle branching pattern was determined by identifying the distal branches and tracing proximally. DGN: dorsal genital nerves; PerN: perineal nerves. (c) A cross-section of the cylindrical FEM tissue model with the pudendal neurovascular bundle embedded in a layer of connective tissue, surrounded by muscle and fatty tissue. (d) Extruded 3D model showing an example of bipolar disk electrode placement within the artery. (a) Adapted from [28], Copyright (2020), with permission from Elsevier. (b) Adapted from [29], with permission from Springer Nature.  [30][31][32][33].
To simplify tissue geometry and reduce computation time, the modeling domain was limited to a cylinder with a finite length and radius. A zero-flux boundary condition was set at the outer boundary of the modeling domain such that electric flux perpendicular to the outer surface was zero. This Neumann boundary condition simulates an infinite system. For all instances tested in this study, Dirichlet boundary conditions were set on the electrode surfaces in a bipolar configuration, with the active electrodes set to 1 V and the return electrodes set to 0 V as ground ( figure 1(d)). The total current flux passing though the active-return electrode pairs was calculated and the computed potential field was scaled accordingly such that the total stimulation current was held constant at 1 mA. A Dirichlet boundary condition was used instead of a Neumann boundary condition to preserve a realistic representation of the current density field around an insulated electrode by assuming non-uniform current through its surfaces. This adjusted potential field was to be further linearly scaled in the NEURON simulation to determine the activation threshold of each axon. The electrodes were modeled as perfect electric conductors with zero resistivity so that tangential electric flux through electrode surfaces was set to zero and the solver did not solve for the internal electric field within the electrodes. An unstructured mesh with tetrahedral elements was created using the built-in multidomain mesh tool of Sim4Life. Mesh element edge length was set to be finer within and near the fascicles and the electrodes, where the potential gradient was expected to be high, and coarser elsewhere.

Outer boundary size sensitivity
Reducing mesh element count increases the error incurred by approximating an infinite tissue volume with a finite simulation domain. The error introduced by the finite volume was investigated by evaluating the differences in simulated potential field within the nerve trunk and in computed axonal responses to the same stimulus in response to changes in the length and radius of the simulation domain. Lengths of 15,20,30,40,50, and 55 mm and radii of 15, 20, and 30 mm of the cylindrical volume were tested, and the models were meshed with the same tissue-specific target element sizes. The positions and sizes of the electrodes remained invariant during these tests and were positioned symmetric to the central transverse plane of the model. Under these conditions, the simulated potential field was expected to be symmetric to the central transverse plane.
To compare the potential fields computed at nodes of different volumetric meshes, an approach modified from a previous modeling study [32] was adopted. Computed electrical potential V e at mesh nodes on seven parallel transverse planes within the enclosed volume of the nerve trunk were isolated. The planes were sampled axially starting from the central transverse plane along the nerve length, with a distance of 1 mm between two adjacent planes, so that the sampled planes roughly encompassed half of the volume of the shortest model length tested, which was 15 mm. V e values on each plane were interpolated with a uniform mesh grid (figures 2(a) and (b)) and the differences between instances of varying model dimensions were quantified by root mean square error (RMSE) [32], where N is the total node number of the mesh grid and V i is the interpolated potential value on node i of instance 1 or 2. Additionally, as the position, fiber diameter and neuron model of each axon remained consistent in different instances, the activation threshold and the location of the first spike of the same axon were also compared between different instances. Increasing the outer dimensions of the simulating domain lowered the finite boundary condition error, but this improvement became insignificant after the outer boundary surpassed a minimal sufficient size. A minimal sufficient outer boundary dimension was defined as a threshold beyond which further increase in the dimension would result in a normalized RMSE (NRMSE, normalized by the span of the computed potential values within the nerve trunk) lower than 5% in terms of individual axon activation threshold, location of first spike, and interpolated potential field on the selected planes. A boundary length of 40 mm yielded a NRMSE less than 5% in all three aspects considered compared to the 50 mm and 55 mm instances and was selected as the minimal sufficient model height (figure 2(c)). Similarly, a boundary radius of 20 mm was selected as the minimal sufficient model radius. The minimal length and radius were adopted in all further instances tested.

Modeling of dynamic axonal responses
Axon seed points were generated randomly from a uniform distribution on a cross-section of the fascicles and extruded to parallel splines. Since we have been unable to identify published data on axon diameters and fascicular axon counts for human pudendal nerve, the axon count for each fascicle was set to be proportional to the surface area of its cross-section so that the smallest fascicle contained 10 axons, resulting in a total population of 431 axons. A bimodal Gaussian distribution with means µ 1 = 4.786 µm and µ 2 = 9.427 µm and standard deviations σ 1 = 0.0714 µm and σ 2 = 0.143 µm was used to approximate the physiological myelinated fiber diameter distribution of the pudendal nerve, based on the myelinated fiber count of the human sural nerve and anterior tibial nerve [32,34] (figure 3(a)). A multi-compartmental cable model of myelinated mammalian nerve fiber, the Sweeney model, was used to simulate the dynamic axonal responses of myelinated afferent fibers. This model assumes perfectly insulating internodal myelin and represents nodes of Ranvier with voltage-dependent sodium channels and voltage-independent leakage channels [35]. The number of the resulting node-internode compartments ranged from 31 to 141 depending on fiber diameter.
To investigate whether endovascular electrical stimulation can recruit the small and unmyelinated C-fibers, a skew Gaussian distribution with mean µ = 0.790 µm, standard deviation σ = 0.613 µm and skewness = 7 was estimated from the unmyelinated fiber size of the cat medial nerve and posterior articular nerve [36] and was used to represent the unmyelinated sensory fiber diameter distribution ( figure 3(b)). The Tigerholm model originally developed for modeling activity-dependent changes of action potential conduction in human C-fibers [37] was used to simulate the axonal responses of the unmyelinated sensory fibers. This model assumes a constant section length of 50 µm and as a result each axon was modeled as 798 sections. Axons were randomly assigned to the myelinated or unmyelinated population so that the incidence of either group was close to 50%.
Sim4Life interpolates the extracellular potential field along each axon spline and directly passes the solution to its built-in NEURON solver (Yale University, CT, USA). A modulating pulse is applied to the static potential field and linearly scales it to vary with time according to the stimulus waveform. The stimulus was a biphasic and symmetric rectangular wave with a pulse width of 300 µs and no interphase delay. The titration function of Sim4Life was used to determine the activation threshold in terms of stimulus amplitude of every axon. Titration introduces a varying titration factor in addition to the modulating pulse to linearly scale the potential field sensed by neurons iteratively until an action potential is detected.

Effect of stent-electrode-array design
As endovascular electrodes are constantly immersed in blood, a major challenge limiting their stimulation efficiency is the less resistive path between electrodes through blood. Blood can potentially shunt away most of the current from the neural tissue nearby and, therefore, result in fewer axons being activated. To account for this issue, endovascular disk electrodes usually have an insulating backing layer to prevent the electrodes from making direct contact with blood. Under the same conditions, a pair of disk electrodes were modeled with or without an insulative silicone backing (conductivity 10 −12 S m −1 ) for comparison. For the insulated electrodes, a 50 µm layer of silicone was added to encase the electrodes, only exposing the surfaces in contact with the vessel wall. The electrodes were modeled to be 8 mm apart along the longitudinal axis of the blood vessel with one serving as the stimulating electrode and the other as ground. In addition, the effect of a thin fluid layer between the electrode surface and tissue was investigated by moving the electrodes 50 µm away from the blood vessel wall.
Previous modeling work has shown that longer inter-electrode distance may reduce current shunted away from the target neural tissue by other tissues near the electrodes [38,39]. To investigate the effect of inter-electrode distance on endovascular stimulation efficiency, a pair of disk electrodes at from 2 mm to 14 mm at 2 mm increments were modeled. Three electrode diameters, 500 µm, 750 µm and 1 mm, were evaluated. A pair of 750 µm disk electrodes directly in contact with the medial side of the pudendal nerve was also investigated in comparison to the endovascular electrodes. For all single-pair electrodes tested, the active and return electrodes were placed symmetrical to the middle transverse plane of the blood vessel. Endovascular electrode orientation in relation to the target neural tissue also has an influence on threshold current amplitude in vivo [22]. The effect of electrode orientation on endovascular stimulation was evaluated by rotating a pair of electrodes at 8 mm distance relative to the longitudinal axis of the blood vessel.
To investigate the effect of electrode geometry and configurations on activation threshold, nine configurations were evaluated in this study as shown in figure 6(a). All configurations were bipolar, and the disk electrodes (750 µm diameter) were divided into two groups of opposite polarity by the middle plane. In addition to the basic single electrode pair (configuration 1), two strategies of disk electrode placement were tested. Electrodes of the same polarity were either positioned vertically along the longitudinal axis of the blood vessel (configurations 2 and 3) or horizontally and around the blood vessel circumference (configurations 4-8). The single pair, vertical configurations and the midline of the horizontal configurations were oriented in the direction of the target nerve. A pair of 1.5 mm wide ring electrodes at 8 mm distance was also evaluated for comparison.

Vagus nerve model 2.2.1. FEM tissue conductor model
A computable anatomical human phantom based on the Yoon Sun model (V3.1) of the virtual population anatomical models (IT'IS Foundation [40]) was used to construct a FEM tissue volume conductor model including the vagus nerve. A cylindrical volume with diameter 128 mm and length 110 mm, mimicking the dimensions of the neck, was created. This volume spanned the vertebral levels T1 to C1 and was used as the model outer boundary. From the human phantom, the internal and external jugular veins, the carotid artery, the cervical vagus nerve, and muscles of the neck surrounding the vagal neurovascular bundle were isolated. The volumes of these anatomical structures bounded by the cylindrical outer boundary were used as tissue volume conductors (figure 4(a)). A previous anatomical study on human cervical vagus nerve morphology has reported the mean effective surface area and fascicle number of the right nerve to be significantly larger than the left nerve [41]. To reduce the computation load, only the tissue volumes of the right lateral half of the neck were included in the model. Inclusion of the vertebra and cerebrospinal fluid or tissues of the left half of the neck only led to less than 1% change in the axonal activation thresholds; therefore, these additional tissues were excluded from the FEM simulation.
To preserve the realistic anatomy of the branching cervical vagus nerve, algorithms were written to map the parent-and child-branch relationship and to populate the branching nerve trunk volume with axoncontaining fascicles. First, with Sim4Life, starting with the proximal, unbranched end of the main vagus nerve trunk, transverse cross-sections of the anatomical structure were taken at a small uniform interval and the centroid coordinate of each cross-section was calculated. On each subsequent transverse plane, new centroids were connected to the closest centroids of the previous plane to create a continuous axial centerline of the tissue volume. Multiple centroids of new child branches below the branching point could be linked to the parent branch using this method. The vertices defining each nerve cross-section contour and the completed mapping of parent-child branches was recorded as a text file and imported into MAT-LAB. The original anatomical blood vessel models of the Yoon Sum model [40] did not separate blood vessel walls from the blood volumes, despite their different conductive properties. Therefore, the same algorithm was also applied on blood vessel models, and each transverse cross-section was radially shrunk by 0.3 mm to create a contour of the vessel wall tissue. With Sim4Life, sequences of cross-sections were connected to create separate blood and blood vessel volumes using the Skin Wires function ( figure 4(b)).
Based on previous morphology studies [41,42], a total of seven circular fascicles were modeled, including three large fascicles with diameter 0.24 mm and four small fascicles with diameter 0.125 mm. Using MATLAB (Version R2021a, MathWorks, MA, US), a fascicle pattern on the terminal cross-section of the main nerve trunk was generated by randomly positioning fascicle seed points within the bound of the nerve trunk while avoiding any overlapping or touching of fascicle contours. The consistency of the fascicle pattern was kept as stable as possible across nerve trunk cross-sections by translating the seed points to account for the change in nerve trunk centroid coordinates, and minimal reshaping was performed algorithmically. On each subsequent nerve crosssection, the locations of fascicle seed points were iteratively checked and optimized to avoid fascicles breaching the boundary of the nerve contour or overlapping. In each iteration, if a fascicle contour was detected not to be fully within the nerve contour, its seed point would be shifted closer to the centroid. This would likely cause fascicles to overlap so an overlap check was performed subsequently. Any pair of partially overlapping fascicles were moved slightly away from each other. The optimization continued until no breaches or overlaps could be detected and the same process was carried out on the following cross-section. At branching points, fascicles of the parent branch were assigned to child branches according to the smallest surface areas of the crosssections along each child branch to avoid clogging a small branch with more fascicles than it could accommodate. The coordinates of fascicle seed points on each cross-section and fascicle parameters were imported back into Sim4Life to create fascicle tissue volumes as splines thickened to the specified diameters ( figure 4(c)). Tissue conductivities in table 1 were also applied to the tissue conductors of this model.
Out of the disk electrode arrangements investigated with the pudendal nerve model, the four configurations that achieved the highest activation level by the same stimulus intensity were adapted to the vagus nerve model and evaluated. A pair of ring electrodes with a length of 1.5 mm was also tested and compared with the disk electrode configurations. All endovascular electrodes were conformed to the inner wall of the jugular vein to ensure proper contact and mesh nesting. A 50 µm-thick layer of silicone insulator was added to engulf all electrodes so that only the electrode surface in contact with the blood vessel was exposed. Boundary conditions similar to the aforementioned were used for the vagus nerve model. A zero-flux boundary condition was applied to the outermost surface of the modeling domain, and bipolar Dirichlet boundary condition was set on electrode surfaces as shown in figure 8.

Modeling of dynamic axonal responses
A total of 375 axons were populated in the fascicle volumes as randomly translated clones of the fascicle centerline bound by the fascicle volume. Large fascicles were populated with 95 axons and small ones with 23 axons to maintain a uniform axon density. A bimodal Gaussian distribution with means µ 1 = 4 µm and µ 2 = 8.5 µm and standard deviations σ 1 = 0.06 µm and σ 2 = 0.3 µm was used to approximate the physiological myelinated fiber diameter distribution of the human cervical vagus nerve [43] (figure 3(c)). The Sweeney model was used to simulate the axonal responses of the vagal myelinated fibers.

Results
For all instances tested in this study, no action potential was successfully evoked in the simulated unmyelinated axons with stimulus intensities up to 20 mA. Therefore, all the activation profiles and activation thresholds reported were yielded from the myelinated axons.

Effect of insulating backing
Simulated neuronal response showed that the insulation backing increased axon recruitment consistently across all stimulation amplitude tested (1-10 mA) with an average increase of 4.0%. In contrast, the inclusion of a thin layer of blood reduced axon activation level by roughly 1.6% ( figure 5(a)). Figure 5(b) shows the effect of inter-electrode distance and electrode size on the population mean threshold current amplitude required to activate each axon. Increasing electrode size resulted in a slight reduction (<5%) in mean activation threshold. Electrodes directly in contact with the target nerve showed higher axon recruitment efficiency at all inter-electrode distances, lowering mean activation threshold by roughly 0.35 mA compared to the endovascular stimulation at sufficiently large interelectrode distances over 8 mm. For all endovascular electrode sizes tested, increasing the interelectrode distance initially caused a rapid drop in mean activation threshold, followed by a slow and slight rise. The minimum activation threshold and thus highest stimulation efficiency occurred at 8-10 mm. Figure 5(c) illustrates the predicted activation profile of the pudendal nerve as the stimulation amplitude increased up to 10 mA. The results predicted that activation of above 85% of the total axon population could be achieved at 10 mA with an inter-electrode distance of 4 mm and above. At 5 mA, the activation level reached above 70% at 6 mm and above.

Effect of electrode size, orientation and inter-electrode distance
The predicted activation level (figure 5(d)) decreased as the electrode-target relative angle increased. At 0 • -15 • , the electrode pair faced the target nerve and most current was directed towards the neural tissue. Orienting the electrodes away from the neural tissue at 135 • -180 • decreased activation level by more than 45% at 10 mA.

Electrode configurations
The resulting mean activation thresholds and axonal activation profiles from the nine electrode configurations evaluated are shown in figure 6(b). Configuration 1, the single pair of disk electrodes 8 mm apart, achieved the highest level of axon recruitment. Additionally, configuration 2, a pair of bar electrodes that resemble the bounding box of the electrodes of configuration 4 (vertical far) with an effective interelectrode distance (distance between the centroids of the active and return electrode groups) of 10 mm, was also tested. Configuration 2 yielded <5% increase in axonal activation level compared to configuration 4, both only slightly underperforming configuration 1. The horizontally arranged electrode configurations (configuration 5-9) showed progressively declined recruitment as the angular distance between adjacent electrodes increased. All disk electrode configurations outperformed configuration 10, the ring electrode pair.

Fascicle selectivity
The single pair electrode configuration was used to investigate the effect of distance on fascicular   the blood vessel lumen as well as axon diameter ( figure 7(b)), well in accordance with the physiological recruitment order observed in electrical stimulation of peripheral nerves [44]. Figure 7(c) shows the predicted activation profiles of the four distal nerve branches of the pudendal nerve, namely the three perineal nerves (PerNs) and the dorsal genital nerve (DGN). The DGN, at the shortest distance from the artery, yielded the highest activation overall. At a moderate stimulation intensity of 3 mA, 89.6% of DGN axons were predicted to be activated, while 29.1% of all axons that comprised the PerNs were predicted to be activated. At 5.5 mA, the DGN was predicted to be fully activated, whereas the PerNs reached a 53.1% activation level. The Kruskal-Wallis test was used to evaluate changes in the normalized activation threshold of every axon of each pudendal nerve branch caused by varying the electrode configuration. The test yielded a p-value of 0.02 for DGN but larger than 0.96 for all other branches across all configurations. Pearson's correlation test showed strong positive correlation between the percentage of activated DGN axons out of all activated axons and the axon population mean activation threshold across all electrode configurations (ρ > 0.93, p < 0.01). Therefore, electrode configurations with lower axon recruitment increased selectivity towards the DGN branch.

Vagus nerve model 3.2.1. Electrode configurations
All disk electrode configurations investigated with the vagus model were positioned in the jugular vein in a way that the longitudinal midline of the electrodes was as close to the vagus nerve as possible. In contrast to the pudendal nerve, configuration 4, the six electrodes horizontally arranged with 60 • angular distance between adjacent electrodes ( figure 8(a)), outperformed all other disk electrode configurations in terms of mean activation threshold and axonal activation levels by roughly 16%, with a population mean threshold of 14.8 mA and 37.1% of all myelinated axons recruited at 10 mA. Assuming that myelinated axons in this model with diameters ranging from 2 to 4 mm were B-fibers, a maximum of 12.7% of B-fibers was activated at 10 mA. The other three disk electrode configurations showed comparable performance. The ring electrodes yielded the lowest level of activation, recruiting only 6.1% of the myelinated axon population at 10 mA.

Fascicle selectivity
A random fascicular map consisting of seven fascicles was generated for the vagus nerve model, as shown in figure 9(b). Fascicles 1-3 were large fascicles with diameters of 0.24 mm whereas fascicles 4-7 had diameters of 0.125 mm. While the electrode pair was positioned above the first branching point at position c in figure 9(a), fascicles 5, 3, and 1 were at the shortest distance from the adjacent internal jugular vein and showed significantly higher fiber recruitment than the rest of the fascicles. Fascicle 5, a small fascicle with a thinner layer of perineurium tissue, showed the highest activation level at 95.6% with a stimulus amplitude of 10 mA. Inferior to the first cervical level branching point, fascicles 3 and 4 extended inferiorly and made up a lateral branch that coursed along the internal jugular vein at a distance less than 2 mm, whereas the other fascicles coursed inferiorly as medial branches at 8-10 mm from the internal jugular vein. At position d, the lateral fascicles 3 and 4 both achieved 100% fiber recruitment at small current amplitude of 1.5 mA; in contrast, no other fascicle showed any activation at the same stimulus intensity. At the further inferior position e, the electrodes were able to recruit more than 86% of the axons of fascicles 3 and 4 at 5 mA, while no other fascicles showed any sign of activation with a stimulus below 9.5 mA.

Discussion
This study provides a theoretical framework for using a FEM-NEURON hybrid model to predict the feasibility of endovascular stimulation of peripheral nerves. The impact of endovascular electrode configuration on axon fiber recruitment has been quantitively evaluated with the computational models. The internal pudendal artery and jugular vein were selected as suitable PNI implantation targets for the pudendal nerve and vagus nerve, respectively. The internal pudendal artery was used since it was closer to the pudendal nerve than the pudendal vein. However, transarterial lead wires have not been chronically implanted in clinical scenarios and further in vivo studies are required to assess the chronic effects of transarterial lead wires or preferably a full wireless approach may be considered [20]. The artery diameter changes during a cardiac cycle due to pulsation and blood pressure variation. However, the change in electrode-nerve distance was estimated to be less than a millimeter and connective tissue in a neurovascular bundle is expected to dampen small relative movements resulting from pulsations. While the activation threshold is not expected to generate noticeable differences, it should be considered when performing in-vivo studies as variability in the thresholds may be explained in part by the pulsations.

Feasibility of endovascular stimulation of peripheral nerves
With an anatomically realistic fascicular map, we predicted that the endovascular approach could efficiently recruit myelinated axons in the pudendal nerve from electrodes placed within the internal pudendal artery. All configurations tested achieved >60% predicted activation of the total myelinated axon population and some configurations led to nearly or above 90% activation at a maximal stimulation intensity of 10 mA. With an optimal interelectrode distance, endovascular disk electrodes were predicted to have comparable axon recruitment efficiency to conventional electrodes directly placed on the pudendal nerve, only increasing the mean activation threshold by roughly 0.35 mA. These results predicted that activation of the pudendal afferents for bladder control is feasible using the endovascular approach. However, the simulation results indicated that it may be infeasible to recruit unmyelinated fibers with endovascular stimulation within the parameters tested. Further in vivo studies with compound action potential recordings are required to verify whether unmyelinated fibers can be recruited via endovascular stimulation within a safe stimulation intensity range, and to investigate the functional implications of the possible lack of C-fiber activation.
We also modeled endovascular stimulation of the cervical vagus nerve with anatomically accurate tissue conductor models and a randomly generated fascicular arrangement. The vagus nerve projects from the brainstem and courses bilaterally in a neurovascular bundle with the internal jugular vein and the carotid artery, before forming an abundance of branches to innervate the viscera. Our model predicted relatively low axonal activation levels of the compound vagus nerve for all endovascular disk electrode arrays tested, ranging from 6.1% to 37.1% of all the myelinated axon population. However, the results also indicated that, when placed distal to the branching points, the endovascular electrodes can achieve a high selectivity in only activating distinct branches in proximity to the internal jugular vein.

Optimization of endovascular stent-electrode-array design
Our model showed that insulated electrodes improved axon recruitment, but a thin blood layer between the electrode and blood vessel surface circumvented this advantage by shunting current away from the neural tissue. This shunting effect may affect stimulation immediately after implantation and prior to incorporation into the blood vessels. However, after approximately 14 d the endovascular electrodes become incorporated into the vessel wall and will be separated from the surrounding blood by tissue [21,45].
A small inter-electrode distance of 2 mm led to a mean activation threshold of more than two times the minimum mean threshold. This was likely due to a significant amount of current being shunted away by the blood vessel wall tissue. This shunting effect with small inter-electrode distances has also been predicted by previous modeling studies [38,46]. As the distance between bipolar electrodes increases, tissue impedance on the direct path between the electrode pair also increases, which may lead to a broader current spread and more current through the neural target. As the inter-electrode distance increases beyond a certain threshold, the configuration starts to resemble a monopolar one, where current spread is not limited to the local region and the potential gradient in the neurovascular bundle is reduced, which conversely reduces stimulation efficiency. This could explain the initial rapid decay and subsequent slow rise in the mean activation threshold in relation to an increasing inter-electrode distance predicted by this model. This study predicted an optimal longitudinal distance of 8-10 mm between the active and return endovascular disk electrodes to minimize the activation threshold.
Surprisingly, varying electrode size did not have a significant impact on stimulation efficiency, contrary to the findings in previous modeling work on other neurostimulation techniques [39,47], but in accordance with the results of a previous modeling study that evaluated the feasibility of deep brain stimulation using endovascular ring electrodes [25]. This could be partially attributed to the small range of electrode sizes tested due to the spatial constraints within a blood vessel. Although large disk electrodes do not significantly alter the stimulation outcome, they provide benefits such as lowered surface charge density and lower susceptibility to stent migration. However, large electrodes may also hinder proper compression of the stent during its delivery through a catheter and will have a larger contact surface area with the inner wall of the blood vessel. Ease of manufacture, biocompatibility of the material and ease of stent delivery and deployment should all be considered when determining the dimensions of endovascular electrodes.
The bar electrode configuration 2 shared the same effective inter-electrode distance 10 mm as the vertical configuration 4. Despite increasing the effective electrode surface area by 4.67 times, configuration 4 only increased axon recruitment by less than 5%. This indicated that inter-electrode spacing and electrode orientation, but not electrode shape and surface area, have the largest impact on axon recruitment. The single disk electrode pair configuration with an effective inter-electrode distance of 8 mm yielded the highest axon recruitment. However, the single pair configuration also had the highest charge density on electrode surfaces at 339 µC cm −2 with current amplitude of 10 mA, well above the recommended limit of 30 µC cm −2 for stimulations involving brain tissue. Nevertheless, it must be noted that the highest current density within the neural tissue under the same condition was predicted to be only 7.65 µC cm −2 due to the distance between the nerve and the electrodes. Indeed, in clinical applications where electrodes are separated from the neural tissue, current densities significantly higher than the established limit can be used [48]. However, at this high focal charge density, electroporation of the blood vessel can become a risk. Future in vivo studies are required to assess the chronic effects on the blood vessel wall and neural tissue through histological examination and nerve conduction velocity evaluation. Other means of improving the charge capacity should be considered, such as iridium oxide electrodes or other novel coatings.
The vertical configurations with sufficient electrode spacing reduced electrode surface charge density without sacrificing efficiency in axon recruitment. However, these configurations are more susceptible to stent misplacement. Currently, there is no established technique to control stent orientation during endovascular stent deployment, and rotation of the electrodes from the optimal position facing the target nerve may significantly reduce stimulation efficiency. The more anatomically realistic vagus nerve model demonstrated the effect of electrode misalignment on axon recruitment. Due to the non-idealized alignment of the vertically arranged electrodes and the vagus nerve, configuration 4, the six horizontally positioned electrodes, outperformed all vertical configurations. The horizontal configurations are less prone to misplacement as the electrodes cover a larger circumference of the blood vessel. More electrodes may be added to cover the entire circumference; however, as the number of electrodes in effect increases, or with an increasing angular distance between electrodes of the same polarity, the performance of a horizontal configuration array may asymptote towards that of a pair of ring electrodes, with a wide current spread that reduces the focal current density at the target nerve. The significantly lower fiber recruitment of the ring electrode pair was consistent between the pudendal nerve and vagus nerve models. Increasing the number of electrodes conforming to the blood vessel circumference may also become a hindrance to stent delivery. The trade-offs of each type of configuration must be taken into consideration when arranging electrodes on a stent.

Limited fascicular selectivity with the endovascular approach
According to histological sections of the compound pudendal nerve and sectioned images of a human female pelvis [28,29], the internal pudendal artery is located ventromedial to the pudendal nerve within the pudendal canal at a distance of about 1.5 mm. Therefore, the fascicles that distally make up the DGN are at a shorter distance to the internal pudendal artery, whereas the fascicles that make up the perineal branches are further away from the artery. Previous histology study has shown that the relative position of the nerve and artery as well as the fascicular pattern remain relatively consistent between individuals [29]. A major constraint on the performance of endovascular neuromodulation is the course of the implanted blood vessel and its distance from the target nerve. It can thus be predicted from this relatively stable anatomy which fascicles and consequently distal branches of the pudendal nerve are accessible via endovascular stimulation.
The pudendal nerve gives rise to perineal and dorsal genital branches within the pudendal canal [29]. Depending on the stimulation frequency and the state of the bladder, stimulation of the pudendal afferents may produce an excitatory or inhibitory effect on bladder contraction [39,49]. This reflex persists after a spinal cord injury, indicating that it most likely only relies on spinal circuitry and not supraspinal control [50], rendering electrical stimulation of the pudendal nerve a promising technique to restore bladder functions in patients with spinal cord injuries.
Although electrical stimulation of the pudendal nerve and its afferent branches was able to evoke robust bladder contraction in animal models [50,51], little success has been seen in stimulation of the compound pudendal nerve in people with spinal cord injuries. This may be due to conflicting excitatory and inhibitory effects on the spinal-bladder circuits, resulting from co-activating functionally distinct pudendal nerve components.
This study predicted that, although constrained by the anatomy of the neurovascular bundle and the fascicular anatomy of the compound target nerve, endovascular stimulation can achieve some degree of selective activation of individual fascicles. The pudendal nerve model predicted that all electrode configurations consistently achieved notably higher fiber recruitment for the DGN branch of the compound pudendal nerve. Therefore, we predicted that endovascular stimulation from within the internal pudendal artery will likely activate DGN branch of the pudendal nerve, a common neurostimulation target in the management of bladder functions, given that the fascicular anatomy of the pudendal nerve is relatively stable between individuals. This passive selectivity may be beneficial in terms of therapeutical efficacy as it prevents possible conflictive stimulation outcomes of the perineal and DGN afferent branches. Our results showed that stimulation through the internal pudendal artery activated the motor branches before DGN fibers activation reached a substantial level. Activation of the pudendal efferent pathways leads to contraction of the external urethral sphincter, an undesired outcome during the voiding phase. Targeted stimulation may be needed for intravascular stimulation. One way to overcome this shortcoming may be to use a combination of proximal low frequency stimulation and distal kilohertz nerve block to activate bladder contractions through the pudeno-vesical reflex while blocking efferent activation [52].
The vagus nerve innervates a myriad of visceral tissues and organs through parasympathetic efferent pathways and conveys sensory information via afferent fibers from visceral organs [42]. Indiscriminate stimulation of all the fascicles of the cervical vagus nerve can lead to adverse side effects generated by off-target stimulation. Although the organ-specific functional organization of the vagus nerve is mostly unknown, evidence showed that some functional outcomes of vagus neuromodulation resulted from localized stimulation within the nerve [53]. This study also predicted that, by changing the relative positioning of the endovascular electrodes to the anatomical branching points of the cervical vagus nerve, selective activation of vagal fascicles and, therefore, potentially selective functional outcomes may be achievable. Placing the endovascular electrodes distal of the vagus nerve branching locations was predicted to impose high selectivity on distinct nerve branches depending on their distance from the internal jugular vein.

Conclusion
Endovascular PNIs provides a minimally invasive alternative neurostimulation modality that preserves the structural integrity of the neural target and mitigates the risk of fibrotic encapsulation. This study expanded the scope of application of the endovascular stent-electrode array to peripheral nerve stimulation and investigated the feasibility of this approach through computational modeling. The models predicted the feasibility of endovascular stimulation of the pudendal nerve and vagus nerve, and that the stimulation threshold is dependent on inter-electrode distance and electrode array orientation in relation to the target nerve. The arrangement of the electrodes on the endovascular stent can be flexibly optimized to lower the stimulation threshold. In comparison to endovascular ring electrodes, a small number of disk electrodes facing the target nerve generally lowered the activation threshold by improving current steering. Overall, these results suggest that the endovascular stent-electrode array is a promising alternative neurostimulation modality for peripheral nerves that would typically require invasive open surgery for electrode implantation and calls for further in vivo investigations.

Data availability statement
The data that support the findings of this study are available upon reasonable request from the authors.