Objective: Peripheral nerve stimulation (PNS) is most effective when specific nerve fiber subpopulations are activated, while minimizing off-target activation, which may cause undesirable side effects. This selectivity depends primarily on electrode design and charge delivery, but it remains unclear how selectivity depends on electrode location within a fascicle and whether intrafascicular field steering can shape recruitment of different nerve fibers. Our objective was to test the hypothesis that selective PNS could be achieved through electrode placement and intrafascicular electric field steering using Longitudinal Intrafascicular Electrodes (LIFEs). Approach: LIFEs were implanted into the tibial fascicle of the sciatic nerve of 17 anesthetized adult rats. We tested whether electrodes positioned at different cross-sectional and longitudinal locations within the same fascicle, together with different electric field-steering approaches, produced distinct activation patterns in the gastrocnemius lateralis muscle. Muscle responses were measured using high-density epimysial electromyography (HD-eEMG). Main Results: Electrodes placed at different locations within the same fascicle activated distinct muscle regions, demonstrating intrafascicular selectivity. Bipolar stimulation recruited nerve fibers differently than monopolar stimulation, showing that electric field steering can further shape the selective recruitment. In both configurations, increasing the stimulation amplitude produced a graded increase in muscle activation. Furthermore, our findings demonstrated that HD-eEMG is an effective tool for evaluating intrafascicular selectivity. Significance: These findings suggest that improving on-target selectivity may support next-generation bioelectronic therapies with better outcomes and fewer side effects, potentially enabling more precise, organ-specific neuromodulation. Using multiple intrafascicular electrodes may provide two complementary strategies for enhancing selectivity: strategic intrafascicular placement to access different fiber subpopulations and bipolar configurations to steer recruitment beyond what a single electrode can achieve.

Data Processing

Monopolar EMG recordings were spatially filtered to decrease the detection volume and isolate the electrical activity of the underlying muscle, as described by Gallina et al. (29) (Figure 1c). A single-differential spatial filter was applied along the direction of the muscle fibers (i.e., the columns) by subtracting each channel from its immediately adjacent channel in the column direction (i.e., Channel _(r,c) = Channel _(r,c) – Channel_(r+1,c), where r denotes the row and c denotes the column). The resulting 26 channels of spatially filtered data were then filtered in the frequency domain using a fourth-order Butterworth bandpass filter with cut-off frequencies of 10 and 1000 Hz. Data from two HD-eEMG channels in one rat were excluded from the analysis because their signals were dominated by recording noise and did not exhibit physiologically meaningful variation.

For each stimulation level, a stimulus-triggered average (STA) of eEMG was calculated for each HD-eEMG channel to construct M-wave responses. To characterize these responses, the root mean square (RMS) amplitude of each STA M-wave was computed across a 16-ms window, beginning with the onset of the evoked response.

To evaluate the spatial distribution of eEMG amplitude (rather than the magnitude) across the electrode grid, we normalized STA RMS values for all eEMG channels for a given electrode configuration and stimulation level to the highest of those RMS values. Spatial heatmaps of stimulation-evoked M-wave amplitude were constructed from the normalized RMS values for data visualization (Figure 1c). Next, eEMG channels were sorted in descending order of their normalized RMS values, which were then plotted against the sorted channel numbers. This created an amplitude rank profile for each stimulation level and electrode configuration for each rat that illustrated the distribution spread of RMS values among the channels. The area under the curve (AUC) for each amplitude profile was calculated using trapezoidal numerical integration (MATLAB, trapz).

The AUC metric was chosen over alternative metrics such as spatial entropy or localization indices as the primary measure of spatial selectivity. It provides a single, continuous index of how concentrated or distributed muscle activation is across the array. Lower AUC values indicate that eEMG amplitude was concentrated among few channels, reflecting greater selectivity, compared with higher AUC values. Higher AUC values indicate that eEMG amplitude was more evenly distributed across channels, reflecting less selectivity. Because AUC quantifies the degree of spatial concentration but not the anatomical location of activation, we used two complementary analyses to capture additional information about the spatial activation, identifying the location of peak activation, which was presented in histograms, and calculating pairwise spatial correlations (Pearson r) between electrode pairs implanted within the same fascicle to assess overlap in muscle activation patterns across LIFEs.

