Tutorial: a computational framework for the design and optimization of peripheral neural interfaces

Peripheral neural interfaces have been successfully used in the recent past to restore sensory-motor functions in disabled subjects and for the neuromodulation of the autonomic nervous system. The optimization of these neural interfaces is crucial for ethical, clinical and economic reasons. In particular, hybrid models (HMs) constitute an effective framework to simulate direct nerve stimulation and optimize virtually every aspect of implantable electrode design: the type of electrode (for example, intrafascicular versus extrafascicular), their insertion position and the used stimulation routines. They are based on the combined use of finite element methods (to calculate the voltage distribution inside the nerve due to the electrical stimulation) and computational frameworks such as NEURON (https://neuron.yale.edu/neuron/) to determine the effects of the electric field generated on the neural structures. They have already provided useful results for different applications, but the overall usability of this powerful approach is still limited by the intrinsic complexity of the procedure. Here, we illustrate a general, modular and expandable framework for the application of HMs to peripheral neural interfaces, in which the correct degree of approximation required to answer different kinds of research questions can be readily determined and implemented. The HM workflow is divided into the following tasks: identify and characterize the fiber subpopulations inside the fascicles of a given nerve section, determine different degrees of approximation for fascicular geometries, locate the fibers inside these geometries and parametrize electrode geometries and the geometry of the nerve–electrode interface. These tasks are examined in turn, and solutions to the most relevant issues regarding their implementation are described. Finally, some examples related to the simulation of common peripheral neural interfaces are provided. Neural interfaces with implantable electrodes are used to modulate and restore function to the peripheral nervous system. Hybrid modeling described in this protocol is used to optimize each aspect of the implantable electrode design and operation.


Hybrid modeling
Modeling the effects of neural stimulation plays a pivotal role in the design of neuroprosthetic devices. It can help reduce the need for preliminary animal experiments and the related ethical and economical costs, thus accelerating the translational process. The extensive use of models could also facilitate accelerating the optimization of neuroprosthetic applications by pre-selecting the most efficient stimulation patterns and thus reducing the number of tests involving patients.
In computational neuroengineering, the use of HMs is a modeling technique implemented in two steps: (i) determination of the electric potential induced by a stimulating electrode in a biological structure (volume conduction problem) and (ii) prediction of its consequences for single-neuron responses (neural response determination).
An HM is normally implemented through finite element modeling (FEM) for the volume conduction problem and detailed neuron models for neural response determination. FEM is a numerical method for solving partial differential equations; in an HM, FEM is used to solve the Poisson equation defined over complex geometry domains with given boundary conditions. For very simple approximations 1 (a point source in an infinite homogeneous, isotropic medium; in infinite homogeneous anisotropic medium; at the interface between two different homogeneous, isotropic media), FEM can be substituted by analytic methods, which will not be covered here. FEM is useful for more complex geometries, where solving the volume conduction problem by integrating the Poisson equation with the needed boundary conditions would be cumbersome or intractable.
The use of HMs relies heavily on the seminal papers by McNeal 2 and Rattay 3 , who provided a first view on the modeling of extracellular electrical stimulation of nerve fibers, and it was defined in the present form by  . The implicit assumptions are that the membrane potential varies only longitudinally with respect to the fiber course, which leads to neglecting transverse excitation modes, and that the fiber excitation does not affect the electric potential in the surrounding area. Additionally, it is generally assumed that quasi-staticity hypotheses hold 7 , namely that all the media are purely resistive and thus that capacitive/inductive reactive phenomena can be neglected.
HMs have been extensively used in the context of spinal cord stimulation [8][9][10][11] , deep brain stimulation [12][13][14][15] and peripheral nerve stimulation (PNS). This tutorial specifically addresses peripheral neural interfaces that directly contact the nerve to be stimulated: the general framework can be easily adapted to different applications (for example, transcutaneous stimulation 16 ).
Hybrid modeling can be used to assess the robustness of a given experimental setting with respect to variations in a subset of partially controllable parameters. For example, surgical insertion parameters for intraneural electrodes (defining the position of the active sites in the nerve section during stimulation) can strongly affect the results of experiments, but they are only partially controllable (they have very large uncertainty). Simulation can help establish bounds on these quantities of interest: if during the experimental procedure care is taken to stay inside these bounds (no matter what the exact values of these quantities is), then the stimulation result should not be affected significantly 24 .
At the same time, we can establish ranges of variation for noncontrollable parameters in our model by comparing the outcome of an experiment with different parametrizations of the given phenomenon. For example, in ref. 21 , the description on encapsulation tissue was investigated by formulating hypotheses on its conductance, and comparing the results with what is expected from stimulation.
Notice that, in the above, the key distinction is between partially controllable parameters (where we assume that the model is a perfect representation of reality and use it to guide experiments) and noncontrollable parameters (where we calibrate our model on experiments and thus modeling is guided by experimental activity). These parameters are listed in Box 1.
These perturbative/robustness analyses must be repeated for a reasonable range of analogous nerve samples, so that we can assure that our results are independent from the single subject we simulated. This requires possessing a statistical description of such nerve samples, so that many repetitions of the simulation routine for slightly different but reasonable samples can be performed.
The presented method has been developed and used by our group in a series of papers in the last ten years. Our HM framework was defined in refs. 18,19,24 and validated in rats 25 . More recent works expanded the framework to modeling of invasive peripheral stimulation for human upper and lower limb sensory feedback 21,22,27 . Moreover, our HM has also been used to optimize intraneural stimulation for vision restoration 28 and hand function restoration (Badi-Dubois, under preparation). Similar HMs have been developed by other teams 20,29 .

Comparison with alternative approaches
The goal of this tutorial is to provide a set of guidelines concerning the development of an HM for PNS based on the use of FEM and detailed neuron models, with attention to all of the issues previously mentioned, particularly parameter selection and model validation.

Sim4Life
Sim4Life combines multiphysics (electromagnetism, acoustics, heat flow and transmission) FEM and neural simulation optimized software with human and animal validated anatomical models and is thus well suited for the solution of very computationally intensive problems. Furthermore, Sim4Life offers customer service which helps the researcher to develop extremely detailed models to simulate specific experiments. Then, when a very specific set up has already been established, Sim4Life can be used as a tool to obtain very accurate simulations for that given set up. Our aim here is somehow complementary to Sim4Life's, as we seek an easy-to-use framework to be able to guide experiment design through the investigation of general principles on hierarchies of approximated models. Our computational framework is explained here in detail, with attention to the construction of modularity and complete customizability. Each block of our workflow is completely independent and provides useful insights into the mechanisms involved in the design and development of a neuroengineering application.
The main drawbacks of Sim4Life are its being a commercial software and the need for dedicated training to manage the complexity of its detailed workflows and methods. Our framework has been developed using MATLAB and COMSOL for the volume conduction problems, and Python and NEURON for the neural computational problem: all of these software packages and languages are either open source (Python and NEURON) or multipurpose (MATLAB and COMSOL), and almost every university owns the necessary licenses. For this reason, almost no dedicated training need be provided. In conclusion, Sim4Life remains an invaluable tool for the solution of very specific and detailed HM problems, whereas our method aims at creating a general and easily usable framework.

PyPNS
PyPNS consists of an efficient neuronal dynamic simulator for HMs. A precomputed electric field is used to extracellularly stimulate compartmental models for myelinated (McIntyre-Richardson-Grill (MRG) model) and unmyelinated (Sundt) fibers, exploiting the quasi-staticity hypothesis to calculate a single static value for volume conduction. Our approach expands this work by adding a completely customizable and modular framework to obtain the extracellular field associated with each simulated fiber (from the construction of representative geometries) and to interpret the results of stimulation experiments in terms of quantities of common use in neuroprosthetics (such as recruitment and selectivity).

Organization
Our review is divided into four main parts, which refer to peripheral nerve (Section 1), nerve fiber (Section 2) and neural interface (Section 3) modeling and to analysis of stimulation results (Section 4). Each subsection contains both a rationale for the introduction of the described entities into our conceptual framework and some implementation notes, which pertain to computational aspects.
The complete workflow can be followed through the schemes in Figs. 1-4. They describe: (1) How to numerically define the geometry of the model.
Histological data are used to draw simplified images of the geometry of the components of the nerve fiber; the dimensions of the chosen electrode are defined; the geometry of the electrode-nerve interface is defined; and changes to the geometry of the nerve fiber components as a result of inserting the electrode are considered. (2) Assigning properties to the different components. This involves specifying the material properties of each component of the nerve and electrode; setting the boundary conditions for simulating zero potential at infinity and active site current injection; and defining a mesh able to resolve the spatial electric potential map caused by the electrode forcing action. (3) Constructing nerve fiber models. Determine the branching pattern for the nerves, and decide the best sites in this branching system for positioning the electrodes; plot the distribution of fiber diameters across the sample; draw the arrangement of fascicles for a given site in 2D space; plot the position of the fiber nodes in 3D space; add values for electrical potential at each fiber node; and specify the waveform for each site. (4) Computational calculations. Figure 4 summarizes the methods that can be applied to analyze and interpret a simulation outcome. Single-fiber activation, population recruitment and stimulation selectivity can be estimated.
These figures represent the natural flow of steps followed when developing an HM simulation. The body of the present tutorial will consider modeling of peripheral nerves, nerve fibers and neural interfaces in a more systematic way, referring to the corresponding workflow scheme when each block is addressed.
As HM encompasses knowledge from many different domains of neurophysiology, applied mathematics and engineering, we refer to Box 2 for a list of the frequently used domain-specific terms that can be useful when dealing with HM.

Section 1: Peripheral nerve modeling
Histological components Nerve fibers are organized into macroscopic structures called nerve fascicles (Fig. 5). The nerve fibers in a fascicle are immersed in a permeable matrix of loose connective tissue called the endoneurium, which is primarily constituted by collagen fibers and supports nerve fibers and capillaries. Each nerve fascicle is surrounded by multiple intermittent cellular and collagen fiber layers that collectively constitute the so-called perineurium. Perineurial cells are connected by tight junctions, which account for the low electrical conductivity that characterizes the perineurial sheath and allows the perineurium to serve as a mechanical and chemical barrier to nerve fiber injury. Fascicles are surrounded by a connective tissue called the epineurium, which is primarily constituted by collagen fibers and adipocytes. Concurrently, the epineurium constitutes a supporting structure for the vasa nervorum, which subsequently penetrate the endoneurium in the form of a network of small capillaries. In larger nerves, fascicles are organized in bundles that are often called (major) branches, since they remain segregated inside the nerve for considerable lengths and ultimately branch into major nerves (for example, in the human sural nerve, the fascicles that will constitute the tibial, common peroneal and sural nerves can be clearly spotted well before their effective branching). The perineurium is constituted by alternating lipidic and collagen fiber layers, which anchor the nerve to the neighboring bodily structures. From a histological perspective, the paraneurium is identical to the epineurium, thus we do not distinguish between their electrical properties in our model 30 .
In the framework of quasi-static electrical stimulation 7,31 , a biological tissue can be characterized by means of rank-2 conductivity tensors. All of the modeled media are homogeneous: the perineurium and epineurium are isotropic media, while the endoneurium is anisotropic, with different transverse and longitudinal conductivities. Table 1 presents the assumed conductivity values along with the experimental references.
Implementation note: A nerve is constituted by four 3D elements that correspond to the endoneurium, perineurium, epineurium and saline bath. The addition of the saline bath addresses two concurrent issues, as it acts as the external medium needed to simulate zero potential at infinity (as explained in more detail below), and it does so by providing a reasonable approximation of the intraoperative extraneural space. The perineurial sheath is assumed to have a thickness of 0.03 times the fascicle effective diameter 32 (the diameter of an equivalent circle).
Since we work in the quasi-static approximation, the conductivity tensor is the only relevant physical quantity (Fig. 1, blocks 1a and 1b). The appropriate charge conservation boundary conditions are applied at material interfaces. The electric potential on the external boundary of the saline bath is set to zero (see 'Boundary conditions' in Section 3).
The 3D model of the nerve is obtained by extrusion of a 2D nerve section, thus no fascicle migration is modeled. This assumption is valid if the hypothesis that the implantation was performed in a region of relative topographic stability is true. The regions of topographic stability, as well as the extent of such stability, should be studied for each nerve (see 'Branching patterns' below). For example, Watchmaker et al. 33 noted that internal branching in the median nerve is infrequent over spans of a length of a few centimeters, so the assumption would hold for this region.

Branching patterns
The study of nerve branching patterns allows determination of which fiber populations are expected in each transverse section along the course of a given nerve. These considerations can be used along with surgical accessibility to choose the most suitable insertion level. A rough description of peripheral nerve branching can be found in classical neuroanatomy textbooks such as ref, 34 : such indications must be supported by quantitative descriptions of the most frequent anatomical variants for the nerve under study 35,36 .
Branching is accompanied by massive internal reorganization of the involved fibers and fascicles (internal branching) to isolate the fibers in the exiting branch [37][38][39] . Consequently, there is an evident trade-off between the need to apply stimulation distally (where functional segregation is maximal), in a region with sufficient topographic stability and in a site along the nerve course which can be easily accessed surgically.
For the median nerve, Planitzer 40 identifies the antecubital fossa as the region where the most distal plexiform  (1) Assignment of material properties to geometric components (epineurium, perineurium, endoneurium and electrode materials); (2) looping the saline bath dimensions to ensure solution reliability: an FEM problem is solved for each saline bath dimension, FEM solutions are compared and the process is repeated until increasing the saline bath dimensions does not substantially affect the solution; (3) grounding the external surface of the saline bath, and fixing the electrode currents: the total injected current (J) is normalized, resulting in active site currents (J i ) exhibiting common time behavior; (4) mesh generation and refinement: meshing is refined (smaller element size at the cost of higher computation time) if the meshing process has problems or the mesh quality is too low; (5) solution of FEM problem in terms of the electric potential calculated for each mesh node.
intermingling occurs and thus indicates the region immediately distal to it as the best place to perform functional electrical stimulation (FES).
For the sciatic nerve, Gustafson 36 claims that the most distal stimulation site should nonetheless be located proximal to the knee joint, as the presence of wires crossing the joint could lead (1) A location along the nerve branching tree is selected; this defines the fiber populations with which we will populate the nerve section; (2) corresponding samples of fiber diameters are sampled from a given distribution or from the distribution deduced from a provided sample (which can be composed of many different subpopulations, relevant for fiber computational modelling-here, in pink and blue); (3) non-homogeneous density and tortuosity can be added; (4) fibers are packed into fascicles; (5) fiber nodes (points along the fibers where the extracellular potential is imposed) are located in the 3D space; (6) the FEM potential is interpolated to fiber nodes and imposed as extracellular potential in the neural computational simulation; (7) the waveform corresponding to each active site is chosen.
to complications. The possibility of some selective muscle activation even for such a proximal location is ensured by the fact that somatotopic organization of the sciatic nerve has been shown along its entire course 41 . Intraneural fiber displacement causes the positions of nodes to differ substantially from a straight line (tortuosity).
Tortuosity is a general characteristic of nerve fibers, as it directly follows from fiber migration linked to the branching process. Moreover, it is exacerbated by the reshaping induced by neural interfaces 42 . To the best of the authors' knowledge, tortuosity is incorporated only in ref. 29 , where it is noted that it primarily affects the recruitment of unmyelinated fibers.  (1) Activation can be established via (1a) the activating function formalism (comparison of the potential nodal second difference with a threshold), (1b) the Warman-Peterson method, through simulation of a linear model of the nerve fiber and calibration via a nonlinear fiber model to establish activation conditions (in the schematics, red crosses represent points in the parameter space that do not lead to fiber activation; a similar plot clearly displays activation conditions as geometrical regions in the parameter space), and (1c) the nonlinear neural response (compartmental model in NEURON); (2) activation data for an entire fiber population can be used to compute recruitment (2a) and selectivity (2b). Recruitment indicates the population-wise number of active fibers at a given stimulation magnitude; selectivity analysis requires the assignment of a target (the muscle that the fiber stimulates or the receptor from which it receives input) to each simulated fiber (2b*) and indicates the ability of an interface to elicit the response of a given fiber group without recruiting the others. The results of a selectivity analysis can be displayed as in block (2b): the group-wise recruitments corresponding to a given stimulation condition are used to compute the selectivity related to that specific stimulation paradigm. From the recruitment characteristics corresponding to many different stimulation conditions and fiber groups, we can characterize the general behavior of a given neuroprosthetic device.
Implementation note: Branching determines the fiber populations corresponding to a given implantation site (Fig. 3, block 1). Tortuosity (Fig. 3, block 3b) is used to determine the fiber node coordinates when interpolating the node-wise extracellular potential and does not need to be accounted for in our fiber computational models. The packing of tortuous fibers assumes zero diameter (the true diameter is used in fiber response simulation but neglected as far as the arrangement of the fibers in the nerve section is concerned); otherwise, very computationally expensive flow-packing algorithms should be implemented to avoid fiber overlaps.

Topography (fascicular morphology)
Nerve topography concerns the quantitative description of nerve cross-sections in terms of the populating fascicles. The size, shape and spatial distribution of fascicles must be investigated. One example of such quantitative analysis can be found for upper extremity nerves in ref. 43 . Of course, information about functional segregation is of paramount importance in predicting the outcome of a stimulation paradigm.
The spatial fascicle density inside a nerve section, or fiber density inside a fascicle, can be imposed by weighting the random location assignment of fascicles and fibers with appropriate probability density functions f XY (x,y) (Fig. 3, block 3b). The function f XY (x,y) specifies the probability that a fascicle or fiber will occupy the point (x, y) inside a given nerve section. This probability function should be specified for each nerve at each point along its path so that, when a suitable nerve and location for an implant have been established, a statistical description of the targeted nerve slice can serve for the generation of a set of reasonable sections for parametric methods. Such probability distributions can also be referred to given subdivisions of the section (for example, quadrant divisions) similarly to the approach used in ref. 33 .
If a sequential step algorithm is used to populate the nerve section (as for the algorithms presented in 'Fiber and fascicle packing' in Section 2), then the requirement for nonoverlapping fascicles and fibers imposes that the distribution function be corrected at each step so that f XY (x,y) = 0 for all (x, y) already occupied. Box 2 | Glossary of frequently used terms General terms Hybrid modeling: modeling technique implemented in two steps: (i) volume conduction evaluation and (ii) neural response simulation. Neural response simulation: prediction of single-neuron responses (mainly in terms of spiking) caused by extracellular stimulation. Volume conduction problem: determination of the electric potential induced by a stimulating electrode in a biological structure.
Terms relating to discussion of the peripheral nerve (Section 1) (Nerve) branches: every subdivision of the nerve encountered when travelling from proximal to distal. Endoneurium: connective tissue surrounding nerve fibers within nerve fascicles. Epineurium: connective tissue surrounding fascicles. Fascicle: bundle of nerve fibers immersed in endoneurial tissue and separated from epineural tissue by a thin perineurial sheath. Perineurium: thin layer of connective tissue separating fascicles from the surrounding epineurial tissue.
Terms relating to the nerve fiber (Section 2) Distribution fit problem: determination of the probability distribution that is more likely to produce a given sample. Here, the distribution of fiber diameters is investigated (which proportion of fibers have a diameter inside a given interval). Extracellular stimulation: neuron stimulation protocol in which the driving action for neural activation is modification of the electric potential on the external surface of the cell through the delivery of a current in the extracellular medium. Nerve fiber: axons of first afferent/efferent neurons; they can be classified into Aα, Aβ, Aδ and C. Aα fibers are large-caliber (diameter 12-20 µm) myelinated fibers, Aβ fibers are smaller (diameter 6-12 µm) myelinated fibers, Aδ fibers are small (diameter 1-5 µm) slightly myelinated fibers and C fibers are very small (diameter 0.2-2 µm) unmyelinated fibers. Motor fibers will be modeled as Aα fibers, while sensory fibers will be modeled as Aβ and C fibers, accounting respectively for somatosensation and noci-and thermoception. Packing problem: determination of the arrangement of a set of objects in a region such that no two objects intersect each other or the boundaries of the given region.
Terms relating to modeling of the neural interface (Section 3) Extraneural electrode: an electrode that does not pierce the nerve structures. In this tutorial, we often refer to FINEs or FINE-like electrodes. These electrodes enhance the selectivity by reshaping the nerve topography through mechanical stress. Intraneural electrode: an electrode that pierces the nerve epineurium; intraneural electrodes can be intrafascicular if the perineurium is also pierced in order to locate one or more active sites inside fascicles. In this tutorial, we often refer to such electrodes as TIMEs or TIME-like electrodes. These electrodes can pierce nerve fascicles and thus attain very high fascicular selectivities. Invasiveness: an invasive medical procedure involves the introduction of instruments (in this case, electrodes) into the body. Invasiveness, then, is a measure of both the difficulty of the implantation and the resulting damage to physiological structures. It is evident that invasiveness is strictly linked to the final position of the electrode within the body.
Terms relating to the analysis of stimulation results (Section 4) Recruitment: the fraction of fibers activated by a given stimulation protocol. Selectivity: the ability to target specific nerve fiber subpopulations without eliciting the responses of other subpopulations that are not necessarily spatially segregated.

Morphological simplification of fascicles
Image segmentation from a histological section generally results in line drawings defining polyshape fascicles. Such polyshapes are often very irregular and non-convex, with a consequent increase in the difficulty of meshing the resulting geometries and nerve section reshaping. The reduction of such polyshapes to convex shapes completely defined through a few straightforward parameters overcomes these drawbacks.
We propose elliptical and circular fascicle approximations, which convert polyshape fascicles to equicentered, equisurface ellipses (five parameters) or circles (three parameters), respectively ( Fig. 1, block 1b). Shape simplification abstracts the relevant information from geometric features where a complete description of each shape would add too much, noninformative complexity to our model. For example, our simplification routines could alleviate the problem of the variability originating from physiological inter-subject variability and from histological manipulation-induced deformations of the nerve sample. The latter variability is always undesirable and should be carefully studied in future research. Guidelines exist, especially for the extraction and treatment of nerve biopsies to study neuropathies in the clinical context 49 , but to the best of the authors' knowledge, no study has quantified the whole deforming action introduced by manipulation of histological samples.
Finally, we note that these kinds of approximations could facilitate the definition and implementation of algorithms for FINE, such as those in ref. 20 , exploiting the fact that there exist transformations such as squeeze mappings that preserve fascicle surface areas and ellipticity.
Implementation note: The circular and elliptical fascicle approximations can be determined as follows: Suppose we have a set of regions ω i f g i representing nerve fascicles. For each region, we compute its centroid C, major and minor axis lengths (a and b), orientation α and surface area A.
We define the effective radius r eff of a region ω i as the radius of the equivalent circle and the effective major and minor axis lengths a eff and b eff of a region ω i as the major and minor axis lengths of the equivalent ellipse having the same proportions as the major and minor axis lengths of ω i , thus we obtain Then, we can choose between a circular and an elliptical fascicle approximation, given by  50 ). b1, Hematoxylin and eosin staining of rat sciatic nerve implanted with TIME (intraneural electrode) during the acute phase (0 d from implantation). b2, Immunohistochemical labeling to distinguish myelinated and unmyelinated fibers in the nerve section. The images in b1 and b2 are taken from ref. 42 .

