Interleaved Multi-Contact Peripheral Nerve Stimulation to Enhance Reproduction of Tactile Sensation: A Computational Modeling Study

Peripheral nerve stimulation (PNS) is an effective means to elicit sensation for rehabilitation of people with loss of a limb or limb function. While most current PNS paradigms deliver current through single electrode contacts to elicit each tactile percept, multi-contact extraneural electrodes offer the opportunity to deliver PNS with groups of contacts individually or simultaneously. Multi-contact PNS strategies could be advantageous in developing biomimetic PNS paradigms to recreate the natural neural activity during touch, because they may be able to selectively recruit multiple distinct neural populations. We used computational models and optimization approaches to develop a novel biomimetic PNS paradigm that uses interleaved multi-contact (IMC) PNS to approximate the critical neural coding properties underlying touch. The IMC paradigm combines field shaping, in which two contacts are active simultaneously, with pulse-by-pulse contact and parameter variations throughout the touch stimulus. We show in simulation that IMC PNS results in better neural code mimicry than single contact PNS created with the same optimization techniques, and that field steering via two-contact IMC PNS results in better neural code mimicry than one-contact IMC PNS. We also show that IMC PNS results in better neural code mimicry than existing PNS paradigms, including prior biomimetic PNS. Future clinical studies will determine if the IMC paradigm can improve the naturalness and usefulness of sensory feedback for those with neurological disorders.


I. INTRODUCTION
T ACTILE sensation allows a human to detect, discrimi- nate, and identify external stimuli and respond to stimuli appropriately [1], [2].Tactile sensation can be lost through disorders such as limb amputation, nerve injury, spinal cord injury, and stroke.Losing tactile sensation creates numerous functional consequences, including impaired performance in tasks that require focus or fine motor skills [3].To mitigate these problems, researchers are developing neuroprostheses that use peripheral nerve stimulation (PNS) to restore tactile sensation.PNS of residual somatosensory nerves in the arm can elicit "artificial" touch sensations in the hand.Current PNS paradigms typically provide one or more sensory percepts on the hand, where input from a force sensor directly scales a PNS parameter, such as pulse amplitude (PA), pulse width (PW), or pulse frequency, to influence the perceived intensity of the evoked percept [4], [5], [6], [7].However, many sensations produced by PNS are described as unnatural tingling or paresthesia, which can be bothersome to participants [8], [9].
One key issue with PNS that may contribute to these perceptual deficits is that the neural activity created by PNS largely does not resemble the complex, dynamic firing patterns in natural touch [10].Recent work has begun to develop biomimetic PNS paradigms that attempt to mimic aspects of the biological neural code of natural touch.These biomimetic paradigms have demonstrated preliminary success in improving sensation naturalness and enhancing object detection reaction times, at least for some participants [11], [12].However, clinical implementation of biomimetic paradigms can be challenging, as it can be difficult to select specific stimulation parameters to evoke important aspects of the natural neural code.Prior studies in both humans and non-human primates have shown that neuron type, neuron location, neural firing rate, and neural population size over time are four important aspects of the neural code underlying the perception of touch stimuli [13], [14], yet no existing biomimetic paradigm mimics all four of these coding properties [11], [12], [15], [16], [17], [18], [19], [20].In addition, computational models of neural activation from electrical stimulation can help to predict neural activation patterns resulting from stimulation and allow us to systematically select stimulation parameters to achieve a given pattern of neural activity [21], [22], [23], [24].
We hypothesize that developing interleaved multi-contact (IMC) biomimetic PNS paradigms will enhance our ability to activate natural tactile codes and ultimately provide useful and intuitive sensory feedback to users of sensory neuroprostheses.We used TouchSim, a model that determines the response of mechanoreceptive afferents to mechanical stimuli applied to the hand [10], to identify the neural code to mimic with PNS (Fig. 1).We used an electrical activation model that predicts activation of peripheral afferents in response to nerve stimulation [42], in conjunction with non-gradient-based optimization algorithms, to select PNS stimulation parameters.We then used this tool chain to create several PNS paradigms constructed under different touch criteria, and used neural response reproduction accuracy as a metric to compare the performance between paradigms.We demonstrate that interleaved multi-contact PNS approaches improve neural response reproduction accuracy over current biomimetic paradigms.The paradigms presented here can be implemented in future clinical trials of PNS to determine how the sensations produced are perceived and utilized in sensorimotor tasks.

II. METHODS
Our goal was to reproduce neural responses to natural touch stimuli using PNS.We first used TouchSim to generate a set of neural responses to touch stimuli applied to the index finger that we aimed to replicate with PNS (Fig. 1a-b).We then coupled optimization with a model of neural activation from electrical stimulation applied through a multi-contact Composite Flat Interface Nerve Electrode (C-FINE) (Fig. 1c).This process resulted in a time series of parameters for C-FINE stimulation, called a playlist, that approximately reproduced the original neural responses (see Fig. 2 for an overview of our approach).