Statistical Analyses

All statistical analyses were conducted using R (version 4.5.0, R Core Team, 2024), with statistical significance set at ⍺ ≤ 0.05. Linear mixed-effects regression models with fixed and random effects detailed below were employed to address our research questions. These models were implemented using the lme4 package utilizing Restricted Maximum Likelihood Estimation. Where applicable, post hoc comparisons were performed using estimated marginal means calculated with the emmeans package, with adjustments for multiple comparisons (Tukey or Dunnett). Model assumptions, including linearity and normality of residuals, were evaluated to ensure the validity of model fitting. The model equations are presented in R syntax rather than in standard linear mixed model notation. Prior to fitting the statistical model, stimulation level values were mean-centered to facilitate interpretation of the regression coefficients. By including mean-centered stimulation level, the model intercept represented the eEMG amplitude response at the mean stimulation level, rather than at a stimulation level of zero.

Effect of stimulation level on eEMG amplitude

To determine the effect of monopolar stimulation amplitude level on eEMG RMS amplitude across channels, we constructed a linear mixed-effects model with fixed effects of Stimulation Amplitude_c (mean-centered), EMG Channel, and the Stimulation Amplitude_c-by-EMG Channel interaction. To fully account for statistical dependencies in the data due to its hierarchical structure, we included random intercepts for Rat, Channel nested within Rat, and LIFE nested within Rat, with a random slope for Stimulation Amplitude_c by LIFE nested within Rat.

(Model 1): RMS ~ 1 + StimAmplitude_c*Channel + (1│Rat)+ (1│Rat:Channel)+ (1+StimAmplitude_c |Rat:LIFE)

Effect of stimulation from different LIFEs on eEMG amplitude

