Modulating individual axons and axonal populations in the peripheral nerve using transverse intrafascicular multichannel electrodes

Objective. A transverse intrafascicular multichannel electrode (TIME) may offer advantages over more conventional cuff electrodes including higher spatial selectivity and reduced stimulation charge requirements. However, the performance of TIME, especially in the context of non-conventional stimulation waveforms, remains relatively unexplored. As part of our overarching goal of investigating stimulation efficacy of TIME, we developed a computational toolkit that automates the creation and usage of in silico nerve models with TIME setup, which solves nerve responses using cable equations and computes extracellular potentials using finite element method. Approach. We began by implementing a flexible and scalable Python/MATLAB-based toolkit for automatically creating models of nerve stimulation in the hybrid NEURON/COMSOL ecosystems. We then developed a sciatic nerve model containing 14 fascicles with 1,170 myelinated (A-type, 30%) and unmyelinated (C-type, 70%) fibers to study fiber responses over a variety of TIME arrangements (monopolar and hexapolar) and stimulation waveforms (kilohertz stimulation and cathodic ramp modulation). Main results. Our toolkit obviates the conventional need to re-create the same nerve in two disparate modeling environments and automates bi-directional transfer of results. Our population-based simulations suggested that kilohertz stimuli provide selective activation of targeted C fibers near the stimulating electrodes but also tended to activate non-targeted A fibers further away. However, C fiber selectivity can be enhanced by hexapolar TIME arrangements that confined the spatial extent of electrical stimuli. Improved upon prior findings, we devised a high-frequency waveform that incorporates cathodic DC ramp to completely remove undesirable onset responses. Conclusion. Our toolkit allows agile, iterative design cycles involving the nerve and TIME, while minimizing the potential operator errors during complex simulation. The nerve model created by our toolkit allowed us to study and optimize the design of next-generation intrafascicular implants for improved spatial and fiber-type selectivity.