A. TouchSim Neural Activation Profiles (TNAPs)
TouchSim was used to generate neural responses to natural touch stimuli applied to the index finger [10].Briefly, the TouchSim model places low-threshold mechanoreceptive afferents in the palmar surface of the hand based on typical human innervation densities, and then predicts how the neurons will fire in response to spatiotemporal mechanical stimuli presented to the skin.We used TouchSim to simulate 1704 afferents in the index finger: 562 slowly-adapting type I (SA), 947 rapidlyadapting type I (RA), and 195 rapidly-adapting type II (PC).
TouchSim data was generated for 1 s indentation stimuli on the index finger that consisted of symmetrical ramp up and ramp down periods, flanking a sustained period (Fig. 1a).Two data sets were generated: 1) stimuli with a fixed depth and a ramp duration that varied from 0.05 s to 0.5 s in increments of 0.05 s, and 2) stimuli with a fixed ramp duration and Fig. 2. Optimization-based biomimetic interleaved multi-contact (IMC) peripheral nerve stimulation for tactile sensation.Blue boxes: TouchSim generated neural activation patterns for a set of mechanical touch stimuli (indentations), which were then split into 1-ms touch neural activation profiles (TNAPs).The TouchSim reachable set is the collection of all unique TNAPs for all modeled touch stimuli.Green boxes: Particle swarm and pattern search optimization were used to select stimulation parameters to reproduce each TNAP.The optimization runs the electrical activation model to generate electrically-stimulated NAPs (ENAPs) from stimulation parameters, calculates the reproduction error, and iterates until a minimum reproduction error is reached.Purple boxes: A lookup table stores the association between each TNAP and the optimal stimulation parameters to approximate it.Lookup tables using a fixed contact for all TNAPs were created using each electrode contact (1C) (n=15), along with sets of two-contact (2C) pairs (n=21).Playlists were generated by ordering stimulation parameters from the lookup tables in series to create a sequence of ENAPs that approximated a given touch stimulus.Playlists were generated for each 1C and the 2C non-IMC paradigms for illustration purposes, though these were not used for further analyses.Orange boxes: Error was compared across non-IMC lookup table results, and the contact, or pair of contacts, resulting in the lowest error for a given TNAP was selected for the 1C and 2C IMC paradigms, respectively.Playlists were generated for the 1C and 2C IMC paradigms by ordering stimulation parameters in series to create a sequence of ENAPs that approximated a given touch stimulus.Rectangles represent process steps, and rounded rectangles represent outputs.
a depth that varied from 1 mm to 6 mm in increments of 0.5 mm.We chose this range of ramp durations to span the possible values of a 1 s stimulus with symmetrical ramp on and ramp off periods.We chose this range of depths based on the values over which TouchSim was validated [10].In total, we simulated 20 touch stimuli -11 with varying depths and 10 with varying ramp durations (one stimulus was shared between the two sets).
TouchSim output the spike times for each modeled neuron (Fig. 1b).The neural response to each stimulus was parsed into 1 ms bins, for a total of 1000 bins spanning the 1 s stimuli.Each bin contains a neural activation profile (NAP), or a row vector with binary values representing whether activation occurred for each of the 1704 modeled index finger neurons (Fig. 1d).Each NAP is considered to be independent from other NAPs, such that no dependencies exist between temporally adjacent NAPs.However, a series of NAPs can be ordered in a sequence called a playlist to reproduce the neural response to a specific stimulus over time.
The unique TouchSim-generated NAPs (TNAPs) among all 20 touch stimuli were found, creating the TouchSim reachable set, or the set of NAPs we must be able to approximate with PNS to recreate the set of touch stimuli (Fig. 2).

B. Electrical Neural Activation Profiles (ENAPs)
To determine the optimal stimulation parameters for PNS, our optimization algorithm used a biophysical electrical activation model to determine which neurons in the median nerve were activated from each PNS pulse.PNS was delivered to the median nerve through a 16-contact C-FINE (Fig. 1c).Of these contacts, one was a designated return path and the other 15 contacts were used for cathodic stimulation.
The electrical activation model consisted of a previouslydeveloped finite element model (FEM) of a human median nerve [22] and a linear approximation method to determine neural firing.The FEM was created in ANSYS Maxwell 3D using ultrasound images of a median nerve that were taken during C-FINE implantation surgery for a previous clinical study [6], [19].The electrical and geometrical properties of the multi-contact C-FINE were simulated in the FEM [27], and voltage fields within the cuff-nerve complex were generated based on the propagation of a current through the tissue from C-FINE stimulation [18].All stimulation pulses were modeled as square cathodal pulses.Neurons were randomly placed within the modeled fascicles based on typical human fiber densities and with fiber diameters selected from normal human distributions, yielding a total of 7334 neurons in the median nerve [43], [44].
A linear approximation method was then used to determine neural activation of these axons in response to PNS [42].Briefly, the linear approximation method is an algebraic function that predicts activation of a myelinated neuron given inputs of stimulation PW, fiber diameter, and extracellular voltage at each node of Ranvier.The approximation was previously generated [42] by fitting a curve to the boundary between neural activation and non-activation produced by hundreds of iterations of the MRG active neural model simulated in NEURON [45].
The electrical activation model determined neural firing due to stimulation pulses delivered through the C-FINE.Each stimulation pulse was defined by the stimulating contact(s), PW, and PA.For a given stimulation pulse, the PW was held constant across all active contacts but the PA could vary across the 15 contacts.The electrical activation model output a vector with binary values representing whether each neuron was activated by the pulse of PNS.
Of the 7334 neurons in the electrical activation model of the median nerve, we selected 1704 neurons grouped within three fascicles to be in the index finger in order to match the number of neurons in the TouchSim model [10].The fascicles selected to contain index finger afferents were assumed to be adjacent Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
to one another based on prior studies [46], [47].Any neuron in the median nerve model that was not assigned to be in the index finger was removed.We then paired each remaining neuron in the nerve model with a neuron in the TouchSim model of the index finger.To produce this nerve-finger mapping, we considered the distributional characteristics of neuron type reported in the anatomical literature [48], [49], [50] and based the location assignment on an assumed somatotopy of the fascicles, such that neurons closer together in the finger were grouped together in fascicles.Since neurons in the TouchSim model had a specific position in the finger and an afferent type (SA, RA, PC), we assumed this same type and position identity for the corresponding neuron in the nerve model.Thus, when the electrical activation model was run, it produced an electrically-activated NAP (ENAP) of neurons in the index finger with the same dimensions as a TNAP, and each neuron in the ENAP had a neuron type, index finger location, and fascicle location (Fig. 1e).

