Towards enhanced functionality of vagus neuroprostheses through in silico optimized stimulation

Ethical statementAll animal protocols and surgical procedures were reviewed and approved by the Animal Care and Use Committee of Feinstein Institutes for Medical Research and New York Medical College.Immunohistochemistry of nerve samplesAfter the humane euthanasia of two pigs (male Yucatan ~40 kg, ~1 year old, Premier BioSource) via Euthasol injection, a meticulous surgical dissection procedure was undertaken, exposing and isolating both the left and right cervical vagus nerves, from above the nodose ganglion down to the termination point in the thoracic vagus. To enable subsequent imaging and tracking of vagal branches, sutures were delicately placed using a 6-0 fine suture loop, affixing them to the epineurium of each branch near its emergence from the trunk. Detailed photographs were captured before and after the nerve extraction process, providing comprehensive documentation of the nerve trunk and its associated labels. Following this, the isolated nerve samples underwent meticulous fixation in 10% formalin for a duration of 24 h. For the immunohistochemistry (IHC) studies, both the left and right vagus nerves (n = 2) were meticulously sectioned at a thickness of 5 microns at four different levels: initially at the nodose ganglion, followed by the upper cervical level, mid cervical level, and lastly, the lower cervical level. From each of these segments, a total of 200 sections were obtained, with each section measuring 5 μm in thickness. Approximately 25 of these sections were then selected for the IHC study. These sections underwent staining for specific markers, including myelin basic protein (MBP), neurofilament (NF), choline acetyltransferase (ChAT), tyrosine hydroxylase (TH), and voltage-gated sodium channels 1.8 (Nav1.8). Additionally, a distinct set of adjacent sections was utilized for H&E staining to identify cellular structures such as endoneurium, perineurium, and epineurium. These staining procedures strictly adhered to standard IHC protocols, encompassing essential steps such as antigen retrieval, blocking, incubation with primary and secondary antibodies, and rigorous washing steps. The imaging of the stained sections was conducted utilizing advanced equipment, including the ZEISS LSM 900 confocal laser scanning microscope and the BZ-X800 all-in-one fluorescence microscope.Segmentation and characterization of single fibers from IHC imagesThe vagus nerve is segmented using a combination of standard computer vision and deep learning techniques from immunohistochemical images16,74. We employed a systematic approach to detect specific anatomical markers in the nerve: neurofilament, myelin, ChAT, and Nav1.8. For each marker, a brightness threshold was empirically selected to classify positive pixels. Blobs or connected components of these pixels were then identified, with noise below the said threshold discarded. For detailed fiber segmentation, the deep learning model Mask R–CNN was initially used to detect axons. It was then complemented with traditional methods to address overlapping neurofilament blobs that the deep learning model missed, achieving a fiber detection accuracy of 96%.After segmenting neurofilaments and assigning myelin to detected fibers, various features, such as fiber diameter, perimeter, and myelin thickness, were extracted. Fascicles were manually annotated to extract features at that level. To generate a comprehensive dataset, manual annotations of Nav1.8-stained images were made to identify the density of C-fibers within Remak bundles, enabling the creation of a linear equation to estimate fiber counts in other bundles. Image pre-processing, like adaptive histogram equalization and color-clipping, was vital to enhance contrast and eliminate noise. Finally, fibers were classified into four distinct categories based on morphology and function: myelinated and unmyelinated efferent and afferent.Creating a realistic 3D model of the cervical vagus nerveFour distinct histological slices (C1 to C4) of the explanted vagus nerve were extracted from the mid-cervical regions of both animals (labeled M1 and M2). C1 was located at 40 mm from the nodose ganglion and served as a reference for the other slices. Toward the thoracic region, these were located at distances of 10 mm (C2), 25 mm (C3) and 35 cm (C4) from C1, spanning a total length of 35 mm.We manually segmented the contours of the epi- and endoneuria of the histological slices using ImageJ with the NeuronJ plugin75,76, yielding the anchor point coordinates of two-dimensional splines. These were imported into Solidworks, respecting the correct distances between the individual slices. All splines were then simplified by applying Solidworks’ Simplify Spline function with a tolerance under 0.005 mm, until no further reduction in spline points could be achieved. This did not appreciably alter the shape of any contour, but significantly improved the rendering efficiency and later simplified the mesh. Next, the fascicular contours at the different cross sections were matched visually to determine the branching and merging paths of each fascicle, which were then reconstructed using the Lofted Boss/Base feature. Using the corresponding cross sections, the epineurium was also lofted. Lastly the feature Interference Detection was utilized to detect any intersections between fascicles as well as fascicles protruding outside of the epineurium, which were manually resolved by adding appropriate guide splines to the lofts.The helical cuff electrode used in the experimental setting was modeled after the manufacturer specification. The electrode is composed of 8 surface active sites radially distributed around the nerve. Each active site is squared and has an area of 0.25 mm2. They are laterally separated by 0.5 mm. The electrode substrate was wrapped around the epineurium using the Wrap function.Setting up a finite element model for vagus nerve stimulationIn order to solve the electrical fields that the active sites induce during stimulation, the 3D model was imported into COMSOL Multiphysics. Here a homogeneous saline cylinder with radius of 35 mm was added to emulate the interoperative environment, with the diameter chosen by a convergence study so that the ground-at-infinity condition was sufficiently approximated30,37,44,45. The perineurium was modeled by adding a thin layer contact impedance to each fascicle boundary, with the surface thickness being set to 3% of the average area over all four cross sections of the corresponding fascicle. This approach leaned on the common assumption that the endoneurium thickness equals about 3% of the fascicle area37,44,46,48,77. Epineurium, perineurium, electrode substrate and saline were assumed to be isotropic, whereas the endoneurium was modeled with anisotropic conductivity, due to the presence of fibers within. The precise conductivity values are based on previous modeling studies44, and listed in Supplementary Table 2. To ensure that the anisotropy is correctly oriented with respect to the fascicular cross sections, and not the global coordinate system, we used COMSOL’s Curvilinear Coordinates feature. In particular, the anisotropic vector field was aligned to the streamlines determined by a diffusion study, in which the fascicle hulls were defined as walls and their flat ends as inlets and outlets. Both nerve models were meshed resulting in about 15 to 20 million tetrahedral elements. For low frequency stimulation, it is adequate to exploit the quasi-static approximation of Maxwell’s equations in order to calculate the electric fields \({{\rm{V}}}{{\rm{ext}}}\)78 (Eq. (1)), where \({{\rm{\sigma }}}\) is the conductivity:$$\nabla \cdot {{\rm{\sigma }}}\nabla {{\bf{V}}}_{\!{{\bf{ext}}}}=0$$
(1)
As this equation is linear, it is possible to determine the electric fields resulting from an arbitrary multipolar stimulation by scaling and superimposing electric fields which were independently determined for all active sites under consideration of unit charge injection44 from surface current sources corresponding to each active site.Populating the model with histologically accurate fiber dataFor both nerve models, precise fiber data was obtained from the immunohistological analysis of cross section C2 in a fascicle-wise fashion. Fiber types and modalities were differentiated by correlating the responses of NF (neurofilament), MBP (myelin basic protein), ChAT (choline-acetyl-transferase) and Nav1.8 stains16. We extracted precise locations and diameters of all myelinated fibers. Due to imaging resolution limits, single unmyelinated fibers could not be labeled. Instead, we identified centroids and areas of clusters of unmyelinated fibers. A supplementary study with higher resolution was conducted, yielding the empirical diameter distributions for both afferent and efferent unmyelinated fibers as well as a linear correlation between the cluster area and count of contained fibers16. Using these results, the fiber count for each cluster was estimated and the diameters were assigned by inverse sampling gamma distributions (Eq. (2)) which had been fit to the empirical diameter distributions (Table 1). The precise locations of unmyelinated fibers were generated by uniformly sampling coordinates within a circle centered on the centroid of the cluster, with area equivalent to that of the cluster:$$p\left(x\right)=\frac{1}{{b}^{a}\Gamma \left(a\right)}{x}^{a-1}{e}^{-\frac{x}{b}}$$
(2)
Table 1 Parameters of the gamma distributions fitted to empirical diameter distributions of unmyelinated afferents and efferentsThe fibers were now integrated into the model by performing fascicle-wise matching between centroids of the modeled fascicles (of cross section C2) and the centroids of all corresponding fibers. Fibers placed outside of the endoneuria were subsequently removed.To reduce the computational costs we subsampled the simulated fibers. A convergence study was conducted to obtain the minimum number of fibers necessary to accurately reproduce the recruitment curves of the full data set. Using the criterion that the mean absolute deviation of recruitment curves should be lower than 5% for each fiber type, we determined that it is suitable to simulate at least 10’000 fibers per nerve (Supplementary Fig. 1). During subsampling, care was taken to ensure that the proportion of the individual fiber types as well as their modalities across the whole data set does not change substantially. The exception was to always sample all Aα fibers, due to their low occurrence, and to remove all myelinated fibers with a diameter smaller than 1 µm or larger than 18 µm due to the expected inaccuracy of the extrapolated MRG model for these values.The fibers were propagated along the nerve based on their placement in the initial cross section as well as the curvilinear coordinates extracted from COMSOL along each fascicle. This further ensured the precise alignment of each axon to the anisotropic vector fields of the endoneurium. Moreover, it ensured smooth fiber trajectories along fascicular branching. Table 2 lists the number of fibers placed and successfully propagated along both nerve models in the context of this paper.Table 2 Population sizes of the different fiber types for both nerve modelsExtracting extracellular potentials along fibersOnce the precise fiber coordinates in 3D space are known, it is possible to extract the extracellular potentials along the fibers resulting from unit charge injection by an active site. For this, the set up COMSOL model is solved for each active site individually by applying a reference current density of 1 A/m2. To obtain the extracellular potential corresponding to a 1 μA charge injection, this solution is then interpolated at discrete points along the fibers and scaled by the electrode area. The potential was extracted at the center of each section of each fiber.Modeling the membrane dynamics of myelinated and unmyelinated fibersTo simulate the response of fibers to extracellular stimulation, or, more precisely, the time-course of their membrane potential, it is necessary to model the ion channels as well as passive components. The dynamics of ion channels are usually formulated as modified nonlinear Hodgkin–Huxley equations. It must further be considered whether a given section is myelinated or not, as the myelin sheath acts as a local insulator. The equivalent electrical circuits for sections of myelinated and unmyelinated fibers are illustrated in Supplementary Fig. 8.Based on a comprehensive review33, where several models of unmyelinated fibers were compared, the Tigerholm (TH) model was selected to model unmyelinated fibers36. It was deemed to reproduce data from single-fiber recordings most accurately, albeit the experimental recovery cycle could not be captured. Furthermore, the TH model has been incorporated into modeling pipelines similar to the one presented in this work26. The TH model was implemented using the previously published files33. However, a discrepancy between the corresponding conductance values of the pump and \({{{\rm{K}}}}_{{{\rm{Na}}}}\) channels and those listed in the original publication36 was observed. The values of the latter were used (see Supplementary Table 3), which we confirmed to accurately reproduce the action potential current traces36.For myelinated fibers the frequently used McIntyre-Richardson-Grill (MRG) model was chosen17,26,37,38,44,45,47,50,58,79,80. It compartmentalizes the region between two nodes of Ranvier into ten distinct sections, the parameters of which are diameter dependent. As the original parameters were presented for only a small set of discrete fiber diameters, it is necessary to inter- and extrapolate their values such that the model can be applied to fibers of arbitrary diameters. In this context polynomial interpolation has previously been performed26, which was also adopted in this project although with slightly deviating results (see Supplementary Fig. 9 and Supplementary Table 4). Furthermore, a rounding imprecision of the published leakage conductance caused a small transient membrane hyperpolarization (\({V}_{{{\rm{init}}}}-{E}_{{{\rm{leak}}}}\)) to occur when the fiber was completely at rest. This was resolved by recalculating the leakage conductance \({g}_{{{\rm{leak}}}}\) at every node such that the leakage current exactly balances all ionic currents \({i}_{{{\rm{Na}}}}\), \({i}_{{{\rm{Nap}}}}\), and \({i}_{{{\rm{K}}}}\) at rest (Eq. (3)), an idea which was also implemented for the previously published TH model implementation33. This causes the leakage conductance to only deviate slightly from its original value:$${g}_{{{\rm{leak}}}}=-\frac{{i}_{{{\rm{Na}}}}+{i}_{{{\rm{Nap}}}}+{i}_{{{\rm{K}}}}}{{V}_{{{\rm{init}}}}-{E}_{{{\rm{leak}}}}}$$
(3)
All NEURON functionality was implemented using the NEURON-Python interface. For both the TH and MRG models, the original HOC scripts were ported to Python while the model description (NMODL) files were not altered.Determining the recruitment threshold of a fiberTo execute a simulation in NEURON the stimulation paradigm must be encoded in a pulse matrix of the form which contains the injected current of each active site in μA along with the time stamp at which a change in stimulation amplitude of any active site occurs. After a simulation is finished the current process is cleared of all remaining NEURON objects by iteratively deleting all detectable distributed membrane mechanisms, point processes as well as sections. This does however not remove the extracellular mechanism, which fundamentally alters the structure of the underlying Jacobian matrix to represent the additional layers in the cable model. In some cases, this might affect a subsequent simulation, which is circumvented by manually resetting the extracellular potential of each node to 0.The method outlined above is utilized as part of a bisection algorithm to determine the recruitment threshold of a fiber, which indicates by what factor the charge injection amplitudes in the pulse matrix must be scaled to elicit an action potential. Our pipeline natively supports execution of this algorithm on a High Performance Cluster (HPC), including automated bidirectional file transfer.Dynamic discretization of unmyelinated fibersWe spatially discretized unmyelinated fibers using sections of varying length, to optimally trade-off the number of sections and simulation accuracy. In this regard, we developed an algorithm (see Supplementary Listing 1) which assigns section lengths based on local estimates which evaluate how relevant a fiber region is for accurately predicting the physiological response of the whole fiber. This idea is motivated by the fundamental result that there is an approximately linear correlation between the subthreshold membrane (de)polarization \({V}_{{\mbox{mem}}}\) at node i and the value of the second spatial derivative of the extracellular potential \({V}_{{\mbox{ext}}}\) at node i along the axon (known as activation function, Eq. 4)55, where \({{{\rm{\rho }}}}_{a}\) is the axoplasm resistivity and \({c}_{m}\) the membrane capacitance:$$\frac{\partial {V}_{{\mbox{mem}},{\mbox{i}}}}{\partial t} \, \, \approx \frac{d}{4{{{\rm{\rho }}}}_{a}{c}_{m}}\frac{{\partial }^{2}{V}_{{\mbox{ext}},i}}{\partial {z}^{2}}$$
(4)
A prerequisite for the dynamic discretization algorithm is that the extracellular potential of each fiber has been extracted from COMSOL by interpolating the solution of the electric field at fine resolution. For our nerve models with a total length of 35 mm a resolution of 10 μm was selected, which yielded about 3500 nodes per extracellular potential per fiber. The dynamic discretization algorithm then proceeds by determining an estimate of the second spatial derivative of the extracellular potential along each fiber. Here, the extracellular potential is first padded by linear extrapolation at both ends of the fiber by 5% of its total length. The 2nd order central differences are now evaluated, which yields a noisy estimate of the second spatial derivative. Smoothing is subsequently performed using a 2nd order forward-backward Butterworth low pass second-order sections (SOS) filter. The cutoff frequency of this filter was set to 0.001. After smoothing the padding is removed. The positive and negative parts of the estimated second derivative are then normalized independently, such that all values fall within the range [0,1]. This renders the estimate invariant to the scale as well as sign of the applied extracellular potential. The fiber is then discretized so that each section length is proportional to the local second derivative within a predefined minimum and maximum section length. Lastly, the extracellular potential is resampled at the newly determined fiber nodes.Pseudocode of the proposed dynamic discretization algorithm for a fixed extracellular potential is given in Supplementary Listing 1. Note that generalization to multipolar stimulation with arbitrary stimulation paradigms and spatial distribution of electrodes is straightforward.Longitudinal truncation of fibersTo determine the truncation threshold for the normalized extracellular potential, based on which the truncation locations are defined, a study is conducted a priori on a sample set of fibers. This study is performed independently for both myelinated and unmyelinated fibers, and the result is then applied to all fibers in the nerve model.For the study, fibers with discretely assigned diameter values are randomly placed across the nerve. The initial recruitment threshold of each fiber is then evaluated using the full-length fiber. Next, a bisection algorithm along the normalized extracellular potential is executed to determine an appropriate value for the truncation threshold, whereas the fibers are truncated to the left and right from the outermost nodes corresponding to this value of the truncation threshold. If the recruitment threshold determined for the truncated fiber deviates by less than 0.5% from the full-length recruitment threshold, the truncation threshold is increased to further truncate the fiber; otherwise it is decreased, elongating the fiber again. The termination condition is that, for two consecutive iterations, the truncation threshold determined by the bisection algorithm corresponds to the same nodes, for which the recruitment threshold deviation is furthermore less than 0.5%. The determined truncation thresholds are grouped by fiber diameters. For both myelinated and unmyelinated fibers, a spline is then fit through the 0.5th percentile of all diameters to obtain a conservative diameter-dependent estimate of the threshold ratio. Pseudocode of the proposed longitudinal truncation algorithm for a fixed extracellular potential is given in Supplementary Listing 2.Increasing the integration time stepThe sweep over different time step values was conducted by evaluating all fibers in M1 under application of a cathodic pulse to active site 1, with a total simulation time of 3 ms and a threshold precision of 0.01 for MRG fibers and 0.1 for TH fibers.Studies for the improvement of computational efficiencyThe results in Fig. 2 and Supplementary Fig. 3 were obtained by stimulating M1 using the following parameters for the baseline model and presented methods (if applicable):Dynamic discretizationThe minimum and maximum section lengths were set to 20 μm and 200 μm, respectively. We discretized the fixed-length baseline models with sections of 20 μm length (yielding ~1750 sections for ~35 mm long fibers).Longitudinal truncationA study with a threshold accuracy of 0.01 for MRG fibers and 0.1 for TH fibers was conducted on an artificial data set. Fiber diameters were sampled from the discrete sets {1, 1.5, …, 18} μm (myelinated fibers) and {0.3, 0.8, …, 4.8} μm (unmyelinated fibers). Each diameter was assigned 50 times, but due to imperfect fiber propagation only 1694 myelinated and 478 unmyelinated fibers were eventually placed in the model.Time step increaseWe increased the time step from 5 μs in the baseline model to 13 μs in the modified model.We set the total simulation time to 3 ms and used a recruitment threshold accuracy of 0.1 µA. In Fig. 2 the nerve was stimulated using a monopolar cathodic pulse with 500 μs pulse width, whereas in Supplementary Fig. 3 the nerve was stimulated using a monopolar anodic pulse with 500 μs pulse width. Relative deviations of recruitment thresholds due to incorporating the three presented methods were evaluated as in Eq. (5):$${\delta }_{\lambda,{\mathrm{mod}}}=\frac{{\lambda }_{{\mathrm{mod}}}-{\lambda }_{{\mbox{orig}}}}{{\lambda }_{{\mbox{orig}}}}$$
(5)
where \({{{\rm{\delta }}}}_{{{\rm{\lambda }}},{\mathrm{mod}}}\) corresponds to the relative deviation, \({{{\rm{\lambda }}}}_{{\mathrm{mod}}}\) to the recruitment threshold obtained with the optimized method, and \({{{\rm{\lambda }}}}_{{\mbox{orig}}}\) to the recruitment threshold obtained with original method. Summary of the improvement in computational time is given in Supplementary Tables 1 and 5.Comparison of models with realistic and linearly extruded fascicle branchingModels with linearly extruded fascicles were created using the cross sections in which fibers are initially placed for both nerve models. This ensures that the same fibers are present in corresponding realistic and extruded models. The electrode was placed on three different locations, at a longitudinal distance of 2.5 mm, 7.5 mm and 12.5 mm from the extrusion cross section in direction of the thoracic region. For each fiber in each model, recruitment thresholds were estimated by independently applying a cathodic pulse to each active site. The relative deviation of the recruitment thresholds of corresponding extruded and realistic models (\({{{\rm{\delta }}}}_{{{\rm{\lambda }}},{{\rm{extr}}}}\)) was evaluated as in Eq. (6):$${{{\rm{\delta }}}}_{{\lambda },{{\rm{extr}}}}=\frac{{{\lambda }}_{{\mbox{extr}}}-{{\lambda }}_{{\mbox{real}}}}{{{\lambda }}_{{\mbox{real}}}}$$
(6)
where \({{{\rm{\delta }}}}_{{{\rm{\lambda }}},{{\rm{extr}}}}\) corresponds to the relative deviation, \({{{\rm{\lambda }}}}_{{\mbox{extr}}}\) to the recruitment threshold obtained with the extruded model, and \({{{\rm{\lambda }}}}_{{\mbox{real}}}\) to the recruitment threshold obtained with realistic model. The means of the distributions were tested by a one-sample, two-tailed t-test.Precise fiber placement can improve selectivity of synthesized stimulation policiesTo create models with uniformly distributed fibers, the locations of all fibers placed in the histologically accurate models were reassigned using Monte-Carlo rejection sampling. The magnitude of clustering of each fiber group was evaluated using silhouette values, which quantify clustering as the difference between the cohesion within the cluster and the separation with other clusters, and ranges from −1 and +156. Curves of the spatial selectivity indices \({\eta }_{{{\rm{group}}}}\) were evaluated using Eq. (7):$${\eta }_{{{\rm{group}}}}={\mu }_{{{\rm{group}}}}-{\mu }_{{\neg} {{\rm{group}}}}$$
(7)
as a cathodic pulse stimulation policy was scaled, where \({\mu }_{{{\rm{group}}}}\) is the recruited fraction of fibers belonging to the target group and \({\mu }_{\neg {{\rm{group}}}}\) is the recruited fraction of non-target fibers. Such metric is commonly used49,67 and is easily interpretable. As other spatial selectivity indexes (as in Eq. (8)), it ranges from −1 (recruitment of all non-target fibers, no recruitment of target fibers), to 1 (recruitment of all target fibers, no recruitment of non-target fibers). For each group, we selected the best active site as the one with spatial selectivity curve with the largest mean value. The presented recruitment curves were also obtained for the same active sites. The functional consequences of fiber clustering were studied by comparing the models with histologically accurate fiber location with four models where the same fiber population was spatially shuffled across the nerve (Fig. 3b and Supplementary Fig. 5).Physiological experiments on new subjectsAfter sedation with Telazol (2–4 mg/kg), three pigs (male Yucatan ~40 kg, ~1 year old, Premier BioSource) were placed on a table in a supine position with body temperature maintained between 38 °C and 39 °C using a heated blanket. Anesthesia was induced with Propofol (4–6 mg/kg, i.v.), and maintained with isoflurane (1.5-3%, ventilation). Then, a 4–5 cm long incision was cut on the cervical region of the neck. After dissecting the underlying tissues, the vagus nerve was identified and isolated. One or two helical cuff electrodes (CorTec, Germany) were placed on the rostral (around 2 cm to nodose) and/or caudal sites with a distance of 4 to 5 cm. The helical cuff included 8 contacts (0.5 mm by 0.5 mm) evenly distributed (1 mm between two contacts) around the circumference and 2 ring-shaped (0.5 mm by 7.4 mm) return contacts. EMG electrodes, which were teflon-insulated stainless-steel wires with de-insulation of 1 mm at the tip, were inserted in the thyroarytenoid (TA) of laryngeal muscle with a needle after blunt dissection of the hyoid bone. All procedures were performed using sterile techniques and approved by the animal care and use committee of Feinstein Institutes for Medical Research and New York Medical College.ECG signals were recorded and amplified (FE238, ADI) by using a 3-lead patch electrode configuration. Blood pressure was recorded from the aorta with catheter (SPR-751, Millar Inc) after amplification (FE228, ADI). All physiological signals were digitized and acquired at 1 kHz (PowerLab 16/35, ADI). Neural signal and EMGs were sent to the data acquisition system through a headstage (RHS 32, Intan Tech). It was filtered (60-Hz notch), amplified, and digitized at a sampling frequency of 30 kHz with an acquisition system (Intan 128 recording/stimulating, Intan Tech).Biphasic stimulus pulses were delivered using an isolated constant current stimulator (STG4008, Multichannel Systems) through each of the 8 contacts of the cuff electrode. Physiological threshold (PT) was detected from stimulus trains of 10 s duration (30 Hz, pulse width 200 μs). There were at least 30 s long waits between successive trains to ensure that physiological measurements had reached a steady state before a new train was delivered.After removal of the direct current (DC) component, the evoked compound action potentials (eCAP) were averaged from each stimulus and classified into different fibers groups based on the conduction speed. The first peak (0.8 ms) shows A-fiber response (conduction speed, 5–120 m/s), and the late peak (~1.5 ms, ~3 ms) shows A-delta/B-fibers (conduction speed, 5–15 m/s). Recruitment maps of different fiber responses (peak-to-peak amplitude) in the recording cuff (8 contacts) were derived as second cuff was stimulated (8 contacts).Model personalization to new experimental subjects and validationWe developed a method to adapt these highly detailed computational models to new experimental subjects for which no histology can be obtained (in vivo). We first simulated the stimulation paradigm used in the experiment (monopolar, biphasic, charge-balanced squared pulse of pulse width 200 μs). We applied a linear regression model with no intercept with input the model recruitment curves of Aα efferents per active site and output the experimental recruitment curves for laryngeal EMG per active site. We then found the scaling factor that applied to the fiber thresholds predicted by the model maximized the R2 of the linear regression model. Since the electrode is radially symmetric, we computed this maximized R2 for all possible shifts of active sites around the nerve circumference (n = 7), and chose the shift of active sites that, together with the scaling factor, predicted best the laryngeal EMG (Fig. 4a, b). We computed the personalization steps for each model (M1 and M2) and experimental subject (S1, S2, and S3) for a total of 6 personalized pairs of models-experimental subject. For each pair, we then compared the final R2 with the R2 of the original non-personalized model with a Wilcoxon signed rank test. We then validated the personalized model by comparing the model recruitment curves for Aα and Aβ efferents with the corresponding fast compound action potential recorded experimentally16, per active site, normalized to [0,1] between minimum and maximum CAP recorded per-subject. We then computed Pearson’s correlation coefficient between recruitment curves (Fig. 4c ii and Supplementary Fig. 6), and we compared the thresholds to obtain 10% to 70% recruitment in model and experiment (in steps of 10%), using a Wilcoxon signed rank test (Fig. 4c iii and Supplementary Fig. 6). The higher boundary of 70% was chosen since Pig 3 did not reach higher recruitment values in the experimental range. Finally, we compared the ranking of active sites between model and experiment by computing the Spearman’s correlation coefficient between the thresholds predicted by the model and by the experiment per active site for each recruitment level (10% to 70%). Finally, to exclude that the result was a result of chance, we performed the same computation of correlation randomly shuffling the active sites in the model 100 times. We compared the two distributions of correlation coefficients with a Wilcoxon rank sum test (Fig. 4c iv and Supplementary Fig. 6).Prediction of hearth rate variation and model-guided reduction of off-target effectsWe predicted the experimentally measured heart rate reduction with the fiber activation estimated by personalized models through multiple linear regression. The model inputs the number of recruited BEff in each fascicle to the predictor (for a certain active site and injected current level), which outputs the predicted heart rate reduction. No intercept was included since there must be no heart rate variation when no fiber is activated. The weights were constrained to be positive to reflect the hypothesis that BEff activation causes heart rate reduction (with positive absolute values), and not increase. The fascicles containing a low number of BEff were excluded using k-means clustering with k = 2. Further, we excluded the fascicles where no BEff was recruited within the current range of the experiments (Fig. 5a, b and Supplementary Fig. 7). We then computed the RMSE between model prediction and experiment (Fig. 5c and Supplementary Fig. 7). From the multiple linear regression we extracted the combination of fascicles responsible for experimental HR changes (weight > 0). Using selected fascicles, we identified for each personalized model the active site with the lowest threshold to obtain 25% recruitment of BEff across the target fascicles. In total, 25% was defined by selecting the experimental charge for targeted HR change (2500 μA16) and by computing and averaging the BEff recruitment for each personalized model. We then computed for that active site the recruitment curves for Aα across the whole nerve (correlated with L-EMG activity) and for BEff in the fascicles which were assigned a non-zero weight by the linear regression. Then, we performed an additional study, simulating tripolar stimulation with a charge-balanced asymmetric waveform. The active site identified in the previous step was set as pseudo-cathode, with a first cathodic pulse with 200 μs width and a second balancing anodic pulse with 20 μs width and one tenth of the amplitude. The active sites on either side of the pseudo-cathode were configured as pseudo-anodes, with opposite sign and each half amplitude than the pseudo-cathode (Fig. 5a ii, d). We computed the recruitment curves for Aα and target BEff as for the monopolar case (Fig. 5d). We then identified the level of recruitment of Aα that corresponded with the 25% recruitment of target BEff for each personalized model and compared the distributions between the monopolar (experimental) and the tripolar case with a Wilcoxon signed rank test (Fig. 5f). Finally, we compared the current thresholds required to obtain the 25% recruitment of target BEff for both stimulation paradigms and compared them with Wilcoxon signed rank test (Fig. 5g).Evaluation of electrode performance due to surgical placement uncertaintyFor each nerve, eighteen models with different electrode placements were created. The three locations were chosen such that the center of gravity of all active sites coincides with a point at a longitudinal distance of −5 mm, 0 mm, or +5 mm from the center of the nerve model. At each location the electrode was rotated by the angles 0°, 60°, 120°, 180°, 240° and 300°. Using a cathodic pulse the recruitment thresholds of all fibers were generated for all active sites of all models. The best spatial selectivity indices of all fascicles were then calculated using Eq. (8):$${\eta }_{i}={\max }_{\left\{A{{\rm{S}}}1,{{\rm{AS}}}2,\,\ldots \right\}}\left\{{\mu }_{i}-\frac{1}{m-1} {\sum}_{j=1,\, j\ne i}^{m}{\mu }_{j}\right\}$$
(8)
where \({\mu }_{i}\) is the fraction of recruited fibers in fascicle i and m is the total number of fascicles44. All fascicles with \(\eta \, < \, 0.7\) were omitted from further analysis. To investigate the effect of electrode rotation, the selectivity indices of all six electrode angles were grouped for each of the three longitudinal locations; to analyze the effect of electrode translation, the selectivity indices of all three longitudinal locations were grouped for the six electrode angles; and to analyze the combined effect of electrode rotation and translation, the selectivity indices of all eighteen models were grouped. For each fascicle, the range spanned by \(\eta\) over all values in each group was determined. These ranges of \(\eta\) were plotted in Fig. 6d–f.Online platformThe platform uses Three.js for in-browser 3D rendering with WebGL. Since it is written in client-side JavaScript, it is fully static and can therefore be hosted on any web server. It requires the offline pipeline to export the components of the model in formats which can be read natively by JavaScript. Geometries of nerve and electrode are saved in STL format, which can be natively rendered by Three.js; numeric data (i.e., fiber trajectories and thresholds) are exported into a plain binary format together with metadata (number of dimensions and size of each dimension), allowing the reading and writing of multidimensional arrays; and model specifications are saved as JSON files.Technical specificationsThe current iteration of the pipeline employs MATLAB 2022a, Python 3.9.12, COMSOL Multiphysics 5.6, and NEURON 8.2.1. Histological slices were segmented with ImageJ 1.52p and NeuronJ 1.4.3. The 3D models were created using Solidworks 2020 SP04. NEURON simulations were performed on the Euler cluster, operated by the High Performance Computing group at ETH Zürich, employing 48 cores of the AMD EPYC 7H12 CPU with 2.6 GHz nominal and 3.3 GHz peak clock speed as well as 98304 MB RAM with 3200 MHz clock speed.Statistical analysisAll statistical analyses were performed on MATLAB R2020a (The MathWorks, Natick, USA). All data were reported as mean values ± SD (unless elsewise indicated). The normality of data distributions was evaluated with a one-sample Kolmogorov–Smirnov test. Significance levels were 0.05 unless differently specified. Asterisks in figures indicate the level of statistical significance, *p ≤ 0.05, **p ≤ 0.01, ***p ≤ 0.001.Reporting summaryFurther information on research design is available in the Nature Portfolio Reporting Summary linked to this article.