Introduction
Recent technological advances have expanded our understanding of how an electronic device interacts with the peripheral nervous system (Chew et al 2013, Christie et al 2017, Koo et al 2018, Zhang et al 2019, Cho et al 2020. The extra-neural cuff electrode (Naples et al 1988, Tyler and Durand 2003, Xue et al 2015, which wraps around the nerve, is a well-established method for interfacing with nerves. This has the advantage of being less invasive than penetrating sharp electrodes but has issues of lower recording signal-to-noise ratio and limited stimulation resolution, owing to the large electrodes placed distant (outside) to the nerve (Ortiz-Catalan et al 2013). Alternatively, newer strategies involve placing microelectrodes inside the nerves (Rijnbeek et al 2018). These arrangements may offer better deviceto-biology interfacing than the clinically established cuff electrodes. A particular implementation of these in-the-nerve neural interfaces is the transverse intrafascicular multichannel electrode (TIME) (Boretius et al 2010, Kundu et al 2014, Badia et al 2016. While the exact implementation varies between research groups, they have the common characteristics of incorporating multiple, micron-scale electrodes in a mechanically flexible substrate. For example, the centimeter-sized, soft array is inserted through the nerve using a removable guiding device (Kundu et al 2014).
A major premise of TIME is that these soft/ flexible electrodes are less invasive than the conventional metallic, penetrating electrodes like the Utah slanted electrode array (Branner et al 2001) while offering similar stimulation and recording capabilities. Another possible advantage of TIME device is improved spatial selectivity during electrical stimulation (Boretius et al 2010). For instance, a comparative study between TIME, longitudinal intrafascicular electrode and multipolar cuff electrodes demonstrated that TIME provide deeper excitation to the implanted fascicle, lower threshold for muscle activation, and higher selectivity for the tested muscles (Badia et al 2011b). Leveraging these possible advantages, the STIMEP portable multichannel stimulator was proposed to optimize selectivity by automatically delivering sequences of complex stimulation patterns for up to four electrodes simultaneously (Guiho et al 2021). TIME has also been shown to provide nearnatural sensory information to an amputee when controlling a dexterous hand prosthesis (Raspopovic et al 2014).
The in vivo and clinical performance of cuff electrodes have been assessed in the vagus nerve (Bucksot et al 2019), arm nerves such as the radial and ulnar nerves (Tan et al 2015), and leg nerves such as the sciatic nerve (Xue et al 2015). Furthermore, computational simulations are well-established methods for optimizing the design of cuff electrodes and for investigating new stimulation approaches (Kent and Grill 2013, Pelot et al 2018a, Eiber et al 2021, Musselman et al 2021. In contrast, less is known about the TIME arrays' performance. In addition, many stimulation waveforms have been reported to improve stimulation efficacy including reduced stimulation charge (Yip et al 2017), improved fibertype selectivity Butera 2014, Pelot andGrill 2020), and minimizing undesirable physiological side-effects (Petruska et al 1998, Ackermann et al 2011. However, the benefits of these waveforms in the high-resolution TIME arrays (compared to conventional cuff electrodes) remains unclear. This study presents a model-based analysis on the performance of TIMEs placed in the sciatic nerve, and more generally in the peripheral nerve, in preferentially modulating different types of fibers and in the context of non-conventional stimulation parameters.
This work began with the development of a software toolkit to simplify and automate the process of creating and using hybrid nerve models with TIME setup in the NEURON/COMSOL environments. NEURON and COMSOL are widely used and thoroughly vetted by the neuroscience and engineering communities. Our toolkit takes as input a nerve model created in the NEURON simulator and automatically transfer the geometry to COMSOL, permitting users to deploy pre-existing, validated, and highquality models in the literature (McIntyre et al 2002, Sundt et al 2015. This eliminates the need to re-create and re-validate models from scratch. The toolkit also computes the extracellular potentials induced by electrical stimulation within the COMSOL environment, then transfers the results back into NEURON for simulating nerve responses. With the aid of this toolkit, we investigated how axons responded to electrical stimulation with different TIME configurations in a sciatic nerve with 1,170 myelinated (A-type, 30%) and unmyelinated (C-type, 70%) nerve fibers, grouped into 14 fascicles. The TIME model incorporated realistic attributes such as an insulating polyimide substrate with embedded platinum electrodes. By testing a variety of electrode placements such as monopolar and hexapolar (Wong et al 2007), and stimulation waveforms (kilohertz stimuli and cathodic ramp), we assessed clinically important considerations such as spatial selectivity of fiber recruitment, preferential activation and blocking of different fiber types, minimizing inadvertent excitation of non-targeted axons, and removal of undesirable onset responses.

Methods
To demonstrate the operation of our software toolkit, here we describe how it automatically created a sciatic nerve, and how it synchronized the extracellular potential data between NEURON and COMSOL following electrical stimulation.
Most of the current approaches for implementing hybrid models requires manually building the same nerve twice, first in NEURON, then again in COMSOL ( figure 1(d)). The COMSOL model provides fine meshes of the tissue environment and calculates the extracellular voltage generated by the stimulator. Such voltage values are exported from COMSOL as a data file, imported into NEURON and interpolated to the corresponding fiber sections. Finally, the NEURON software computes the fiber responses based on cable equations. Such processes are laborious and error-prone, due to significant operator interventions between steps. Simplifying and increasing the speed of modeling nerve fibers using our automated toolkit. (a) Nerves are structurally complex, with constituent tissue having different biophysical properties. To achieve high accuracy and to reduce computation time, a hybrid simulation approach is common. This combines the results of Hodgkin-Huxley (HH) cable-equation models (b) and finite element method (FEM) models involving Maxwell's equations (c). (d) Conventionally, the hybrid approach requires repeated model construction and manual data synchronization between platforms. For example, the extracellular voltage is calculated and exported from COMSOL and imported back to NEURON. (e) Our toolkit automates such processes of model reconstruction and data imports/exports.
Our toolkit eliminates these issues through custom software combining NEURON's Python interface and COMSOL's MATLAB interface (figure 1(e)). The toolkit uses these back-end interfaces automatically, thereby shielding users from the underlying programming complexity. Specifically, once users build the NEURON model, the software automatically extracts the neuronal information via NEURON's Python interface, writing them into COMSOL's MATLAB interface to generate an equivalent model in COMSOL. Then, the COMSOL model is meshed and a finite-element model with boundary conditions appropriate to the stimulation electrode configuration is solved to compute the extracellular voltage values. These are automatically exported from COMSOL and imported back to NEURON as a data file. Using Python's interpolation tool, the electric potential is calculated for each fiber section. Users can then apply NEURON's fast computation to simulate fiber responses but at the same time with an accurate representation of the extracellular environment.
Below we begin by describing our core modeling steps. We then elucidate the model parameters used in our demonstration nerve. The computer used in this work has two 8-core Xeon processors running at 2 GHz, and with 128 GB RAM. We used NEURON version 8.0.2, COMSOL version 5.6, MATLAB version R2021a and Python version 3.9.8.

Step 1: NEURON model construction
The starting point for the process is a fascicle/ nerve model already implemented in NEURON (figure 2(a)). This is rarely a barrier with conventional neurophysiological or neural engineering studies. For users in possession of experimental data, for instance, via anatomical staining and microscopy studies conducted in their labs, these data would first be imported into NEURON using the software's built-in capabilities. For users without experimental data, many detailed models of fascicles/nerves are available via curated, online databases such as ModelDB (Marenco et al 1998). To illustrate the latter case, our sciatic nerve model simulated A and C fibers based on the published MRG model (McIntyre et al 2002) and T-junction model (Sundt et al 2015). We chose rat as the modeled species as it is well-studied and mostly used in preclinical peripheral nerve stimulation studies. To better match the characteristics of A and C fibers in rat sciatic nerve, we modified the following parameters from both models: the A fiber diameter of the MRG model was changed from 10 to 5.7 µm (Rich and Brown 2018), and the C fiber involved only the linear peripheral axon of Sundt's model , whose diameter was changed from 0.8 to 0.43 µm (Smith et al 2012). The length of both fibers was set to 1 cm. The number of compartments was set to be 221 for A fibers, the same value as the MRG model, and 100 for C fibers. The sciatic nerve model was created using NEURON's .hoc scripting language. It was constructed based on the toluidine-blue-stained, cross-section (figure 2(a)) of a sciatic nerve obtained at mid-thigh of a rat (Burks et al 2014). The cross sections of nerves and fascicles were approximated to be circles, with cross-sectional area manually made close to the size of the original histological data to generate the NEURON morphology. The sciatic nerve contained three sets of fascicles-3 large (1 mm diameter), 5 medium (0.6 mm diameter), and 6 small (0.3 mm diameter). Each fascicle involved both A (30%) and C (70%) fibers to conform with the anatomy (Castro et al 2008), resulting in a total of 1,170 fibers bundled together to form the whole sciatic nerve. The placement of each fiber in the transverse plane was randomized using values drawn from a uniform distribution to account for biological variability. While we are not able to determine the precise location of each fiber from the light-microscopy-based image, histology data based on electron microscopy will readily allow for exact replication of fiber positions. Nevertheless, our approximated positioning is inconsequential for the toolkit's operations.
NEURON provides efficient computation of HH-based ionic channel kinetics by approximating equation (1) (Carnevale and Hines 2010): where s is the arc-length distance along the neuron, ρ i is the intracellular resistivity of the nodal elements, V i is the intracellular potential, β is the reciprocal of the nodal radius, C m is the membrane capacitance per unit area, V m is the membrane potential, and J ion is the total ionic current density. For multicompartmental models, equation (1) is repeated for each section to associate with its ionic properties. Segment resolution thus plays an important role in determining neuronal response. Generally, the higher the resolution, the closer the simulation results approach the underlying biology. For extracellular electrical stimulation, the membrane potential of each section is calculated using equation (2): where V e is the extracellular potential. The extracellular potential arising from electrical stimulation using, for example, a simple, monopolar disk-like electrode, is commonly computed in NEURON using equation (3): where I is the extracellular current delivered by the stimulating electrode, r s is the transfer resistance of the electrode, r is the radius of the disk, a and z are the radial and axial distance between the neuronal section and the center of the electrode respectively. One can readily observe that equation (3) implies a homogeneous extracellular resistivity because the extracellular potential V e simply scales inversely with the radial distance from the stimulating electrode. While easy to compute, such an approximation suggests that NEURON is incapable of representing a non-homogeneous extracellular environment. This is the primary reason why we and others use COMSOL, a finite element modeling (FEM) solver, to calculate a more accurate version of V e (also see equation (4)) for the non-homogeneity of the extracellular environment (Guo et al 2014, Zelechowski et al 2020.

Step 2: extracellular stimulation environment configuration
Our toolkit accepts and incorporates into the COMSOL extracellular environment not available in NEURON's .hoc files (and hence not captured by HH-based, cable-equation neuronal models; figure 2(b)). Users can customize various model parameters including the radius and length of the nerve and fascicles, the resistivity of the interfascicular epineurium and endoneurium, as well as the geometry and electrical properties of the electrode(s). Unless noted otherwise, the conductivity of the endoneurium was set to 0.0571 S m −1 , and the conductivity of the internal epineurium was set to 0.01 S m −1 (Pelot et al 2018b). The geometries of the TIME used in this study is described in section 2.4 of Methods.

The toolkit engine
After having specified an initial, user-provided model in NEURON, our toolkit took the foregoing NEURON nerve model as input, and automatically recreated a COMSOL-based nerve model with added extracellular space (fascicles, nerve bundle) and electrode(s) for finite element analysis (figure 2(c)).

2.3.1.
Step 3: auto-generate COMSOL model from NEURON model Specifically, NEURON's .hoc file allows users to access individual axonal sections, and configure their geometry, morphology, and electrical properties. A multi-compartment model, such as a myelinated fiber with nodes of Ranvier and segments covered by Schwann cells, can be built by stringing together these sections. On the other hand, in COMSOL, models are generally created using the graphical user interface (GUI). The GUI provides users the freedom to design model's shapes and physics. For example, one can use small-diameter cylinders to represent the node of Ranvier and large-diameter cylinders to represent the portions with Schwann cells. Each cylinder is then assigned with the corresponding material properties (resistivity in particular). While simplifying usage, this precludes automation. As an advanced feature, COMSOL provides a MATLAB (also available for Visual Basic and Java) interface, enabling programmatic creation of models. Since NEURON (from version 7.0 onwards) comes with Python language support, our toolkit uses a Python program to orchestrate the creation of a COMSOL model from the existing NEURON model. In our sciatic nerve example, the toolkit extracted the 3D position, radius, height, and conductivity of each nodal and internodal section from the .hoc file. These pieces of information were automatically written as a MATLAB script that described the equivalent model in COMSOL, comprised of cylindrical elements with the same properties as those in the NEURON model. Subsequently after executing the MATLAB script, the equivalent model was re-created within COMSOL.

2.3.2.
Step 4: auto-mesh COMSOL model and calculate electric potential The generated COMSOL geometry from step 3 is automatically meshed into tetrahedral elements. The COMSOL meshing grid resolution is configurable.
High resolution ensures simulation accuracy but tends to consume more computational resources (e.g. simulation time and memory, see also figures 5(c) and (f)). In our nerve model, the meshing resolution was set to be 'finer' (average element size of 21.13 µm) for each fascicle and 'extremely fine' (average element size of 8.06 µm) for the A, C fibers due to their extremely large radius-to-length aspect-ratio. The average element size on the electrode boundaries where stimuli were applied was 8.47 µm. The extracellular potentials were computed by solving Laplace's equation for each meshed element using the conjugate gradients iterative solver. We applied the default Neumann's continuity condition on interior boundaries. Furthermore, we made the simplification of not considering any reactive components (especially, capacitive components) in the neural tissue, thereby assuming only ohmic (steady state) behavior following stimulus delivery. This allowed quasistatic approximation of Maxwell's equations within the nerve volume (Raspopovic et al 2017).
Laplace's formula can thus be used to solve for the voltage by viewing it as an electromagnetic question (equation (4)): where ρ n is the resistivity of the nodal elements, and V e is the extracellular voltage. The computed voltage provided a direct visualization of how fibers experienced the extracellularly applied stimuli. Figure 2(c) shows that when moving further away from the electrode, the transversal voltage decayed due to the resistance of the epineurium and endoneurium.

2.3.3.
Step 5: auto-export extracellular voltage After obtaining the voltage generated by the electrical stimulus over the non-homogeneous tissue media using COMSOL, the toolkit automatically exports the resulting data as a text file (figure 2(d)). The extracellular voltage over three-dimensional space is saved as a matrix of values. This text file is used by the toolkit again in the next step. In addition to the voltage text file, the engine also outputs the COMSOL morphology file and MATLAB for user to view and run the model in COMSOL.

Step 6: auto-import and interpolate extracellular voltage
The voltage matrix computed by COMSOL is automatically imported back to NEURON by reading the spatial coordinates and corresponding voltage values. The voltages computed by COMSOL may not correspond to the exact locations at which NEURON solves the HH-based cable equations. Thus, these data points are interpolated to the corresponding NEURON sections using Python's LinearNDInterpolator() function from the scipy module by our toolkit. Leveraging NEURON's extracellular mechanisms (the xtra.mod file and stim.hoc file), the voltage values were converted into transfer resistances which coupled the membrane to the extracellular stimulating current source (i.e. the electrode) (figure 2(e)). We rely on the principle of reciprocity, which assumes that the intervening tissue can be treated as linear (i.e. ohmic, as we did in the auto-generated COMSOL model of step 3). Suppose a stimulus current of amplitude I s , applied to a particular configuration of extracellular electrode(s), produces a potential V ext (x, y, z) at location x,y,z, then the transfer resistance r xtra between the electrode(s) and location x,y,z is simply (equation (5)): For an electrode placed at a particular location delivering current I s , equation (5) converts the COMSOL-computed results (V ext ) into a form (r xtra ) used by NEURON when solving each fiber response to the extracellular stimuli.
As a side note, because of the one-way coupling of physics between COMSOL and NEURON, it is possible to avoid re-simulating the COMSOL extracellular voltage values even when the stimulation current amplitude changes. Due to the ohmic-representation and reciprocity principle (equation (5)), any changes to the extracellular voltage V ext when the stimulation current I s changes (from the original COMSOLsimulated current amplitude) can be easily calculated by scaling the transfer resistance r xtra by a matching scalar factor.

Simulating fiber responses of the NEURON model using the interpolated voltages
Once the extracellular voltage values are imported and interpolated into NEURON as transfer resistance, they are used to convert the user-specified stimulus current waveform into extracellular stimulation voltage values. A stimulating electrode delivering current I produces a voltage V e (equation (6)): It should be noted that computing the extracellular voltage matrix in COMSOL (step 4) is by far the most computationally involved step (see also figure 5). A key advantage of our automated hybrid approach is that these voltage values only need to be computed once, via a time-invariant stimulation in COMSOL (e.g. 1 µA in our nerve model). The COMSOL-computed extracellular voltage are embedded into NEURON via equation (1) and are reusable across varying stimulation current amplitudes by simple linear scaling, as noted previously in step 6. This allows users to vary the stimulation waveform amplitude without needing to re-construct and re-simulate the COMSOL extracellular environment again. Therefore, fiber responses can be solved much more quickly than in COMSOL, using cable equations. This remains true as long as the stimulation configuration is not changed.

Stimulation electrode geometry
We simulated two types of TIME arrangements: monopolar and hexapolar. The monopolar configuration had one disk stimulating electrode (100, 200, 400 and 800 µm diameters) with a large remote return electrode. The hexapolar configuration had one stimulating electrode surrounded by six return electrodes (80 µm electrode-to-electrode distance and 40 µm diameter). All the electrodes were recessed 5 µm deep into a polyimide substrate. The polyimide substrate was drawn from COMSOL's built-in material 'P25N Polyimide-based No Flow Prepreg' , with 1000 µm width, 4000 µm length, and 10 µm thickness (Lago et al 2007, Mueller et al 2017. The TIMEs were inserted from three different angles (horizontally, vertically, and with 45-degree; also see figure 12 in appendix B) into the nerve model.
The stimulating and return electrodes were all with 100 nm thickness, and 9.43 MS m −1 conductivity. The electrode positions were fixed. The stimulating electrode was set to emit 1 µA current in all configurations, according to equation (7): where Ω represents the electrode boundaries and n is the surface normal. It applies a constant voltage to the boundaries such that the total normal electric current density J equals the total current I 0 . Each return electrode in the hexapolar setup was set to sink one-sixth of the current injected by the stimulating electrode.

Extracellular electrical stimulating parameters
We tested fiber responses under different TIME diameters and configurations using biphasic pulses. The biphasic stimuli were cathodic-first pulses, with 1 ms pulse width and 0.5 ms inter-pulse interval. The pulse amplitude was increased from 20 to 45 µA for A fibers, and from 10 to 150 µA for C fibers due to different activation thresholds. We also simulated kilohertz sinusoidal waveforms (Patel and Butera 2014) with various amplitudes and frequencies. The stimulation amplitudes were 1-9 µA for A fiber, and 10-500 µA for C fiber due to different activation thresholds. The stimulation frequencies ranged from 1 to 9 kHz. The largest stimulation frequency (9 kHz) required a sampling frequency of at least 18 kHz during simulation according to the Nyquist theorem. In practice, we chose the NEURON simulation time-step to be 5 µs, corresponding to a sampling frequency of 200 kHz, to ensure sufficient time-resolution for the partial differential equation solver to converge correctly in all stimulation conditions. The kilohertz sinusoidal stimuli were applied for 100 ms.
When preferentially activating C fibers in the sciatic nerve, it is desirable to use the smallest stimulation amplitude and the lowest frequency capable of activating C while blocking A fibers, to minimize the stimulator's power consumption and to avoid electrode degradation. We chose a 4 kHz, 100 µA, and 100 ms duration stimulation waveform for all analyses in preferential excitation modelling. To better visualize the voltage produced by the TIME setup, we adjusted the voltage range displayed to cover solely 0-10 mV.
Additionally, we devised a new waveform using cathodic DC ramp combined with anodic kilohertz ramp as extracellular stimuli (more details in section 3.8 of Results). This waveform consists of three phases (more details in figure 9): a cathodic DC ramp phase that eliminates onset response, an anodic kilohertz frequency stimulation (KFS) ramp phase that brings the waveform back to zero DC offset, and the KFS phase that maintains the blockade of A fiber(s) afterwards. We gradually increased the cathodic DC ramp amplitude from 0 to 3 µA with an increment of 1 µA. The 0 µA cathodic DC ramp waveform was essentially KFS and was used as the control. The amplitude of KFS was set to be the blocking threshold of the A fiber (2.45 µA). The cathodic DC ramp time was 50 ms. The anodic KFS ramp time was increased from 0 to 250 ms with an increment of 10 ms.

Selectivity index (SI)
When simulating population-level fiber responses, to compare C fiber selectivity performance between the monopolar and hexapolar configurations, we used the SI as defined by equation (8): where A and C are the number of recruited myelinated and unmyelinated fibers within each configuration. SI C aims to maximize targeted C fiber activity while minimiz non-targeted A fiber activity.

A computation model of rat sciatic nerve implanted with TIME electrodes
We used our toolkit (figure 1(e)) to automate the construction of a sciatic nerve model implanted with a . It incorporated user-specified fascicles, nerve bundle and TIME configuration to generate the extracellular stimulation model for a sciatic nerve in COMSOL, which was used to compute the extracellular voltage along fibers during electrical stimulation through finite element analyses. The toolkit synchronized data transfer between NEURON and COMSOL, so that fiber responses to the extracellular stimuli could be efficiently computed in NEURON using cable equations. The implementation details and step-by-step demonstrations are available in sections 2.1-2.3 of Methods.

Validating the automatically created model in NEURON
For model validation, we began by comparing the action potentials (APs) elicited in our modeled fibers against published models following a 0.1 ms suprathreshold intracellular pulse injection. The MRG model exhibited a depolarizing afterpotential (DAP) of ∼15 ms and a hyperpolarizing afterpotential (AHP) of ∼80 ms (McIntyre et al 2002). As shown in figure 4(a), our A fiber response waveform closely matched the published data by producing a DAP of 16 ms and AHP of 82 ms. The small differences were due to the fiber diameter change (see step 1 in Methods). C fiber APs also closely matched published data (Sundt et al 2015) (figure 4(b)). We then verified the spike propagating conduction velocity (CV) of both fiber types. In figures 4(c) and (d), we applied the intracellular current pulse at the proximal node of both fiber types and recorded their responses proximally, medially, and distally. APs were observed at all locations following brief propagation delays, indicating successful spike propagation. The CV of the A fiber was 25.25 m s −1 and of the C fiber was 0.24 m s −1 , consistent with the propagation speed of these fiber types (Harper and Lawson 1985). The APs also retained their general shapes, as expected of regenerative spikes created by active ionic conductance. We also validated our toolkit for extracellular stimulation. We used our toolkit to recreate a published myelinated axon model, then simulated its responses to extracellular stimulation. We then compared the responses so generated against the responses produced by the original model (Reilly 2016) based on NEURON. The simulation electrode was positioned 2500 µm distally from the myelinated axon, consistent with the existing literature, and delivered a biphasic square wave with a phase duration of 5 µs. Figure 4(e) shows that the APs generated by our toolkit closely mimicked the shapes exhibited by the NEURON-generated APs described in the literature. The slight deviation can be attributed to improved representation of extracellular potentials using our finite-element approach in COMSOL. Specifically, the published model used a single equation to approximate extracellular resistance, while our toolkit captures tissue resistance variations in intricate details, across different tissue types (e.g. endoneurium verses epineurium), and at high spatial resolution.

Computational resources for a large-scale nerve modeling
Simulating complex biological structures such as the nerve in figure 3 requires a careful balance between computational resources, time, and accuracy. To provide some guidance on these considerations, we measured several usage-related issues for our toolkit: simulation time, RAM usage, and auto-generated file size. We examined how these factors are influenced by the number of fibers modeled, the size of extracellular environment (i.e. nerve volume), and the FEM mesh resolution. The FEM mesh resolution was quantified by averaging the sizes of individual meshed elements across all COMSOL domains.
Because our hybrid modeling approach uses both NEURON and COMSOL, the analyses were conducted separately for these packages. For the COMSOL portion, the simulation time and RAM usage refer to the time elapsed and physical memory used to complete steps 3-5 of Methods, which involve autogenerating the COMSOL model (in the form of a MATLAB script), automatic meshing, automatic solving the FEM model, and automatic export of the computed voltage data. The auto-generated file size refers to the size of the MATLAB script. For the NEURON portion, the simulation time and RAM usage refer to the time elapsed and physical memory used to complete step 6 of Methods, which includes importing and interpolating the voltage data (provided by COMSOL) in NEURON. The autogenerated file size refers to the size of text file describing the computed extracellular voltage profile.
For the COMSOL part, simulation time, memory requirement and file size scaled linearly with the number of fibers ( figure 5(a)). The volume of extracellular environment/nerve does not readily influence any of these variables ( figure 5(b)). Finally, mesh resolution had an exponential relationship with simulation time and RAM usage, while having no impact on the auto-generated MATLAB file size ( figure 5(c)).
For the NEURON part, simulation time and autogenerated file size scaled with the number of fibers in the model. But fiber count had relatively little impact on the RAM usage ( figure 5(d)). The volume of the extracellular environment/nerve had no influence on any of these variables ( figure 5(e)). The FEM mesh resolution was found to have an exponential effect on NEURON's simulation time and the auto-generated text file size, while having limited impact on the RAM usage (figure 5(f)). For our nerve model, with 1,170 fibers in total, 139 mm 3 model volume, and 'extra fine' mesh resolution, the total COMSOL simulation time was around 4 h. The RAM usage was 18.3 GB. The autogenerated MATLAB file size that described the model was 24.7 MB. The auto-generated voltage data file size was 95.7 MB. While automatic generation of a fiber with uniform shape and conductivity (e.g. C fiber) can be handled using 'for loops' in MATLAB, our toolkit avoids 'for-loop' usage to consider nonuniform fiber shape and conductivity (e.g. A fiber), at the expense of increasing the auto-generated file size. Importantly, 1 kB of MATLAB file size corresponds to roughly 20 lines of code. Hence, the 24.7 MB file corresponds to approximately 494 000 lines of MATLAB code that were automatically generated by the toolkit (figure 13, appendix C). This demonstrates how the toolkit saves significant amounts of time and effort in constructing the nerve model. It is also crucial to observe that the NEURON portion of the simulation is completed significantly faster than the COMSOL portion (e.g. 400 s vs 4 s for 100 fibers), due to the simplicity of cable equations. Indeed, this is a key advantage of the hybrid modeling approach. Furthermore, once the COMSOL voltage field is computed for a particular combination of nerve information and electrode placement, one can run the NEURON-side analyses without having to repeat the COMSOL simulation again (see Methods sections 2.3.4 and 2.3.5), thereby saving considerable time when exploring the effects of various neuromodulation approaches.

The toolkit can be used with a large variety of nerve geometries and electrode configurations
While we focus on the sciatic nerve implanted with monopolar or hexapolar TIME electrodes in this paper, the toolkit is highly adaptable. It can be applied to any complex, cable-like neuronal structures, such as the three nerve arrangements in figure 6. It can also incorporate various user-defined TIME electrode arrangements, such as the four demonstration electrode configurations in figure 6. As shown in the appendix C, creating new nerve models and changing the stimulation electrode configuration are easily achieved through only a small amount of programming efforts.
3.5. TIME diameter affects spatial confinement of voltage and fiber responses in the nerve TIME was reported to be able to elicit more spatially focused fiber activation compared to conventional settings such as the cuff electrodes, the flat interface nerve electrodes (FINE), and the longitudinal intrafascicular electrodes (LIFE) (Badia et al 2011a). However, no experiment/clinical evidence has yet provided a comprehensive mapping about how different TIME electrode sizes (diameters) influence population fiber responses, at the resolution of individual fibers.
We simulated the voltage generated by different electrode diameters ( figure 7(a)). It was found that the electrode with the smallest (100 µm) diameter produced the sharpest voltage peaking at around 1.2 mV and with the fastest voltage amplitude decay over space. In contrast, the electrode with the largest (800 µm) diameter caused the flattest voltage peak, at around 0.2 mV, and with the slowest voltage amplitude decay over space. These results suggest that smaller electrodes tend to induce more spatially confined voltages than the larger electrodes.
Next, we examined how such differences in voltage confinement across electrode sizes affect fiber responses. Charge densities of each electrode under all current amplitudes are shown in figure 7(b). For small stimulus amplitudes (<32 µA), smaller electrodes recruited similar number of, or very slightly more, A fibers (figure 7(c)) comparing to the larger electrodes. But for large stimulus amplitudes (>32 µA), larger electrodes recruited more A fibers than smaller electrodes.
Repeating the analyses for C-fiber recruitment, as a function of electrode diameter and current amplitude, we observed a similar, but more pronounced, trend (figure 7(d)), that smaller electrodes recruited comparatively more fibers than the larger electrodes when the amplitude is small (<90 µA). However, this trend reversed as the stimulation amplitude increases further. This held true across all three TIME insertion angles ( figure 12, appendix B). In summary, our simulation suggested that smaller (i.e. ⩽ 400 µm in our study) electrodes can recruit less fibers given the stimulus amplitude is above a certain threshold compared to large electrodes, and thus provide more spatially focused fiber activations.

TIME arrangement affects spatial confinement of voltage and fiber responses in the nerve
Electrode configuration also affected spatial confinement. The hexapolar configuration (figure 8(a) bottom) was proposed to create a more confined electric field in retinal stimulation (Habib et al 2013, Matteucci et al 2013, Italiano et al 2022. However, it remains unclear whether hexapolar configuration also offers more focal neural recruitment in peripheral nerve stimulation. We found that the hexapolar setup induced an extracellular voltage peak at a smaller amplitude (1.3 mV) compared to the monopolar setup (1.8 mV) ( figure 8(b)). The extracellular voltage near the return electrodes was sunk to a negative value at around −0.1 mV, which contained the extracellular voltage profile nearby. Comparing A and C fiber recruitment between these two setups (figures 8(c) and (d)), hexapolar setup was found to activate less A and C fibers consistently under all stimulus amplitudes. Therefore, we concluded that the hexapolar configuration can provide a more spatially confined extracellular voltage profile as well as fiber responses.

Kilohertz sinusoidal stimuli for preferential C fiber recruitment
Ideally neuroprostheses should provide spatiotemporal selectivity as well as fiber-type selectivity. A major confounding issue of fiber-type selectivity in peripheral nerve stimulation is that A fibers have considerably lower activation threshold than C fibers during extracellular stimulation (Parker et al 2017). It is generally difficult to recruit C fibers selectively without undesirably co-activating large numbers of neighboring A fibers. A possible solution for activating C fibers with minimal A fiber responses is through conduction blockade of A fibers using KFS waveforms (Patel and Butera 2014, Pelot and Grill 2020, Neudorfer et al 2021. This method uses finely tuned frequency and amplitude combinations to inactivate the A fibers (Chang et al 2022), thereby preventing generation of further APs. However, the effectiveness of this strategy across a complete nerve, with structural inhomogeneity throughout, has not yet been investigated. Herein, we study the combined effects of (1) a µm-scale TIME placed near the target fibers, and (2) delivering kilohertz stimuli for improved C fiber selectivity.
The A and C fibers responded differentially for a given sinusoid stimulus ( figure 9(a) and (b)). These observations are consistent with studies in the vagus nerve (Chang et al 2022). As such, A and C fibers demonstrated disparate parameter spaces of excitation and blocking, offering the possibility to recruit each type preferentially through proper choices of stimulation parameters.
Next, to achieve preferential excitation/blocking, we used stimulation frequencies and amplitudes that maximally differentiate between the A and C fibers in the following analyses (figures 9(c) and (d)). At a 4 kHz stimulus, the A fiber started firing at  1.5 µA. Its spike rate increased as the current amplitude increased and reached it maximum at 25 µA. However, as the stimulus strength increased further (>25 µA), the spike rate decreased substantially and eventually ceased firing after an initial burst of spikes (i.e. onset responses) at ⩾ 40 µA, creating a nonmonotonic response profile as a function of stimulus amplitude. Conversely, the C fiber started firing at a much higher current (170 µA). Its spike rate increased as the current amplitude increased.
As shown in figure 9(d), by modulating the A fiber into the inactivated state using supra-threshold current amplitudes, it becomes possible to recruit the C fiber selectively, while simultaneously avoiding further A fiber responses. But it should be noted that kilohertz stimulation causes the A fiber to exhibit an initial onset burst response, which is a brief period of repetitive firing (Bowman and McNeal 1986, Kilgore and Bhadra 2004. Depending on neural interfacing applications, this may not be desirable. A solution to this problem is presented in the next section.

Cathodic DC ramp combined with KFS eliminated A fibers' onset responses
Our results from figure 9(d) agrees with previous studies-before nerve conduction blockade occurs with kilohertz stimulation, there always exists an undesirable onset response (Patel and Butera 2014). It has been verified through both simulations and in vivo rat sciatic nerve experiments that an anodic amplitude ramp was ineffective in eliminating these onset responses (Miles et al 2007). However, a cathodic ramp was proposed to be effective in eliminating onset response (Ackermann et al 2011). But it has the drawback of needing two electrodes and two stimulators in the setup, one for delivering the cathodic DC ramp, and the other for delivering the KFS waveform. For medical neurostimulators, this would increase circuit complexity and device size.
To overcome these shortcomings, we proposed a new one-electrode strategy incorporating cathodic DC ramp (figure 10(a)) to eliminate onset response in kilohertz-stimuli-induced conduction block. With an increasing cathodic DC ramp amplitude, we found that A fiber's onset responses were gradually reduced and eventually eliminated ( figure 10(b)). However, too large of a DC ramp amplitude (>3 µA) without also increasing the anodic KFS ramp time would induce offset responses-spikes occurred during the second, anodic KFS ramp phase. Such offset responses could be eliminated by increasing the anodic KFS ramp time from 200 ms to 240 ms (figures 10(c)-(d) summarizes the effects of cathodic DC ramp amplitude and anodic KFS ramp time on eliminating A fiber's onset and offset responses). In summary, a sufficiently large cathodic DC ramp amplitude paired with a long anodic KFS ramp time could eliminate both onset and offset responses during KFS-induced conduction blockade.

Preferential excitation of C fibers inside the whole nerve
In figures 9(c) and (d), by choosing the stimulation current amplitude and frequency appropriately, it is possible to preferentially activate A and C fibers adjacent to the stimulation electrode. However, it is not clear how fibers further away from the electrode would behave. In particular, C fiber activation needs high current amplitudes. This in turn creates a large extracellular voltage field, possibly causing unintended fiber activation distant to the electrode. To assess such possibilities, we examined the responses of all fibers in the nerve during electrical stimulation.
As shown in figures 11(a) and (b), the extracellular voltage gradually decayed with distance. Moreover, the electrode array's polyimide substrate acted as an insulating barrier, concentrating most of the extracellular voltage along only half of the nerve. The slight 'spill-over' into the opposite half was due to the current flowing over the edge of the polyimide, on the bottom side of the substrate. This 'hemisphere-like' shape of voltage resulted in most of the fibers being activated or blocked in one half of the nerve. Figure 11. Preferential excitation of C fibers while blocking A fibers; PE: preferential excitation. (a) Due to the decaying voltage distribution produced by monopolar setup, activating C fibers causes undesirable co-activation of A fibers. (b) In hexapolar setup, the induced voltage was comparatively more confined, resulting in less mis-activated A fibers. (c) The average percentages of C fiber activation, A fiber blockade, and A fiber mis-activation for monopolar and hexapolar TIME across three TIME electrode insertion angles. (d) The average selectivity index for C fibers for monopolar and hexapolar TIME across three electrode insertion angles.
The electrical stimuli activated several C fibers next to the electrode. Moving away from the electrode, the extracellular voltage weakened and could no longer activate C fibers. At these distances, A fibers could still experience depolarization block, thereby ensuring suppression of A fiber responses. However, further away, the extracellular voltage became sufficiently low such that A fiber blockade was no longer possible.
Because hexapolar electrodes gave rise to more spatially confined electric field ( figure 8(b)), we asked whether mis-activation of A fibers can be reduced using such an electrode configuration ( figure 11(b)). Quantitatively comparing the two configurations, the monopolar setup could preferentially activate 0.51% C fibers while blocking 2.39% A fibers, but misactivating 7.01% A fibers on average among the 1,170 fibers in the nerve across three TIME electrode insertion angles. The hexapolar setup, by creating a more confined electric field, recruited fewer C fibers (0.48%), blocked fewer A fibers (1.34%) near the electrode, but also mis-activated less A fibers (6.01%) at the nerve periphery ( figure 11(c)). The hexapolar electrode setup improved spatial localization by reducing the number of unintended A fibers activation. Summary plots for the horizontal and vertical TIME electrode insertion setup can be found in the appendix D.
We used the SI (see Method section 2.6) to quantitatively assess the differences between monopolar electrode and hexapolar electrode in activating C fibers. The average SI produced by monopolar and hexapolar setups were 0.073 and 0.08 respectively ( figure 11(d)), suggesting moderately improved selectivity using the hexapolar electrode. However, the main benefit of using the hexapolar arrangement is improved spatial focality of the electrical stimuli, as described previously.
In summary, at the scale of the whole nerve, preferential C fiber recruitment is generally only achievable in a localized area using monopolar TIME electrodes. Some mis-activation of distant A fibers are often unavoidable. However, the extent of A fiber misactivation can be reduced by using a hexapolar TIME.

Discussion
Several in vivo animal investigations and pilot clinical studies have been done for TIME arrays (reviewed in Raspopovic et al 2021). However, quantitative assessments of TIME-evoked responses across an entire sciatic nerve, particularly in the context of unconventional stimulation waveforms, remains unexplored. The sciatic nerve is structurally complex ( figure 1(a)), with densely packed sub-microndiameter axons (Giuffre and Jeanmonod 2022). Some axons are myelinated while others are not. The axons are bundled into fascicles. These anatomical intricacies preclude large-scale, high-resolution electrophysiological investigations involving thousands of axons simultaneously and at the resolution of individual axons. Computational models of nerve stimulation are indispensable tools for the development of better neural interfaces (Raspopovic et al 2011). There are several modeling approaches that are commonly used. Simulation environments such as NEURON (Hines and Carnevale 1997), based on the HH formalism (Hodgkin and Huxley 1952) and cable equations, can efficiently model neuronal/axonal electrical responses. On the other hand, FEM software such as COMSOL Multiphysics (COMSOL, Inc. Palo Alto, CA, USA) are used to capture the effects of a non-homogeneous tissue environment due to myelination and intrafascicular space under extracellular stimulation. However, the latter method is computationally intensive. As such, a new class of 'hybrid' methods have emerged (Zelechowski et al 2020) to leverage the advantages of each approach. These models offer the computational efficiency of cable equations while providing accurate extracellular space representation through FEM. But with few exceptions (Eiber et al 2021, Musselman et al 2021, existing hybrid models are tailor-made for specific nerves and particular electrodes. They also require manual transfer of results between disparate simulators. These shortcomings limit model scalability and generalizability to other nerve types, requiring significant computational biology skills and considerable operator effort.
We developed an automated toolkit to simplify the creation and usage of complex nerve models. Using this toolkit to build a peripheral nerve model, we showed that it was generally easier to achieve focal activation with small TIME electrodes when the stimulus amplitude is large. We found that hexapolar electrodes further improved the spatial confinement of fiber activation. We also explored the current ranges and frequency ranges required for preferential activation of A and C fibers using kilohertz stimulation. Second, we proposed a new stimulation waveform that could eliminate undesirable onset responses during A fiber blockade by only using one stimulation electrode. Third, our population-based simulation suggested that the fiber-type selectivity offered by kilohertz stimuli is only achieved near the stimulation electrode. Nevertheless, fiber-type selectivity can be further improved by using hexapolar configuration.

Small TIME electrodes may not be necessary in practice
Our model suggests that while TIME electrode with smaller diameter produced a sharper extracellular voltage profile, they do not always provide more focal neural recruitment than the larger electrodes ( figure  7). Indeed, when the stimulation current amplitudes are low, small electrodes can at times recruit marginally more fibers than the larger electrodes. The benefit of small electrodes for confined activation only became consistently so when the current amplitudes are large. Therefore, if only small current levels are used, as might typically be the case to reduce stimulator power consumption and avoid electrode degradation, the added complexity of manufacturing small electrodes might not be warranted.

Hexapolar TIME electrodes improve spatial confinement and fiber response confinement
We also encourage using hexapolar TIME configuration to improve spatial localization of fiber responses. Our simulation (figure 8) showed that hexapolar setup produced better voltage profile and fiber responses confinement compared to the monopolar setup. Such improvement is consistent throughout all stimulus amplitudes. This agrees with previous findings that hexapolar electrode configuration possesses a more confined current distribution than the monopolar setup in nerve stimulation (Jolly et al 1996, Matteucci et al 2013, Barriga-Rivera et al 2017. The current emitted by the center stimulating electrode was sunk by the surrounding return electrodes, preventing dispersion of electric field to the nerve periphery (Joarder et al 2011).

Preferential excitation of C fibers through kilohertz stimulation is effective only in a local region
At the single fiber level, we examined the selective stimulation between A and C fiber. We verified that with sufficient current amplitudes, large fibers (A) exhibited depolarization block before the excitation of small fibers (C) (figure 9), agreeing with recent findings reported in rat vagus nerve (Chang et al 2022). However, in a complete nerve bundle with a large fiber population, the level of selectivity reduced due to the non-uniformly distributed fiber locations relative to the electrode placement (figure 11), confirming the challenge proposed by a previous study (Patel and Butera 2018). In particular, A fibers near the electrode (<800 µm) could be blocked while those further away may not be blocked, due to the decaying voltage amplitude over space. In general, although electrode setup can be arranged to optimize spatial activity, no configuration can achieve true selectivity.

Cathodic DC ramp combined with KFS can block A fiber without inducing onset responses
Significant efforts have focused on developing stimulation techniques to minimize side effects caused by onset responses (Petruska et al 1998, Miles et al 2007, Ackermann et al 2011. We developed a new waveform that is capable of eliminating A fiber onset responses, which normally occur prior to conduction block. This waveform has the advantages of requiring only one electrode, compared to the two-electrode setup (Ackermann et al 2011). While it reduces invasiveness, we also envision such a waveform to be useful in many other scenarios. For instance, one can enforce unidirectionality of extracellular stimulation by silencing fiber activity in the unwanted direction. Also, one can inhibit fibers completely, rather than partially, to aid disease treatment. The next step would be to verify this technique at a population level, for example, in a complex nerve. Our toolkit provides researchers the opportunity, in an accurate and efficient manner, to observe population-based neural responses under newly designed devices that aim to investigate this issue.

Automated NEURON-COMSOL interface allows high-speed, high-accuracy nerve modeling
Computational models of nerve electrical stimulation can be significantly improved by accurate FEMbased simulation of the extracellular environment while maintaining the computing efficiency of cableequation-based simulation (Frankemolle et al 2010, Raspopovic et al 2017, Pelot et al 2018a, Ge et al 2020, Zelechowski et al 2020. This study presents an efficient and accurate nerve modeling toolkit for the NEURON/COMSOL environment. Conventional hybrid approaches require the same nerve to be built twice in the two environments separately, followed by manual integration of simulation results. This burden is prohibitive for large-scale studies, involving complex nerves (Miocinovic et al 2006, Schiefer et al 2008, Song et al 2020. As a solution, our toolkit automatically wrote out the code for the modelled nerve. Furthermore, by automatically exchanging extracellular voltage data between the two platforms, our toolkit frees users from doing repetitive work and significantly speeds up the simulation process. Recent nerve modeling has given more attention to automated simulation pipeline to efficiently exchange information between finite element analysis and NEURON environments. Compared to a recently developed automatic computational tool of visceral nerve ViNERS (Eiber et al 2021) using MATLAB and relying on additional external dependencies GMSH (an open source 3D FEM generator) (Geuzaine and Remacle 2009) and EIDORS (Electrical Impedance Tomography and Diffuse Optical Tomography Reconstruction Software) (Vauhkonen et al 2019), followed by fiber simulation in NEURON, most practitioners will likely find our approach considerably more familiar and adaptable, due to the dominance of NEURON and COMSOL environments. This significantly reduces the possibility of algorithmic bugs and improves usability. Indeed, we suspect that many potential users will find our codebase to be easily integrated into their existing in-house modeling work. However, comparing to open-source FEM software, COMSOL being a commercial product might be cost-prohibitive for some users. More recently, a COMSOL-NEURON pipeline ASCENT (supported by Python and Java) was developed to characterize population-based nerve excitation and blocking thresholds with different clinical cuff electrode (Musselman et al 2021). Our toolkit, on the other hand, focused on testing µm-scale devices User-defined 1 T-NEURO module 2 Membrane potential 3 Single unit field potential 4 Evoked compound action potential (e.g. TIME) implanted within the nerve. Therefore, our toolkit can be seen as complimentary to ASCENT, by enabling studies using non-cuff-like electrodes.
Our toolkit further extends this ability by allowing definition of new electrodes beyond cuff and TIME electrodes. The other major differences between our toolkit and ASCENT are the granularity of anatomical representation and the computational speed. Comparing to ASCENT, which requires 98 s to simulate 100 MRG fibers when operating on a 16-core CPU, our toolkit takes 400 s. The four-fold time difference is due to the anatomical details that we consider, such as the Schwann cells and the unmyelinated internodes along each A fiber in a fascicle. ASCENT, on the other hand, treats each fascicle as homogeneous tissue, without considering the constituting fibers and possible anatomical differences between fiber types. Another popular commercial platform SIM4LIFE (https://zmt.swiss/sim4life/) (Zurich, Switzerland) also allows users to simulate arbitrary electrodes, fibers, nerves, and even nerve trajectories. While such flexibility supports complex and accurate neuronal simulation, it is not intended for scripted automation of, for instance, incorporating published models available from ModelDB repositories. A detailed comparison between similar tools is presented in table 1.

Limitations
A limitation of our toolkit is that it only incorporate neuronal properties captured in NEURON's .hoc object, such as sectional length, diameter, and resistivity. Biological structures like blood vessels are not represented. These structures may have significant presence in larger nerves, contributing nontrivially to extracellular potential during electrical stimulation. Similarly, our toolkit's ability to accurately represent the structure of nerves in the longitudinal direction is limited when using a simple extrusion method from a cross-section. In future updates, the toolkit can leverage NEURON's capability to describe axons with varying segment diameters using multiple point3d() statements and use these additional information to create fibers with changing diameter and/or trajectory over their path. Likewise, the COMSOLequivalent model can be built by dividing it into small segments and modifying the diameter, location, and trajectory step by step, thereby improving longitudinal accuracy.
Finally, our population-based model did not consider the full diversity of nerve fiber subtypes mediating different functions (Smith et al 2012, Rich andBrown 2018). Nevertheless, our toolkit can readily accommodate additional fiber types once verified electrophysiological data and morphologies become available, for example on ModelDB, allowing users to integrate them into the model.

Conclusion
In summary, we developed a toolkit that can automatically create the equivalent nerve fibers in COMSOL using nerve structures implemented in NEURON, and export the extracellular voltage computed in COMSOL back to NEURON for accurate and efficient simulation of non-homogeneous nerve tissue. This toolkit is easy-to-use, saving users significant amounts of time and effort when designing large population-based peripheral nerve models. Using this toolkit, a generic peripheral nerve model containing A and C fibers was built. With micrometer-scale TIME array placed inside the nerve, our simulation showed that smaller electrode size and hexapolar configuration offered more prosthetic control over the number of fibers activated. At the population level, we found that KFS can selectively stimulation C fibers at least locally. We also devised a new waveform that combines cathodic DC ramp and KFS to eliminate undesirable onset responses during KFS-induced conduction blockade.

Data availability statement
The data that support the findings of this study are openly available at the following URL/DOI: https:// modeldb.science/267618. Some code has been written to specify fiber distributions in the Fascicle.hoc file. In this example, we place ten A fibers randomly inside a 400 µm fascicle with uniform distribution. After storing all the .hoc and .mod files in the NEURON_In folder, the engine will automatically read from these files when pipelining. If users want to model other fibers or distributions, simply program and replace with new .hoc and .mod files that describe those.

A1.2. COMSOL input files
The COMSOL input file contains a Python script specifying the parameters of the extracellular stimulation environment.
param.py involves variables read by the engine, such as nerve size, fascicle size, electrode configurations etc. In param.py, it has been provided a full list of explanation on what each variable means. Below are the parameters for this example.  Once we store all the parameters in the COMSOL_In folder, the engine will automatically read from these variables when pipelining. If users want to modify the stimulation environment, simply change the variables accordingly in param.py and re-run the engine.

A2.
Step 2: running the engine Next, we run the script autoToolkit.py to initiate the pipeline. While users can treat this as a black box, here is a brief description of what it does. The engine calls NEURON2COMSOL_auto_conv.py to read from the NEURON_In and COMSOL_In folders and convert our NEURON nerve model to an equivalent COMSOL nerve model. It simulates a TIME electrode to extracellularly stimulate the nerve, meshes and studies the COMSOL model. The engine calls COMSOL2NEURON_auto_conv.py and setrx.hoc to transfer the extracellular voltage from COMSOL back to NEURON as transfer resistance.
We can setup this engine via a three-step process. Here is how: (1) open a new Python script in the working directory, include neuron and autoToolkit as headers, (2) call our NEURON fiber model in the NEURON_In folder, and (3) call autoToolkit's pipeline() function.
J. Neural Eng. 20 (2023) 046032 Y Xie et al`À n example for pipelining a simple fascicle between NEURON and COMSOL using the automated toolkit`# Step 1: import headers for (i) NEURON's .hoc object and (ii)  Voltage.txt contains the extracellular voltage along fibers as a matrix of values (a snippet below):

A3.2. NEURON output files
The NEURON_Out folder contains the Rx.txt. This is the interpolated transfer resistance along fibers derived from Voltage.txt following Ohm's law. It is to be read by NEURON's Xtra.mod when simulating extracellular stimulation. Below is a snippet: This brings us to the end of using the toolkit. Users can proceed to simulate fiber responses in an extracellular stimulation environment in NEURON. This can be done by specifying a stimulation waveform and loading setrx.hoc in the engine folder. For instance, the following stim.hoc code loads the 10-fiber fascicle in this example, then inserts to each fiber the transfer resistance corresponding to the extracellular voltage inside the nerve, then simulates the neve's responses to the extracellular stimulus (e.g. a 4 kHz, 2 µA sinusoidal waveform) delivered using NEURON. .v(0.5)", 1, 1) run() Users can run this script by typing the following command in the terminal: nrngui stim.hoc It will generate the response of fiber#0 in the fascicle: Appendix D. Preferential excitation in horizontally/vertically-inserted TIME setup Figure 14. Preferential excitation of C fibers while blocking A fibers in horizontal and vertical TIME insertion setup; PE: preferential excitation. (a) A and C fibers population response. (b) The average percentages of C fiber activation, A fiber blockade, A fiber mis-activation, and the average selectivity index for C fibers for monopolar and hexapolar TIME across three TIME electrode insertion angles.