C. Optimization
We used optimization to minimize the reproduction error between the TNAPs in the TouchSim reachable set and the ENAPs created with PNS.We defined the reproduction error between a TNAP and an ENAP, or the objective function of the optimization problem, to be the mismatch of four critical elements of the neural code of tactile sensation: neuron type, neuron location, firing rate, and population size.We calculated this mismatch as follows: First, each neuron modeled by TouchSim and the electrical activation model had a type (SA, RA, PC) and location in the index finger (proximal, middle, distal segment).Thus, there were a total of nine possible type-location combinations, and each neuron fell into one of these nine categories.Each NAP could then be represented as a 9-dimensional vector, where each component contained the count of active neurons in each of the 9 type-location categories.The firing rate was implicit because of the timing of the TNAP within the stimulus, and the population size was represented by the overall error vectors.
To quantify the error (E) between a TNAP and an ENAP, we took the absolute value of the difference between the 9-dimensional vectors generated from the TNAP (T ) and that generated form the ENAP (M), and then summed the components, as in Equation ( 1).
We used a hybrid of particle swarm and pattern search optimization algorithms to determine optimal stimulation parameters that are used to approximate the TNAPs via PNS (Fig. 2).For each TNAP being approximated, particle swarm optimization was run first, and then pattern search optimization was run using the particle swarm optimization results.We limited PW to 0-0.2 ms and PA to 0-0.5 mA.Our optimization results were verified for a subset of TNAPs with an enumeration method.See [51] for more information on the optimization parameters and enumeration method.
We then repeated this process for each TNAP in the reachable set.Although some TNAPs might appear multiple times within a given stimulus, or multiple times across all modeled stimuli, we only ran the optimization once per TNAP in the reachable set to reduce total computation time.We stored the results of a given optimization in a "look-up table," such that each entry included the target TNAP and its associated stimulation parameters that resulted in the closest matching ENAP (Fig. 2).

D. Lookup Tables
We repeated the process described in Section C with four sets of constraints, resulting in four stimulation paradigms.First, we allowed the optimization to select PA and PW for each pulse, but required that the same single C-FINE contact be selected for each TNAP in the reachable set.We repeated this for each of the 15 C-FINE contacts, creating 15 singlecontact lookup tables.
Second, we generated an interleaved pattern of stimulation we call the "IMC" paradigm in which only a single contact could provide stimulation at a time, but the active contact could switch on a pulse-by-pulse basis.The single-contact IMC lookup table was created by selecting among the 15 singlecontact lookup tables to find the contact (and stimulation parameters) that led to the lowest error for each TNAP (Fig. 2).
Third, we generated two-contact stimulation paradigms in which two contacts could be active simultaneously.We generated 21 non-IMC two-contact lookup tables using all two-contact pairs of the following C-FINE contacts: M1, M2, M3, M4, M13, M14, and M15.These non-IMC patterns used the same pair of electrodes for all TNAPs in the reachable set.We chose these contacts because they resulted in the lowest error (see Results).
Fourth, we generated a two-contact IMC lookup table using an interleaved pattern in which two contacts could activate simultaneously, but the active contact pair could switch on a pulse-by-pulse basis.The two-contact IMC lookup table was created by selecting amongst the 21 two-contact lookup tables generated as described above.
In summary, these lookup tables were generated (Fig. 2): • 15 single-contact non-IMC lookup tables (one for each of the 15 cathodal C-FINE contacts) • 1 single-contact IMC lookup table that switches between the 15 C-FINE contacts on a pulse-by-pulse basis • 21 two-contact non-IMC lookup tables (one for each pair of M1, M2, M3, M4, M13, M14, and M15) • 1 two-contact IMC lookup table that switches between the selected 21 two-contact pairs on a pulse-by-pulse basis