and
C ω e Figure 6 shows the manual segmentations from three histological nerve sections from ref. 50 . It can be seen that different sections from the same nerve sample retain the gross organization into major branches, destined to the main nerves into which the sciatic nerve will branch (tibial, peroneal and sural nerves). Nonetheless, although the number of fascicles does not undergo strong oscillations, it is evident how the fascicular topography is not conserved. A good manual segmentation from histology requires some care, mainly due to the thickness of the tissue. We do not discuss these aspects here since our main purpose is to employ these segmentations to demonstrate our simplification routine. Figure 7 shows the results of the approximation and geometrical reshaping (see 'Reshaping effect' in Section 1 for the specifics of the algorithm used) of the third section. It can be observed that such simplification does not disrupt the shape or relative position of fascicles and that it produces a visually convincing result.
For the three nerve sections under analysis, we compared the meshing process of true nerve geometries (the geometries from ref. 50 have been manually segmented and are shown in Fig. 6) with the elliptical approximation (lower degree of approximation) to quantify the minimum gain in computational time, while the variation in the fascicle shape is introduced into the circular approximation (higher degree of approximation) to quantify the maximum simplification relative error in the representation accuracy.  true geometry. m ε is the median of the simplification relative error corresponding to the true-to-circular geometry change, IQR ε is the interquartile range of the simplification relative error corresponding to the true-to-circular geometry change, t true is the meshing time for the true geometry, t ell is the meshing time for the elliptical geometry, ξ is the ratio of the true-geometry to elliptical-geometry meshing time and n fasc is the number of fascicles in the nerve section.
From these data, it can be seen that the median superposition mismatch between the true fascicles and circular ones is between 7% and 14%, and the calculations take between six and eight times less time for the circular simplified data.