To evaluate the effect of monopolar stimulation delivered by different LIFEs at different stimulation levels, the model from the previous section (the effect of Stimulation Amplitude_c-by-EMG Channel on eEMG RMS amplitude) was used again. Our experimental design precluded inclusion of LIFE as a fixed effect, given that LIFE placement was not anatomically standardized across animals (i.e., LIFE #1 in one rat was not implanted at a location equivalent to LIFE #1 in another rat). Therefore, we examined the influence of different LIFEs by comparing the model variances due to LIFE, Stimulation Amplitude_c, Rat, and Channel in relation to the global average RMS amplitude. To do so, we extracted the corresponding variances of the random intercepts and slope. This allowed us to evaluate which factor caused the most variation in eEMG amplitude (i.e., using different LIFEs, using different stimulation amplitudes, recording from different channels, or recording from different rats).

Effect of stimulation  level on the spatial distribution of eEMG amplitude

To determine the effect of stimulation level on the spatial distribution of eEMG (quantified as AUC of the eEMG amplitude profiles, we constructed a linear mixed effects model with fixed effects of Stimulation Amplitude_c, random intercepts of Rat and LIFE nested within Rat, and** a** random slope of Stimulation Amplitude_c by LIFE nested within Rat. After visualizing the AUC data vs. stimulation level, we also included polynomial terms in the model, with fixed effects of Stimulation Amplitude_c² and Stimulation Amplitude_c³ and random slopes of the polynomial terms by LIFE nested within Rat.

(Model 2): AUC ~ 1 + StimAmplitude_{c + StimAmplitudec^2 ^}+ StimAmplitude_c^3 ^+ (1│Rat) +(1+StimAmplitude_c+StimAmplitude_c²+StimAmplitude_c³│Rat:LIFE)).

We also examined whether the shape and variability of the ranked amplitude profiles changed with stimulation level. At each stimulation level, we ensemble averaged the ranked amplitude profiles from all LIFEs and rats. We characterized the shapes of the resulting average profiles using polynomial regression and examining the coefficients across stimulation levels. The average amplitude profile for each stimulation level was fitted to a cubic polynomial using least squares in the following form:

RMS_Norm (x)=β_₀+β₁x+β₂x⁰+β₃x³

where β₀ is the intercept term, β₁-β₃ are the coefficients of the polynomial curve, and is the rank of the channel (e.g. 1^st, 2^nd… rather than the channel number).

To examine the inter-animal variability of the ranked amplitude profiles at each stimulation level, we ensemble averaged the amplitude profiles from all LIFEs for each rat and stimulation level. We then calculated the standard deviation of the mean profile across all sorted channels. This resulted in a standard deviation profile by sorted channel, the AUC of the standard deviation profile (trapz, MATLAB), and a calculated mean standard deviation among channels.

Effect of stimulation type on eEMG amplitude

To evaluate the effect of stimulation type (monopolar or bipolar) on eEMG RMS amplitude across channels and at different stimulation levels, we constructed a linear mixed effects model with fixed effects of Stimulation Type, Stimulation Amplitude_c, EMG Channel, and their three-way interaction.

To fully account for statistical dependencies within the data, we included random intercepts for Stimulation Type nested within Rat, Channel nested within Rat, and a random slope for Stimulation Amplitude_c by Stimulation Type nested within Rat. A random intercept of Rat was included initially but then removed because it accounted for zero variance and created a model singularity.

(Model 3):RMS ~ 1+ StimTypeChannelStimulation Amplitude_c+(1│Rat:Channel)+(1+Stimulation Amplitude_c│Rat:StimType)

Effect of stimulation type on eEMG spatial distribution

We modeled stimulation type as having three levels: bipolar stimulation with two LIFEs, monopolar stimulation with one of the two LIFEs (labeled MS1), and monopolar stimulation with the other of the two LIFEs (labeled MS2). To determine whether spatial distribution of the eEMG (quantified as the AUC of the eEMG amplitude profiles), differed among these three stimulation configurations and stimulations levels, We constructed a linear mixed effects model with fixed main effects of Stimulation Type, Stimulation Amplitude_c, and their** two-way** interaction. We included a random intercept for Rat and a random slope for Stimulation Amplitude_c by Stimulation Type nested within Rat. To be consistent with Model 2 above, we also included polynomial terms in the model, with fixed effects of Stimulation Amplitude_c² and Stimulation Amplitude_c³ and random slopes of the polynomial terms by Stimulation Type nested within Rat.

(Model 4): AUC ~ 1 + StimTypeStimulation Amplitude_c+ StimTypeStimulation Amplitude_c²+Stim TypeStimulation Amplitude_c³+(1│Rat)+(1+Stimulation Amplitude_c+Stimulation Amplitude_c²+Stimulation Amplitude_c³│Rat:StimType).*

Effect of stimulation level on the HD-eEMG channel with the highest eEMG amplitude

For each stimulation level, the channel with the highest RMS value for each LIFE was identified, and the number of occurrences per channel was counted to build histograms and perform the analyses. Analyses were performed separately for the two types of stimulation (monopolar and bipolar) at each stimulation level.

Effect of electrode placement and stimulation level on spatial correlation of eEMG activity elicited by LIFEs.

To quantify overlap between muscle activation patterns produced by different LIFEs, pairwise spatial correlation (Pearson r) was calculated between electrode pairs within the same fascicle at the same stimulation level. Because stimulation levels did not align exactly across LIFEs, normalized RMS values for each channel were interpolated onto a common set of stimulation levels (xTH), restricted to each LIFE’s tested range without extrapolation.

The normalized RMS response at each of the 26 HD-eEMG channels elicited by each electrode in the pair were compared on a channel-by-channel basis (i.e., for Rat “5”, Channel 1 of LIFE 1 with Channel 1 of LIFE 2 at 1.2xTH) and a single Pearson correlation coefficient was computed across all 26 paired values, resulting in one correlation coefficient per electrode pair per stimulation level. This was repeated across all stimulation levels to assess whether spatial activation patterns between electrode pairs varied with amplitude. Comparisons at a given stimulation level were excluded if, for both LIFEs, every EMG Channel’s RMS value fell below 5%of that LIFE’s own maximum RMS, to avoid computing correlations on signals that were predominantly noise.