E. Playlists
To demonstrate the implementation of the IMC approach in conveying touch stimuli via PNS, we created PNS pulse trains, called a playlist, to reproduce a typical touch stimulus encountered in everyday tasks: an indentation stimulus with a depth of 3 mm and a ramp duration of 250 ms (Fig. 1a).Whereas a lookup table holds the C-FINE stimulation parameters for the entire TouchSim reachable set, a playlist holds the pattern of C-FINE stimulation parameters that, when applied to the nerve in sequence, mimic the time course of this indentation stimulus as closely as possible.Thus, each playlist contains an ordered sequence of 1000 sets of stimulation parameters, one for each pulse in a 1000 Hz pulse train to create the 1000 TNAPs in the 1 s touch stimulus.We constructed four versions of this playlist by pulling the stimulation parameters from each of the following lookup tables: 1) single-contact non-IMC lookup table constructed for contact M1 (Fig. 3a), 2) singlecontact IMC lookup table (Fig. 3b), 3) two-contact non-IMC lookup table constructed for contacts M1 and M2 (Fig. 3c), and 4) two-contact IMC lookup table (Fig. 3c-d).Contacts M1 and M2 were chosen for the non-IMC playlists based on their frequency of occurrence in the IMC lookup tables (Fig. 3e) [51].Note that the two-contact IMC paradigm had the option to use single-contact PNS to create an ENAP if it was determined that single-contact PNS would result in a lower reproduction error than two-contact PNS.The optimization algorithm selected single-contact stimulation for 39 % of ENAPs (Fig. 3d).Note that only the two-contact IMC playlist was used for further analyses; the non-IMC and single-contact IMC playlists are only depicted to explain the differences across paradigms.

F. Analyses 1) Comparison of Non-IMC and IMC Methods:
We compared reproduction accuracy between the non-IMC and IMC PNS paradigms to determine any benefits of the interleaved multicontact approach.We hypothesized that using IMC PNS will improve reproduction accuracy compared to non-IMC PNS.
For this analysis, we compared the single-contact IMC lookup table to each of the 15 non-IMC lookup tables.We compared lookup tables rather than playlists so that each TNAP had equal weight in the error calculation since each TNAP is sampled only once in a lookup table, whereas a single TNAP could occur multiple times within a playlist.Thus, the reproduction error across a lookup table more accurately represents the ability of the stimulation paradigm to replicate TNAPs associated with touch in general, whereas a playlist more accurately represents the ability of the paradigm to replicate specific tactile stimuli.
We subtracted the reproduction error of every ENAP in each of the 15 non-IMC lookup tables from the error of the corresponding ENAP in the IMC lookup table to compare magnitude of reproduction error difference between IMC and non-IMC paradigms.Then, we calculated the prevalence of ENAPs that produced a smaller, larger, or equal reproduction error with the IMC lookup table compared to each of the 15 non-IMC lookup tables.We also performed Wilcoxon rank sum tests to see if there was a significant difference (α < 0.05) between the IMC lookup table reproduction errors and each of the 15 non-IMC lookup table reproduction errors.

2) Comparison of Single-Contact and Two-Contact IMC PNS:
We then compared the reproduction accuracy between the single-and two-contact IMC PNS paradigms to determine any benefits of using two-contact stimulation.We hypothesized that using two-contact IMC PNS will improve reproduction accuracy compared to single-contact IMC PNS.
For this analysis, we compared the single-contact IMC lookup table and the two-contact IMC lookup table.We subtracted the reproduction error of every ENAP in the single-contact IMC lookup table from the error of the corresponding ENAP in the two-contact IMC lookup table to compare magnitude of reproduction error difference between single-contact and two-contact IMC paradigms.Then, we calculated the prevalence of ENAPs that produced a smaller, larger, or equal reproduction error with the two-contact IMC lookup table compared to the single-contact IMC lookup table.We also used the Wilcoxon rank sum test to test for a significant difference (α < 0.05) between the reproduction errors of the single-contact and two-contact IMC lookup tables.
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