Reshaping effect
In case of extraneural FINE-like electrodes (for example, FINE and SC-FINE 51 ) as well as TIME-like electrodes (for example, TIME and SELINE 52 ), the implantation of neural interfaces (see Section 3) affects the intraneural topography both in the location of fascicles in the nerve cross-section and in their Step 1: polyline fascicles Step 2A: circular fascicles Step 2B: elliptical fascicles Step 2: simplified nerve section Step 3: reshaped nerve section cross-sectional dimensions according to various as yet incompletely understood mechanisms (Fig. 1, block 5).
FINE-induced reshaping. Extraneural FINE-like electrodes modify the fascicle geometries through the application of stress on the external surface of the nerve. Fascicles tend to displace inside the extrafascicular epineurial matrix and align such that stimulation selectivity is enhanced. Schiefer 20 proposed a purely geometrical reshaping algorithm that displaces fascicles without altering their surface or shape so that they do not overlap in the deformed intraneural space (Fig. 7, Step 3). In ref. 53 , a fibrotic tissue layer was observed at the interface between the nerve and electrode. This tissue can be easily added to our model as a fifth solid entity inside the nerve (see 'Histological components' in Section 1), with a conductivity of 0.1 S m -1 (refs. 21,54 ).
TIME-induced reshaping. Intraneural TIME-like electrodes induce an inflammatory response that promotes endoneurial swelling and fibrosis formation. A proposed mechanism has been discussed in ref. 42 , where data relative to such phenomena in the rat sciatic nerve can be found. It is shown that only the pierced fascicles (local endoneurial environment) undergo swelling. We propose the following two-step procedure to simulate intraneural electrode-induced reshaping and apply it in the simple case of the rat sciatic nerve, employing Wurth's data. From the time curves of the local endoneurial environment and of the fibrotic tissue area, the amounts of overall swelling and fibrosis can be obtained. Then, local endoneurial environment fascicles are gradually swollen one at a time until the found surface area is reached, and the other fascicles are displaced according to Schiefer's algorithm 20 until no intersection is obtained and the swollen fascicle surface corresponds to the data of the assumed model. Then, when defining the FEM geometry, a fibrotic tissue layer (σ = 0.1 S m -1 ) can be introduced inside the local endoneurial environment fascicles, which will alter the electric potential field inside the fascicles. The complete procedure can be seen in Fig. 8.

Section 2: Nerve fiber modeling
Overview of cable-like models Nerve fibers are the axons of the afferent or efferent neurons that directly contact the periphery. Here, we only employ lumped-parameter cable-like electrical equivalent models of nerve fibers 2,55 . These models belong to the family of multicompartment models, in which each neuron is divided into compartments that are simulated independently. Cable-like multicompartment models include longitudinal resistive paths that connect units, representing different transverse layers of neuronal structures. Typically, in cable-like models, only a portion of the cell axon is modeled, and activation is estimated from the excitation and propagation in that portion 20,21 .
The simplest cable-like models are single-cable models where the nerve is considered as a single unit; for this reason, they are currently the preferred choice for modeling unmyelinated fibers 29, 56 and have also been used in the past for myelinated fibers 57 . The relevance of periaxonal conduction phenomena in myelinated fibers requires the introduction of multi-cable models 58 .

The MRG fiber model
The current gold standard for myelinated fiber modeling is MRG 59 . The MRG model is a double-cable conductance-based compartmental model based on an architecture first introduced by Blight 60 and incorporating the paranodal ultrastructure 61 . Non-nodal compartments are modeled as passive compartments, whereas nodal compartments contain active elements in the form of Hodgkin-Huxley-like conductances 62 . The outer  cable represents the myelin sheath and is composed of passive compartments (see ref. 63 for alternatives). The characteristic equations for nodal ion channels can be found in ref. 59 and, for humans, in refs. [64][65][66] . The modularity of the defined fiber model can be readily adapted to other animals (for example, in ref. 24 , the rat motor axon nodal dynamics from refs. 67,68 was used) and other fiber types (for example, in refs. 69,70 , the sensory axon nodal dynamics from ref. 71 was used).