3) Comparison of IMC Paradigm and Existing PNS
Paradigms: Finally, we sought to compare our novel IMC paradigm to existing PNS paradigms reported in the literature.We compared each existing PNS paradigm to the two-contact IMC paradigm, as this IMC paradigm performed the best in our previous analyses (see Results).For this analysis, we compared playlists rather than lookup tables, as existing PNS paradigms in literature are presented in terms of their ability to replicate touch stimuli rather than their ability to replicate TNAPs in general.We hypothesized that the IMC paradigm would more closely mimic the neural code of tactile sensation, resulting in lower reproduction error, than existing PNS paradigms.
For the non-biomimetic PNS paradigms, we selected stimulation parameters from Table I.The fixed-pulse paradigm [4], [6], [8] had fixed stimulation parameters across the entire stimulus duration (Fig. 5b).We fixed PW at PW 25% = 0.05 ms and PA at PA TNAP = 0.091 mA.The force-based paradigm [5], [7] varied PA based on the indentation depth throughout the stimulus (Fig. 5c).Thus, to reproduce our target stimulus, the PA ramped up linearly with the force exerted on the finger, then held constant, then ramped down.We fixed PW at PW 25% = 0.05 ms.We set the upper limit of PA to PA TNAP = 0.091 mA and the lower limit to PA th = 0.07 mA.The sinusoidal paradigm [6], [7], [8] varies PW over time based on the amplitude of a 1-Hz sine wave (Fig. 5d).We fixed PA at PA 25% = 0.125 mA.We set the upper limit of PW to PW TNAP = 0.035 ms and the lower limit to PW th = 0.025 ms.
The two biomimetic paradigms [11], [18] attempted to match the neural population size of the TNAP with electrical stimulation, but in different ways.Biomim 1 chose stimulation parameters by directly matching the overall neural population size for each stimulation pulse, while Biomim 2 linearly scaled the PA of the stimulation based on a percentage of the maximum neural population size.For Biomim 1 [18], we fixed PA at PA 25% = 0.125 mA.We then used the electrical activation model to find a PW that activated an ENAP with the closest match in neural population size to that of the TNAP at each millisecond of the neural response (Fig. 5e).
For Biomim 2 [11], we specifically mimicked the neural population size component of HNM2.We first calculated the neural population size of the TNAP at each millisecond and found the maximum neural population size (NPS max ) of any TNAP across the stimulus, which was 58 neurons for the chosen indentation stimulus.We then used Equation ( 2) to calculate PA for PNS at each millisecond of the stimulus.We fixed PW to PW 25% = 0.05 ms, and we set PA max in (2) to PA TNAP = 0.091 mA.The PA calculated by (2) was rounded to three decimal places to match our 0.001 mA step size used in the other paradigm simulations (Fig. 5f).

P A =
N P S N P S max P A max (2) As described above, the IMC paradigm uses optimization algorithms to select optimal PWs and PAs for each TNAP in a stimulus over the range of 0 ms to PW max = 0.2 ms and 0 mA to PA max = 0.5 mA, respectively (Fig. 5a).Note that the IMC paradigm is the only paradigm that varied PW, PA, and contact on a pulse-by-pulse basis.The other paradigms only vary at most one stimulation parameter at each millisecond and use the same contact for the entire stimulus (Fig. 5b-f).
To calculate reproduction error, we first used the electrical activation model to generate neural responses for each existing PNS paradigm (Fig. 5).We then calculated the reproduction Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I STIMULATION PARAMETERS USED FOR PNS PARADIGMS
error of each existing paradigm at each millisecond of the stimulus, using the error metric as in Eqn 1.We also calculated the reproduction error in number of neurons broken down by afferent type (SA, RA, PC) for each paradigm.
Then, we subtracted the reproduction error of every ENAP in each of existing PNS paradigm playlists from the error of the corresponding ENAP in the IMC playlist to compare magnitude of reproduction error difference between the IMC paradigm and existing paradigms.We also determined the prevalence of ENAPs that had a smaller, equivalent, or larger reproduction error when using the IMC paradigm compared to each of the existing paradigms.Finally, we used the Wilcoxon rank sum test to test for a significant difference (α < 0.05) between the IMC paradigm and each existing paradigm.

A. The IMC Method More Accurately Reproduced Neural Activity Than the Non-IMC Method
First, we compared the reproduction error of each of the 15 single-contact non-IMC lookup tables to the single-contact IMC lookup table.While all paradigms in this analysis used our modeling and optimization approach to select stimulation parameters, only the IMC paradigm involved interleaving stimulation contacts on a pulse-by-pulse basis.The single-contact IMC lookup table had significantly lower reproduction errors than all 15 single-contact non-IMC lookup tables (Wilcoxon rank sum test, p<0.001) (Fig. 4a).The single-contact IMC reproduction error was less than or equal to the 15 singlecontact non-IMC reproduction errors for all ENAPs (Fig. 4c).Thus, we conclude that IMC PNS outperformed non-IMC PNS and used IMC PNS for all remaining analyses.
While no single-contact non-IMC ENAPs had lower error than the IMC ENAPs, non-IMC paradigms with contacts M1-M3 and M13-M15 had a higher proportion of ENAPs with error equivalent to that of the IMC paradigm (16-57 %), compared to contacts M4-M12 (1 %) (Fig. 4c).That is, contacts M1-M3 and M13-M15 generally reproduced TNAPs more accurately than contacts M4-M12.This is why these contacts were selected to generate the two-contact IMC approach.

B. Two-Contact IMC PNS More Accurately Reproduced Neural Activity Than Single-Contact IMC PNS
After determining that IMC PNS outperformed non-IMC PNS, we then compared the single-contact and two-contact IMC lookup table reproduction errors to determine the influence of field shaping on ENAP reproduction accuracy.The two-contact IMC lookup table had significantly lower reproduction error than the single-contact IMC lookup table, (Wilcoxon rank sum test, p<0.001) though the numerical difference in error was very small (Fig. 4b).The two-contact IMC paradigm had smaller, equivalent, or larger error than the one-contact IMC paradigm for 22%, 77%, and <1% of ENAPs, respectively (Fig. 4d).Thus, we conclude that two-contact IMC PNS more closely mimicked touch neural responses than single-contact IMC PNS, and we used twocontact IMC PNS for all remaining analyses.

C. The IMC Paradigm More Accurately Reproduced the Neural Code of Tactile Sensation Than Existing PNS Paradigms
We compared the two-contact IMC paradigm playlist to three non-biomimetic paradigm playlists (fixed-pulse, forcebased, and sinusoidal) and two biomimetic paradigm playlists (Biomim 1 [18] and Biomim 2 [11]) for a representative indentation stimulus delivered to the index finger (Fig. 1a).The stimulation parameters for each paradigm and the neural activation pattern resulting from each paradigm are shown in Fig. 5a-f.The natural neural activation for this touch stimulus as predicted by TouchSim is shown in Fig. 1b for comparison.
We calculated the reproduction error for each paradigm and found that across all paradigms, the reproduction error was highest for RA neurons, intermediate for SA neurons, and lowest for PC neurons (Fig. 5, right column).As expected, this demonstrates that the relative contribution of each afferent type to the overall reproduction error correlates with its proportion in the neural population [1].We then compared the reproduction errors of each of the existing non-biomimetic and biomimetic PNS paradigms to the IMC paradigm.The twocontact IMC paradigm had significantly smaller reproduction errors than all five existing paradigms (Wilcoxon rank sum test, p<0.001) (Fig. 6a).While the two-contact IMC paradigm had lower error than any existing paradigm in the majority of all ENAPs, the existing biomimetic paradigms performed better than the existing non-biomimetic paradigms (Fig. 6b).The two-contact IMC paradigm had less reproduction error in a higher proportion of ENAPs when compared to the non-biomimetic existing paradigms (98-100 %) than when compared to the biomimetic existing paradigms (41-59 %).In addition, the two-contact IMC paradigm was more likely to have an equivalent reproduction error to the existing biomimetic paradigms (41-59 % of ENAPs) compared to the existing non-biomimetic paradigms (0-3 % of ENAPs) (Fig. 6b).Still, the two-contact IMC ENAPs rarely had reproduction errors larger than the existing paradigm ENAPs (<1 %).In fact, this only occurred for one ENAP in one existing paradigm (Biomim 2).Thus, we conclude that the two-contact IMC paradigm performed better than the existing PNS paradigms we simulated, including the previously reported biomimetic paradigms.