Activation threshold estimation methods
The simulation of the nonlinear membrane dynamics of a large number of nerve fibers in different stimulation conditions is currently the true bottleneck of the optimization process. This has been made clear, for example, in ref. 20 , where the simulation of a number of fibers on the order of millions (corresponding to the high number of stimulation parameter sweeps for a single nerve) required~90 d.
Alternatives exist that estimate fiber activation without solving the complex, highly nonlinear membrane dynamics or that deduce with a single run of the nonlinear solver some kind of calibration curve that can be used to determine from direct comparison whether a fiber will be active under a given stimulation. Such methods can be divided into activating function methods (Rattay's activating function 3 and Joucla's mirror estimate 72 ), passive-cable methods (Sweeney's passive cable 73 ) and hybrid methods that employ weighted forms of activating functions deduced from passive-cable theory and that employ one run of the nonlinear fiber model to establish an activation threshold 74,75 (the transition from Warman's to Peterson's model is justified in ref. 76 ). Activating function methods (see Fig. 4, block 1a) are essentially algebraic. They consist of thresholding specific quantities that possess a direct, easily computable relation with the extracellular potential (its second differences and its opposite for the activating and the mirror functions, respectively), and they are thus very computationally efficient; nevertheless, they provide only relative recruitment order. On the other hand, the Warman-Peterson method (see Fig. 4, block 1b) requires the determination and simulation of equivalent linear models subject to the given extracellular potential: it thus requires a one-time simulation of the preferred nonlinear model describing the nerve fiber in order to calibrate the intercurrent relation between the nonlinear model and its linear approximation with respect to a number of varying parameters that could modify such correspondence.

Fiber diameter distribution
It is well known that different fiber types are characterized by different typical diameters 77 , which relate to the conduction velocity 78 and fiber excitability 79,80 . Development and physiological alterations affect fiber diameters, consequently changing their properties. Consequently, both clinical practice and biomedical engineering would benefit from further studies of fiber diameter distributions, which can characterize the physiological state of a subject and help them act accordingly under FES.
Samples are usually obtained from electron microscopy (EM), which allows invasive high-resolution imaging of small regions of nerves. An entire branch of research has been developed for the automatic segmentation of EM images; we cite here three open-source libraries for segmentation: the method described in ref. 81 , AxonSeg 82 and AxonDeepSeg 83 .
Methods for noninvasively measuring fiber diameters from diffusion magnetic resonance imaging (MRI) were introduced in ref. 84 . Actually, inference procedures quantifying microstructural neural features from diffusion MRI have received attention to noninvasively investigate the central nervous system, but, to the best of the authors' knowledge, no further developments or guidelines have been published to investigate peripheral nerves.
The efforts in fitting the diameter probability distribution function have moved both in the direction of choosing distribution families that are theoretically more appropriate than a simple Gaussian (gamma, beta or extremum value distribution 85,86 ) and toward the development of automatic routines for the individuation of the best mixture distribution 87 representing the sample (see ref. 88 , where the Bayesian information criterion (BIC) 89 was used).
When the optimal number of mixture components and the optimal parameters for each component have been found, different distribution family assumptions can be compared through the use of goodness-of-fit techniques such as qq-plots and sample density versus theoretical density plots. Finally, the underlying fiber type can be recovered since different fiber types have in general different typical diameters 77 . Our proposed workflow is shown in Fig. 9; results with respect to the dataset shown in Fig. 10 are presented in Table 3 and Fig. 11.

Fiber and fascicle packing
The spatial arrangement of fibers inside a nerve affects their activation since it directly determines their extracellular potential 2 ; in turn, it depends on the distribution of fascicles and their subdivision in functional areas.
It is common practice to populate fascicles obtained by manual segmentation of histological sections. Nonetheless, we may be interested in the ex novo generation of plausible fascicular geometries (Fig. 1, block 1a).
Both the problems of populating a nerve with fascicles of given dimension and fascicles with appropriately distributed fibers are packing problems, which are challenging to tackle in their general form. We describe two algorithms that produce acceptable suboptimal results for the two processes, while requiring reasonable computational time (on the order of seconds for naturalistic structures) and memory resources for the cases at hand (Fig. 12).
A priori check for intersections (ACI) algorithm. Initialization. We build a square grid with step ε completely covering Ω. Fascicles are ordered in decreasing diameter order to enhance the probability of successful termination.
Step. Suppose that we have already located some fascicles, having centers {C j } j and diameters {d j } j . The next fascicle (diameter d) is centered at a random point C among the grid points that are at least (d i + d)/2 away from each C i and at least d/2 away from each point of the nerve contour.
Termination. If a fascicle cannot be placed, then the procedure restarts with a different placement of the first fascicle. If the procedure fails in more than a given number of attempts, the diameters are resampled. Since termination cannot be ensured a priori, we studied the fraction of failed attempts and the packing time with respect to fascicular shape (the results are shown in Fig. 13).
Minimum distance between fascicles. A minimum distance δ between fascicles can be specified to improve meshing. The steps above remain unchanged if each fascicle is supposed to have a diameter of d + 2δ.
Multiscale square grid algorithm. Nomenclature. We define a sequence of lengths l n f g S þ 1 n ¼ 0 , where S is the desired number of grid refinements (scales) to be performed such that: and a partition D n f g S n ¼ 0 of the set of fiber diameters such that D n ¼ d i jl n þ 1 d i l n f g : We denote by Q n ð Þ i the ith square cell of a uniform grid with step l n .
Initialization. We build a square grid with step l 0 completely covering Ω. We then obtain the set of M square cells Q Step. For k 2 0; S ½ , we perform the following actions: We performed packing of fiber diameter samples from the distribution fit results displayed in Fig. 11 Fig. 9 | Complete workflow of the fiber diameter distribution fitting routine. For each distribution family (in this paper, we studied normal, gamma and beta distributions, as we find that they are particularly appropriate for fiber diameter samples), the optimal number of components is selected through the BIC. Then, through qualitative goodness-of-fit (GOF) assessment techniques, the best fitting mixture is selected. Such distributions can be used as sampling distributions for fiber diameters in a given simulation, and the corresponding samples can be packed in the relative nerve fascicles. From the characteristics of each of the obtained components, the corresponding fiber type can be deduced, which determines the most suitable computational model for each modeled fiber.

NATURE PROTOCOLS
fibers from 1,000 to 10,000 (the actual axon count per millimeter for a subject over 20 years old is~30,000 (ref. 90 )). The packing was always admissible and was always accomplished in <0.5 s. It is noteworthy that multiscale square grid (MSG) produces almost random bidimensional packings for reasonable packing fractions (see Table 4). We assumed that statistically random packing was the strictest requirement to test the quality of the result of our algorithm and that any other nonrandom assignment could be easily managed through adequate density functions.

Section 3: Neural interface modeling
PNS is intrinsically linked with the concepts of selectivity and invasiveness 91 . Selectivity is the ability to target specific nerve fiber subpopulations without eliciting the responses of other subpopulations that are not necessarily spatially segregated. Invasiveness corresponds to the regime position of the electrode, which indicates both the difficulty of the implantation and the resulting expected damage to physiological structures. As a general rule, electrode selectivity increases with invasiveness. Comprehensive reviews of the currently used peripheral nerve interfaces can be found in refs. 91,92 . As a first approximation, electrodes are bi-domain objects: they include active sites, which inject independent currents into the nerve, and an electrode body, which is usually developed using an isolating material.
In FEM, electrode geometries must be simplified, and a minimal set of geometrical parameters must be defined. We propose to build a set of proper electrode geometry simplifications in terms of simple primitives to obtain a good trade-off between the fidelity of the model to the represented objects and the generalization capability, which will help with rapid adaptation of the studied geometries to unexplored contexts and applications.
Proposals for the parametric simplified models in terms of a predefined set of a few straightforward parameters for the TIME 93 and the FINE 94 are outlined in the next subsection.

Electrode parameterization
In the following, we introduce a parametrization for the simplified 3D geometries of TIMEs and FINEs. This allows the automatic generation of electrode geometries given a few interpretable geometrical parameters. Examples of the resulting simplified geometries can be seen in Fig. 14. TIME. The geometry of the TIME is described in the electrode reference frame S E by symmetric extrusion for a length h b /2 along the Y axis of the 2D polygonal shape x : z : The electrode is composed of a body (* b ) and a tip (* t ): the body is a parallelepiped with dimensions [l b , h b , w b ], while the tip attaches to a w b × h b face and has a parallel face with dimensions w t × h b and a length l t .
Circular active sites (diameter d as ) are built, whose center coordinates are where n represents the identification number for each active site. Each active site extends for a length h as from the electrode body surface.
Fixing the insertion geometry of an intraneural electrode corresponds to determining its pose relative to the nerve reference frame (six degrees of freedom). We assume that our stimulation occurs in the mid-plane of the nerve. Additionally, we impose that the electrode main shaft should be normal to the XY plane. This reduces the dimensionality of our problem to three. The insertion is fixed by the position of the center of one particular (leading) active site and the angle of rotation between the electrode and the nerve reference frame (Fig. 1,  block 4).
FINE. The body of the electrode is constituted by a parallelepiped with side lengths [(a + t), (b + t), l] with a parallelepiped hole with side lengths [a, b, l]. Active sites are in positions where l cc ¼ a= n as 2 þ 1 À Á and n as is the number of active sites. Each active site has diameter d as and depth h as .
Between the electrode internal surface and the nerve section, a saline bath can be inserted as in ref. 21 , where the FINE stressinduced reshaping effect can be neglected. When moderate reshaping effects are taken into account, as in ref. 20 , it can be assumed that the epineurium completely fills the internal space of the electrode. The large mismatch between the conductivities of the epineurium (σ = 0.083 S m -1 ) and saline (σ = 2 S m -1 ) affects the intraneural electric field distribution, and thus, for each implant, the best hypothesis should be carefully chosen and tested.

Boundary conditions
The potential map solves the Poisson equation with adequate boundary conditions. The most widely used boundary conditions are those of Dirichlet and Neumann, which specify the value of the electric potential or its flux, respectively, relative k: number of components; P: number of used parameters; ln L: log likelihood of the resulting distribution; BIC: Bayesian information criterion value, which is lowest for the best fitting distribution; ΔBIC: difference between BIC values for the corresponding distribution and bestfitting distribution. The best-fitting mixture has k = 2.
to selected subdomains 21,95 . Typically, the Dirichlet condition of zero potential (ground) is imposed at infinity, while the Neumann conditions specify the value of the current density on the active site surfaces (Fig. 2, block 3).
Since FEM discretization does not allow the Dirichlet ground condition to be imposed at infinity, it is normally imposed on the outer surface of a large saline bath that surrounds the involved structures. The dimensions of such a solid can be found while imposing the condition that the potential field inside the nerve with larger saline baths does not vary substantially 21 (see Fig. 2, block 2, and 'Measures of field similarity' in Section 4).
Neumann conditions at the active sites can be substituted with Robin boundary conditions, mentioned in ref. 96 and thoroughly investigated in refs. 95,97,98 . They can be used to take into account the surface conductance of electrodeelectrolyte interfaces but in general must be estimated case by case. They can be simplified by assuming that the potential drop at the interface is high, in which case non-homogeneous Neumann boundary conditions can be used, which require only knowledge of the magnitude of the stimulation current and the surface of the active site of stimulation 97 .
In ref. 95 , special boundary conditions representing nonlinear nerve fiber dynamics are presented in the framework of a full FEM approach to FES.

Stimulation waveform
With the assumption of quasi-staticity, the extracellular medium is purely resistive, and the extracellular potential at each node follows the same time-dependent waveform. Exploiting this hypothesis, the action of more electrode active sites injecting current according to the same waveform (or waveforms that differ only by a multiplicative constant) can be obtained through a single FEM simulation (Fig. 2, blocks 3-5). In general, the number of FEM simulations needed corresponds to the number of independently driven active sites during stimulation.
The NEURON framework constituted by the extracellular and xtra mechanisms allows the complete definition of custom stimulation waveforms. Concurrently, they rely on heavy approximations, and the choice of the stimulating waveform should be made with care. Specifically, it should be clear that our framework assumes electromagnetic quasi-staticity, which critically depends on the waveform frequency content 7 . Finally, safety limits must be taken into account when developing new In the first step, all fibers between d max and d max /2 are packed; in the second step, all fibers between d max /2 and d max /4 are packed; and in the last step, all remaining fibers are packed; right: illustration of the ACI packing algorithm for a three-fascicle nerve section. (Bottom) left: packing of 10,000 fibers (myelinated Aβ and Aδ, and unmyelinated C fibers, following the empirical distributions described in ref. 90 ) in a fascicle of radius 0.5 mm; right: packing of 18 fascicles in a nerve section of a 3-mm radius (dimensions are compatible with the lower median nerve described in ref. 43 ).
stimulation strategies in silico in order to maintain the translatability to in vivo experiments. An up-to-date treatment of such limits can be found in ref. 99 .

Section 4: Analysis of stimulation outcomes
Measures of field similarity The shape of the FEM electric fields alone provides some information that can offer low computational burden hints regarding the results of stimulation. Figure 15 shows FEM potential maps that highlight the spatial information provided by the FEM electric field shape. For example, the volume of activation methods, most notably the activating/mirror function formalism, rely only on the resulting potential maps and underscore the need to define proper field similarity measures.
We can compare two fields through the indices introduced in ref. 100 and then used in ref. 101 for electroencephalography/ magnetoencephalography reconstructed source comparison, which have been used in refs. 21,24 in the context of the nerve stimulation volume conduction problem solved by FEM. We refer here to Fig. 5 from ref. 21 , which highlights the very different electric field shapes generate by a FINE and a TIME.
One alternative is to use the above measure of similarity to compare the activating functions (second differences of the field) instead of the fields. This method offers the advantage of employing a quantity with a well-established neurophysiological meaning.

Recruitment curves
Recruitment curves (see Fig. 4, block 2a) represent the fraction of active fibers for a given magnitude of the stimulation 102 and are used mainly with three concurrent issues in mind: threshold identification, recruitment steepness and recruitment order.
Recruitment is directly linked to motion or sensation, and thresholds for palpable muscle contraction are generally set to a recruitment of 10% of the fibers in a population 103 . On the other hand, it has been shown that the recruitment of a single sensory fiber can in some cases be perceived by subjects 104 .
Recruitment steepness (gain) pertains to the ability to control recruitment with varying stimulation. We identify low maximum gain and high linearity as two desirable characteristics of a recruitment curve, which facilitates controllability and simplifies the design of the relative controller 105,106 .
Recruitment order concerns the dependence of activation on the fiber diameter. Specifically, natural recruitment of motoneurons follows a size principle 107,108 in that small-calibre fibers are recruited first. Artificial inverse recruitment (caused mainly by extraneural stimulation 79,102 ) leads to muscle fatigue and does not allow fine muscle control 109 . The recruitment order can be determined through the comparison of recruitment curves (or thresholds) for populations of fibers with different diameters. From ref. 25 , we refer here to Fig. 9, where differential recruitment with respect to fiber diameter is shown, and Figs. 6-7, where the validity of hybrid modeling for invasive PNS was studied by comparing modeling and experimental recruitment curves for a rat sciatic nerve stimulation application.

Selectivity indices
Selectivity corresponds to the ability of a given system to elicit the response of a given subpopulation of fibers without affecting the remaining population (see Fig. 4, block 2b). We employ two indices to estimate the selectivity of a neuroprosthetic device 24 .
Suppose that our fibers are partitioned into m groups and that during a stimulation procedure, a fraction 0 ≤ μ i ≤ 1 of fibers in the ith group have been recruited. Thus, we define the group-wise selectivity of the given routine through the index which lies in the range [−1, 1] and is 1 if μ i = 1 and μ j = 0 8j ≠ i. This index penalizes coactivation of different fiber groups, which can result, for motor fiber stimulation, in muscle fatigue or in the impairment of the movements to be restored through activation of antagonist muscles; for sensory fiber stimulation, co-activation is equally detrimental since it leads to paresthetic sensations. In ref. 21 , the above selectivity index was coupled with the selectivity index where n i is the number of active fibers from the ith fiber group, which can be used to ensure a sufficient activation level in a specific fiber group. By combining the selectivities corresponding to stimulation procedures s performed with the same implant through the index we can characterize the selectivity of an implant (through the maximal group-wise selectivity) and compare it with different implants employing the same electrode (for example, to test the best implantation site for the electrode) or with implants employing different electrodes (for example, to compare the selectivity of different electrode families or designs). The most adequate combination of selectivity indices must be chosen for each situation: for example, one might seek to combine, for a number of specific stimulation paradigms, the indication of specific activation of a given group i of fibers through Sel i with the indication of sufficient activation of that group through S i to ensure that the stimulation sufficiently stimulates the fiber group to transmit a reliable sensation/ motion. Concurrently, one might seek to study the selectivity ability of an entire implant: a high value of S T , in fact, ensures that different stimulation paradigms can be found that can elicit a selective response from all fiber groups (in an average sense).

Our results
We presented a modular, customizable and expandable framework for hybrid modeling of invasive PNS interfaces. For each step, we provided a short review of the biological evidence, along with a list of several different conceptual and computational approaches.
The problems of fiber diameter distribution fitting and consequent type separation have been tackled. Such problems have never been fully addressed but have become relevant with the introduction of different fiber types into the HM framework. Two algorithms have been proposed, allowing almost random packings of fibers inside fascicles and fascicles inside a given nerve section. They require several seconds for each run, even for naturalistic packing fractions, and thus do not substantially increase the solving time for the corresponding model.
A set of possible approximations for nerve geometries has been introduced to provide the right level of abstraction for each kind of scientific question. Exploiting such approximations, we define a completely general framework that can be used to provide the adequate statistical power required by in silico investigations, allowing us to produce (owing to the implementation of ex novo generation from statistical characterization of fascicles) and solve for (owing to computational time reduction) a great number of comparable naturalistic geometries.
Geometrical reshaping algorithms have been introduced in our workflow, which will allow the integration of new data and models regarding the neural tissue response to neural interfaces implanted both intra-and extraneurally.
Most notably, we introduced a hierarchy of methods in which different approximations can be used to answer different kinds of scientific questions ranging from the planning of single implants to the optimization of a new neural interface and general research questions regarding neural tissue physiology.
Below we briefly describe the limitations of our model, future developments and future research routes. While limitations refer to inherent constraints of our framework, future developments address the analysis of a number of aspects that can be easily introduced in the present framework but that our group has never openly addressed. Finally, the 'Future research routes' subsection investigates applications of neuroprosthetic interest directly following from future developments of the model.

Limitations
Physics assumptions. Several limitations and caveats must be underscored when choosing an HM framework. Strong assumptions are made with respect to the physical characteristics of materials, most notably linearity, homogeneity, frequency independence (quasi-static application of Maxwell's equations) and time invariance. An HM is based on the strong assumption of decoupling between the volume conduction and neural response modeling, which inherently excludes all phenomena involving ephaptic coupling and electrodiffusion.
2.5D modeling. The proposed model does not aim at a full 3D reconstruction of nerve structures: instead, our nerve geometries are obtained through extrusion of single histological sections. Addressing 3D reconstruction of nerves would be a major improvement of our framework, but a number of developments are required that have never been tackled by our group. For example, (semi)automatic segmentation and relative alignment of histological sections would be needed, as well as an automatic way to detect and manage internal branches while maintaining a tractable computational burden of FEM.
Further validation data. There is a published example of a general methodology for hybrid modeling validation for invasive PNS for rat sciatic nerve stimulation 25 . The approximations described in this paper were not explicitly validated: their validity was inferred from comparison with the outcome of the most detailed variant. A similar approach has also been used in ref. 27 , where stimulation via microneurographic needles and via TIMEs were compared without an explicit experimental validation. Specifically, a future validation of HMs in humans would be indirect and can only confirm very general claims in terms of global selectivity (see ref. 25 on generalization), as it requires individual nerve morphologies, which currently can only be determined invasively. However, in the future, new imaging techniques could provide additional data, at least in animal models.
Ways in which the model could be further developed Evaluting the effects of steering and multisite stimulation. We refer to Fig. 8 in ref. 21 , where it is shown that bipolar (anodecathode) stimulation allows modulation of the recruitment of different nerve fascicles, thus allowing selective targeting of a higher number of fascicles. A full evaluation of the advantages of combinations of steering, bipolar and multisite stimulation has not yet been carried out but would be helpful to provide an idea of the true selectivity capabilities (for example, spatial selectivity enables leveraging of the multimodal potential field maps that can be elicited via superposition) of an implant, which can be highly underestimated without carrying out an extensive search in the space of controllable parameters.
Varying the waveform. The influence of stimulating waveform parameters on neural excitation has not been studied so far: optimizing such parameters could pave the way for selective activation of fiber populations of different types/ diameters and is one of our targets in the next future.
Including C fibers in the model. Small unmyelinated C fibers, responsible for thermo-and nociception, can be incorporated into HMs. Our group has not studied them, but ref. 29 , for example, showed that the C fiber model in ref. 56 subject to extracellular stimulation had a much higher threshold with respect to the myelinated MRG fiber model 59 .

Future research routes
In the spirit of several initiatives begun in recent years, such as ModelDB 110,111 or the Human Brain Project (https://www. humanbrainproject.eu/en/explore-the-brain/share-data/), we advocate the creation of a shared knowledge open database to collect and unify data on the morphology of peripheral nerves and nerve fibers. In the present paper, we suggested a series of assumptions and approximations that ease the computational analysis of PNS: such assumptions could provide the fundamental backbone for such a project, outlining the set of information readily obtainable from morphology studies which would be useful in a computational setting. For example, an analysis of the expected selectivity for an implant in a given location of the median nerve could be performed once location and effective diameter distributions of the fascicle for a sufficiently high number of histological samples are known. From the probability distributions, a great number of simplified circular-fascicle stereotypical nerve sections can be built, on which selectivity analysis can be performed. The selectivity results corresponding to the different stereotypical nerve samples (extruded from the sections) can be aggregated to provide a hint of the expected selectivity which is robust with respect to inter-subject variability.
Finally, we here list a few fields of potential application of HM, which would greatly benefit from the expansions of the model listed in the previous sections and would definitely allow HM to be regarded as a fundamental tool in the neuroprosthetics research field: 1 Addition of new fiber types to simulate, for example, the response of autonomic unmyelinated nerves for bioelectronic medicine applications 112 . Moreover, investigating unmyelinated fibers could pave the way for more naturalistic stimulation protocols exploiting the interoceptive action of such fiber types 113 . 2 Characterization of the changes in the response in patients suffering from neuropathies that affect modeled characteristics of fibers. For example, demyelination in diabetes 114 can be simulated by modifying the fiber geometry in the MRG model, and small fiber neuropathies 115 can be simulated by modifying type-wise fiber diameter distributions. 3 Addition of physical-physiological grounded models of the tissue response to neural implants 116 to obtain reliable geometries from the application of a given neural interface to a nerve of interest whose section has been specified, for example, through histological techniques. 4 Development of new biomimetic encoding strategies 117 for comparison of the peripheral activation corresponding to each strategy and to provide a first offline benchmarking. 5 The possibility of linking the model to higher neural pathway processing models to simulate the modulating action of higher centers (for example, spinal ganglia 10 ). 6 Use of reciprocity to develop recording models from the current stimulation framework [118][119][120] . We believe that the use of HMs could allow the development of more effective PNS neural interfaces and therefore facilitate their translational use for clinical applications.

Data availability
The authors declare that all data needed for the simulation examples in this tutorial can be found within the paper and its references. No new experimental data have been used for the writing of this protocol.

Software availability
Illustrative code for the COMSOL/MATLAB setting of HMs, for fascicle shape simplification and for fiber and fascicle packing is available at https://github.com/s-romeni/PNS-HM.