D. Discussion 1) Interleaved Multi-Contact Stimulation Most Closely
Approximates Natural Neural Activity: We showed that IMC PNS more accurately reproduced the neural response to touch stimuli than non-IMC PNS, demonstrating that interleaved multi-contact biomimetic stimulation outperforms non-interleaved single-contact stimulation.Importantly, the non-IMC PNS paradigms were also biomimetic, since they were generated using the same optimization approach to approximate natural neural activity as the IMC PNS, yet they still resulted in higher error than the interleaved approach.We believe that the critical factor was that the IMC paradigm was able to recruit numerous different populations of neurons throughout the stimulus, whereas traditional biomimetic (i.e.non-interleaved) approaches can only scale the size of a single recruited population.Prior interleaved stimulation approaches to reduce fatigue in motor neuroprostheses were similarly founded on the concept of activating smaller sub-groups of neurons at different times, although these prior approaches did so by cyclically repeating a set of contacts rather than attempting biomimicry [36], [37], [38], [39].While a prior study reported an interleaved surface stimulation approach for sensory feedback [41], this approach was based on the idea of modulating the activation threshold of a single neural population by applying stimulation through two contact pairs with sub-millisecond delays, rather than activating multiple sub-groups of afferents by successive pulses.
Our results also show that the two-contact IMC paradigm resulted in significantly smaller reproduction error than the single-contact IMC paradigm.This is likely because the twocontact approach produced field shaping to enable selective activation of even more specific groups of neurons than can be achieved with single-contact stimulation [25], [26], [27], [28].While no prior biomimetic PNS paradigm applied stimulation through multiple contacts simultaneously, these approaches are beginning to be evaluated for biomimetic brain stimulation using groups of up to four microelectrode contacts simultaneously [35].Assuming the initial clinical results from intracortical microstimulation will extend to PNS, multi-contact biomimetic PNS may be able to improve force discriminability and dynamic range of evoked touch percepts.However, one drawback to multi-contact stimulation is that it is likely to require higher power consumption, since more electrical charge could be delivered throughout the stimulus.
Given that our paradigm benefited the ability to apply interleaved stimulation via contacts distributed around the nerve, as well as the ability to stimulate with multiple contacts at a time, further studies should examine the impact of contact number, spacing, and size on the implementation and outcomes of IMC biomimetic stimulation.These findings could influence future multi-contact electrode designs.
2) The IMC Paradigm Outperforms Prior PNS Paradigms: We compared the performance of our novel IMC paradigm to several existing PNS paradigms previously reported in the literature.We found that the IMC paradigm outperformed every simulated existing paradigm, including both biomimetic and non-biomimetic approaches, by producing lower reproduction error.
In natural touch, tactile sensation produces a dynamic neural response such that the size and composition of the recruited neural population varies throughout the stimulus [1].We found that prior biomimetic paradigms [11], [18] resulted in smaller reproduction error than the non-biomimetic fixed-pulse, force-based, and sinusoidal paradigms.The prior biomimetic paradigms have better performance than the nonbiomimetic paradigms because they are attempting to replicate some aspect of the neural code, such as neural population size and firing rate.However, they do not consider neuron type or neuron location, which likely contribute to their lower performance as compared to the IMC paradigm, which did Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
account for these neural codes in addition to neural population size and firing frequency.In addition, the IMC paradigm allows pulse-by-pulse switching of the stimulating contact(s), as well as the use of multi-contact stimulation for a given pulse, to enable selective recruitment of specific neural populations throughout the stimulus, while the prior paradigms only involved a single electrode contact for a given stimulus.
However, the importance of lower reproduction error in creating more natural or more informative sensation has not yet been confirmed in clinical tests.Further, it is unknown how low reproduction error needs to be for the evoked neural activity to be "close enough" to natural activity to yield a natural feeling sensation.Future clinical testing of the IMC paradigm and direct comparisons to other paradigms would allow a better understanding of the relationship between reproduction error and sensation experience.Additionally, our results are based on neural activation patterns predicted from a biophysical model, which has not yet been validated.Neurophysiological studies utilizing neuron-specific recording techniques, such as microneurography, are needed to confirm that each stimulation paradigm produces the predicted neural firing patterns in vivo.
3) Future Improvements to the Optimization Approach: In the single contact analyses, electrode contacts M1-M3 and M13-M15 had smaller errors than contacts M4-M12 and were preferentially selected for the IMC paradigm.This might be explained by the proximity of these contacts to the fascicles of interest.Interestingly, the optimization algorithm selected near-contacts rather than far-contacts even though our error metric did not penalize the recruitment of off-target neurons (i.e.neurons in the median nerve that innervate other hand regions other than the index finger).Future research should examine the effect of the inclusion of a penalty for off-target neural recruitment in the optimization function [51].
Another factor impacting contact selection was the bias of the optimization algorithm toward lower numbered contacts.In times where reproduction error was equivalent between one or more contacts, the optimization algorithm would select the lower numbered contact.For example, if an ENAP had the same error when using contacts M1 and M15, contact M1 would be selected by the algorithm.In scenarios where multiple contacts output the same reproduction error for an ENAP, the contact selection process could be improved by adding a secondary optimization term after reproduction error is minimized, such as selecting the contact that produces the minimum charge rather than the contact with the smallest numerical label.
Regarding the use of field shaping in the IMC approach, we hypothesize that enabling the optimization algorithm to select more than two contacts at a time (i.e. three or more simultaneous multi-contact PNS) has the potential to further reduce reproduction error.Our optimization algorithm can easily be expanded to account for additional stimulating contacts, but with a tradeoff in computation time.Another benefit of the field shaping approach is that the algorithm can still choose to use fewer contacts for a given pulse, and thus may select a single contact if it would result in a lower reproduction error.
4) Adapting the IMC Paradigm for Real-World Use: Ultimately, our goal is to implement the IMC paradigm into real-time sensory encoders to be used to deliver sensory stimulation to augment functional task performance for people with neurological disorders or disabilities.For example, an upper limb prosthesis user could receive IMC stimulation to convey naturalistic, informative haptic feedback from pressure sensors on the prosthetic fingertips.
First, we will need to test the IMC paradigm in human participants to characterize its perception.To implement the IMC paradigm in clinical studies, we should develop patientspecific electrical activation models based on the participant's fascicular structure, neuroanatomy, and electrode configuration.We would then couple neural activation predictions from the novel IMC paradigm with clinical assessments of the perceived location of sensations produced by single-contact stimulation to identify target fascicles in the nerve corresponding to specific hand locations.We would also need to scale the stimulation parameters in the model to real stimulation parameters based on clinical assessments of sensory detection threshold and maximum limits.To use the IMC paradigm during functional task performance, we would need to reconfigure the optimization approach to select stimulation parameters on a pulse-by-pulse basis in real-time based on an incoming sensor signal.In contrast, in the approach reported here, optimization was used offline to select stimulation parameters for an idealized indentation stimulus.The optimization algorithm would need to be sped up and must be able to account for sensor noise or artifact.Moreover, it would be ideal if the optimization algorithm could select stimulation parameters and electrode contact in a single step, rather than sequentially as reported here (see Fig. 2).This would increase the difficulty of solving the optimization problem, increasing computational complexity and time.
In addition, we must map TNAPs to parameters of natural stimuli that can be measured with real-time sensors.To approach this problem, we could use our existing dataset to correlate TNAPs with specific depths in the ramp-up, hold, and ramp-down phases of the indentation stimuli.This would be challenging, however, as multiple TNAPs produce the same depth in each phase of our TouchSim-generated neural responses.The most important features, or in this case neurons, that contribute to a specific sensation must be selected and reproduced from all TNAPs that produce the same depth at a stimulus.The TNAPs that came beforehand must also be considered, as current TNAPs depend on previous TNAPs rather than being independent at each millisecond.

IV. CONCLUSION
In conclusion, we present a new approach to biomimetic stimulation.The IMC paradigm is based on biophysical models of neural activation and can approximate the neural response to tactile stimuli with high accuracy.The IMC paradigm's key innovation is that it includes pulse-by-pulse updating of electrode contact and stimulation pulse parameters to selectively recruit different small sub-fascicular afferent populations throughout the stimulus time course.The IMC paradigm can also include multi-contact field shaping approaches, which leads to better performance than single contact stimulation.Overall, the IMC mimics critical neural coding properties better than existing PNS paradigms, including previously reported biomimetic PNS paradigms.Future studies will test the IMC paradigm in human subjects to determine whether it leads to improvements in the perceptual experience and information content of generated sensations.If human subject testing yields improved perceptual performance, the IMC paradigm will then be adapted for real-time use in sensor-based neuroprosthetic devices.We envision that the IMC paradigm will provide natural sensory feedback to those who have lost tactile sensory capabilities, greatly improving their functionality and overall quality of life.

Fig. 1 .
Fig. 1.(a) Example indentation stimulus with a 3 mm depth and a 250 ms ramp duration.The star ( * ) denotes the timing of the example TouchSim neural activation profile (TNAP) displayed in panel d.(b) Raster plot of the TouchSim-generated neural response to the touch stimulus in panel a.Color indicates afferent type.(c) 16-contact Composite Flat Interface Nerve Electrode (C-FINE) placed around a nerve.(d) Layout of all low-threshold mechanoreceptive afferents in the index finger (grey) showing neurons active in the starred TNAP (purple).(e) Cross section of a portion of the nerve with surrounding C-FINE showing a subset of fascicles.Neurons active in the example electricallyactivated NAP (ENAP) are highlighted in purple (grey neurons are inactive).

Fig. 3 .
Fig. 3. Stimulation parameters selected by the optimization algorithm for example non-IMC and IMC paradigms.(a) PA and PW for the first 50 ms of a single-contact (1C) non-IMC playlist for the example touch stimulus in Fig. 1a.Note that for non-IMC PNS, the same contact is used for every pulse.(b) PA and PW for the 1C IMC playlist for the same stimulus.Color denotes the stimulating contact for the pulse (legend in panel e).Note that for IMC PNS, the stimulating contact can change on a pulse-by-pulse basis throughout the stimulus.(c) Contact usage for an example two-contact (2C) non-IMC playlist and the 2C IMC playlist.With 2C non-IMC PNS, the same two contacts are used for each pulse.With 2C IMC PNS, two contacts can be active during each pulse, and the two-contact pair can vary across pulses.(d) Selection frequency of 1C (39 %) and 2C (61 %) across pulses in the 2C IMC playlist for the touch stimulus.(e) Breakdown of contact usage in the 1C non-IMC, 2C non-IMC, 1C IMC, and 2C IMC playlists shown in panels a-c.

Fig. 4 .
Fig. 4. Reproduction error comparisons for lookup table ENAPs across optimization-based biomimetic paradigms.Significant differences (Wilcoxon rank sum test, p<0.001) are denoted with a star ( * ).Note that outliers were removed from boxplots.(a) Comparison of the magnitude of reproduction error between each of the 15 single-contact (1C) non-IMC paradigms and the 1C IMC paradigm.Positive error indicates that IMC error was greater than non-IMC error, while negative error indicates that IMC error was less than non-IMC error.(b) Comparison of the magnitude of reproduction error between the 1C IMC and two-contact (2C) IMC paradigms.Positive error indicates that 2C error was greater than 1C error, while negative error indicates that 2C error was less than 1C error.(c) Percentage of ENAPs in the 1C non-IMC look-up tables that had error greater than, equal to, or less than the corresponding IMC ENAPs.Note that the IMC paradigm never had a larger reproduction error than any of the non-IMC paradigms for any ENAP.(d) Percentage of ENAPs in the 1C IMC look-up table that had error greater than, equal to, or less than the corresponding 2C IMC ENAPs.

Fig. 5 .
Fig. 5. Stimulation playlists reproducing the touch stimulus shown in Fig. 1a.Left: Stimulation parameters selected for each PNS paradigm.PA is shown in grey on the left axis and PW is shown in purple on the right axis.Middle: Raster plots of neural responses to each stimulation pulse train.Color denotes neuron type (SA, RA, PC).Right: Reproduction error magnitude for each paradigm broken down by neuron type (SA, RA, PC).Each row is a different stimulation paradigm: (a) Two-contact IMC PNS.(b) Fixed-pulse PNS.(c) Force-based PNS.(d) Sinusoidal PNS.(e) Biomimetic PNS Approach 1 (Biomim 1).(f) Biomimetic PNS Approach 2 (Biomim 2).

Fig. 6 .
Fig.6.Comparison of reproduction error between the interleaved multi-contact (IMC) approach and existing (EP) of peripheral nerve stimulation.(a) Reproduction error magnitude between the two-contact IMC playlist and each EP playlist.Positive error indicates that IMC error was greater than EP error, while negative error indicates that IMC error was less than EP error.(b) Percentage of pulses in the playlist in which the IMC error was equal to, less than, or greater than the EP error.Significant differences (Wilcoxon rank sum test, p<0.001) are denoted with a star ( * ).Note that outliers were removed from boxplots.EP key: B: Fixed-pulse PNS.C: Force-based PNS.D: Sinusoidal PNS.E: Biomimetic PNS approach 1 (Biomim 1).F: Biomimetic PNS approach 2 (Biomim 2).