Introduction

Muscle force generation and movement originate from commands issued by the central nervous system. These commands are transmitted through descending cortical, subcortical, and spinal pathways that con-verge on α-motor neurons in the anterior horn of the spinal cord, which transmit the neural signals to their respective innervated muscle fibers, forming individual motor units (MUs). The recruitment and derecruit-ment of MUs are key neural mechanisms of force modulation and are traditionally described by Hennemańs size principle, which states that MUs are recruited in an orderly fashion from the smallest to the largest (Hen-neman (1957); Henneman et al. (1965, 1974)), primarily determined by the intrinsic input resistance to the motoneurons. The firing activity of recruited MUs is further described as an “onion-skin” scheme, where earlier-recruited MUs maintain systematically higher firing rates than later-recruited MUs (De Luca et al. (1982); De Luca and Hostage (2010); Masakado et al. (1995); Hu et al. (2014)). These principles provide an efficient mechanism for graded force production and precise control particularly for small force levels and are thought to minimize fatigue (Bawa et al. (2014); Conwit et al. (1999); Stein et al. (2005)). They have been supported in a variety of muscles and tasks, including slow isometric contractions (Milner-Brown et al. (1973); Thomas et al. (1986, 1987); Riek and Bawa (1992); Feiereisen et al. (1997); Scutter and Türker (1998); Conwit et al. (1999)), different contraction speeds (Desmedt and Godaux (1977); Calancie and Bawa (1985)), and changes in muscle length (Garland et al. (1996); Stotz and Bawa (2001)). However, studies have reported some variability in recruitment when considering different functional tasks in the form of task-related selective recruitment (Butler et al. (2005); Saboisky et al. (2006); Oßwald et al. (2025)) and changes in RO within groups of MUs in the non-human (Marshall et al. (2022)) and human biceps brachii (Ter Haar Romeny et al. (1982); Tax et al. (1989)) and FDI muscle (Desmedt and Godaux (1981); Thomas et al. (1986)). These studies suggest that the central nervous system can optimize MU recruitment thresh-olds depending on task demands. Limited by the concurrent sampling of only small numbers of MUs with single-channel intramuscular EMG electrodes as used in these studies, the neural determinants of such RO flexibility currently remain unclear. Desmedt and Godaux (1981) speculated that the altered RO they ob-served in FDI during index finger abduction versus flexion may arise from differences in the distribution of descending cortical drive. Notably, the FDI serves distinct mechanical and functional roles in these two movements. During abduction, it is the sole muscle with a moment arm in that direction and therefore acts as the prime mover. However, the FDI also has a flexion moment arm and functions as a synergist to the long extrinsic flexors during index finger flexion, contributing to metacarpophalangeal joint stabilization together with the first lumbrical muscle (Masquelet et al. (1986); Liss (2012); Ranney and Wells (1988)). In addition to further insight into the mechanisms of MU control, understanding task-related changes in recruitment and rate coding of MUs, particularly for hand and digit actions, has practical implications. Real-time MU identification has gained traction in myocontrol and human–machine interfaces, where reliability and flex-ibility of recruitment patterns are critical (Campbell et al. (2025); Farina et al. (2017); Chen et al. (2021); Capsi-Morales et al. (2025); Grison et al. (2024, 2025)).

In this study, we investigate the variability of MU recruitment, discharge rate and common input in the FDI muscle by comparing within and across isometric contractions of the index finger in abduction and flexion direction, which corresponds to two different functional roles of the FDI muscle as described above. Using state-of-the-art HD-iEMG electrode arrays, we are able to identify large populations of concurrently active MUs and track their activity across contraction directions. By matching the global level of FDI muscle activity across contractions, we can observe the recruitment activity of MUs on a larger scale and in more controlled conditions than previously possible. The increase in bandwidth of MU activity correlation measures in larger populations of identified MUs (Negro and Farina (2012)) allows us to apply intramuscular coherence analysis to explore potential neural origins of the neural drive received by FDI MUs that might underlie the observed recruitment flexibility.

Materials and Methods

Participants

Six healthy volunteers (mean age 28 ± 3.8 years; 2 male, 4 female), all right-handed, participated in the study. None had a history of neurological or neuromuscular disorders and fulfilled standard inclusion and exclusion criteria in neuromuscular studies in healthy participants. The study protocol was approved by the ethics committee of the Friedrich–Alexander University Erlangen–Nürnberg (approval number 23-271-B) and was conducted in accordance with the Declaration of Helsinki. All participants provided their informed written consent prior to participation.

Experimental Setup

2D Index Finger Dynamometer

A custom-built dynamometer was developed to measure isometric forces of the index finger in all possi-ble contraction directions. Detailed images of the dynamometer and the 2D force feedback are shown in Supplementary Figure S1A. The setup is adjustable to different arm lengths, hand sizes, and handedness of participants. During the experimental protocol, the arm is positioned in the dynamometer with an elbow flexion angle of approximately 100^◦ and the arm pointing forward. The hand and digits are arranged in a “pointing” posture, with the palm facing medially and the index finger pointing forward, parallel to the ground. The thumb position is fixed with velcro straps to prevent effects of thumb posture on the flexion moment arm of the first dorsal interosseous (FDI) muscle during contractions (Hudson et al. (2009)).

The index finger is attached to a three-axis load cell (ME Messsysteme K3D40, nominal force 50 N) via a 3D-printed ring-shaped attachment positioned just proximal to the proximal interphalangeal (PIP) joint (Fig. 1A). To immobilize the wrist and eliminate force generation via wrist flexion or extension, the hand is secured on both the palmar and dorsal sides using two adjustable restraints. These restraints can be extended or retracted for fixation, release, and adaptation to different hand sizes.

Figure 1:

The load celĺs x and y outputs are connected to an EMG amplifier system with BNC connectors. Real-time visual feedback of both force amplitude and contraction direction is provided via a custom recording interface, displaying the force vector as a dot within a circular target similar to a clock face. The distance from the center indicates force magnitude either as absolute force in N or as a percentage of maximum force. The angular position represents the contraction direction, with 12 o’clock corresponding to abduction, 9 o’clock to flexion, 6 o’clock to adduction, and 3 o’clock to extension. This mapping matched the actual orientation of the participant’s hand in the setup and can be mirrored for measurements of the left hand.

EMG Setup

We recorded High-Density EMG signals from the FDI muscle with multi-channel intramuscular electrode arrays designed for acute recordings (Fig. 1A) (Muceli et al. (2015); Negro et al. (2016)). The arrays contain 16 platinum-coated electrodes on a polyimide thin-film structure linearly distributed along a distance of 15 mm with an interelectrode distance of 1 mm. The electrode arrays are implanted using a 25G cannula that is removed after the insertion, with only the thin film containing the electrodes remaining in the muscle for the duration of the experimental session. In each participant, one electrode array was implanted in the FDI muscle of the right hand by a medical professional. Prior to insertion, the skin was cleaned, disinfected and the insertion point marked with a skin marker. The insertion was performed in an oblique direction from the distal third of the muscle towards its proximal third. To guarantee safe implantation, avoiding blood vessels and nerves and to ensure the placement of all 16 electrodes within muscle tissue, a desktop ultrasound device (Telemed ArtUs EXT-2H with L15-7H40-A5 probe) was used to guide the implantation in real-time (see Fig. 1A). EMG signals were recorded using a multi-channel amplifier (Quattrocento, OT Bioelettronica, Turin, Italy) at a sampling rate of 10,240 Hz with 16-bit resolution in a monopolar configuration, with the reference electrode placed on the wrist.

Experimental Protocol

Participants were in a seated position with their arm and hand placed in the dynamometer setup as described above. As hand dominance is linked to differences in the common synaptic input to upper-limb MUs (Maillet et al. (2022)), all participants used their dominant hand. The participants performed a series of isometric contractions at constant direction of force application in index finger abduction and flexion. As the DR to MUs can have a strong influence on the estimation of coherence when comparing between conditions (Negro and Farina (2012); Negro et al. (2009); Farina et al. (2014); De la Rocha et al. (2007)), we controlled for comparable net neural drive to the FDI muscle in both abduction and flexion through a dedicated calibration protocol prior to the study protocols. The calibration protocol consisted of isometric ramp contractions at force levels of 2, 3, 4, 5, 6, 7, and 8 N in each contraction direction. For each participant, linear regression was used to determine the force levels in index finger flexion that produced the same global FDI EMG amplitude as in abduction for contractions corresponding to 10% and 20% of the participant‘s maximum force in abduction. EMG amplitude was calculated as the mean root mean square (RMS) value during the plateau phase of the contractions averaged across the 16 iEMG channels. This approach of participant-specific target force selection was chosen because there is no evidence that the FDI receives similar neural drive during abduction and flexion at equivalent percentages of their respective maximum forces across individuals. In pilot data collected with the same dynamometer setup, we observed a wide spread of FDI EMG activity levels during flexion relative to abduction in contractions at matched relative force levels, as well as at matched absolute force levels. More details on the corresponding methods and results can be found in Supplementary figure S2. Across participants, the calibration yielded abduction forces between 2.0 and 3.0 N and corresponding flexion forces between 1.8 and 6.2 N for the 10% maximum abduction force level to be used in the following study protocol.

The study protocol consisted of two separate parts to allow for optimized protocols for each part of the data analysis. Recruitment protocol: Participants performed a series of six isometric ramp contractions for each force level in abduction and flexion with plateau phases of 3s per contraction. Each combination of force level and contraction direction was separately recorded, from here on referred to as “trial”. In both contraction directions, concentric and eccentric phases of the ramps had durations of 5 and 10 seconds for the two force levels respectively, to ensure for equal FDI EMG RMS/s increase and decrease across force levels and direction. Neural input protocol: Participants performed two isometric ramp contractions for each force level in abduction and flexion with plateau durations of 20s per contraction and concentric and eccentric phases of 2.5s and 5s for low and high force levels. In both protocols, each trial was repeated once, resulting in two recordings with identical parameters. During data analysis only the repetition with the better initial decomposition result and sufficient force tracking accuracy was used for each trial.

The 2D real-time force feedback visualizing the force amplitude and exact contraction direction proved to be challenging for first-time users. Consequently, each participant underwent a familiarization session in the days prior to their experimental session, where they practiced isometric ramp contractions at different force amplitudes and contraction directions until sufficient precision in amplitude and direction was achieved. An initial analysis of the study data confirmed high directional accuracy of the applied force across participants during both protocols (Supplementary figure S1B).

Signal Processing

Motor Unit Identification & Tracking

The monopolar EMG signals were bandpass filtered between 100 and 4400 Hz and decomposed into single MU activity through convolutive blind source separation (Holobar and Zazula (2007); Negro et al. (2016)) using the MUedit tool (Avrillon et al. (2024a)) with an extension factor of 40 and a silhouette threshold of 0.9. Subsequent visual inspection and manual editing was performed to identify and correct missed or falsely identified discharges (Del Vecchio et al. (2020)). For the data obtained from the recruitment pro-tocol, particular attention was given to accurate identification of the first firing in each contraction, as this is critical for determining the corresponding recruitment threshold. If the first firing of a MU for a partic-ular contraction could not be identified reliably, all identified firings in that contraction were consequently discarded.

As all contraction levels and directions were decomposed separately, MUs were tracked across force levels within the same contraction direction by concatenating the decomposed signals and application of the MU filters obtained from the decomposition from the lower onto the higher contraction level signals, and vice versa. This approach was described in previous studies with high-density surface EMG data (Frančič and Holobar (2021); Avrillon et al. (2024b)) and proved highly effective in our data (Fig. 1B). The likely reason is that intramuscular EMG provides higher temporal and spatial resolution with less low-pass filtering of the signals by volume conduction through muscle and other tissue. This increases the uniqueness and therefore separability of the action potential shapes compared to MUs data obtained from surface EMG signals. After removal of resulting duplicate MUs in cases where the same MU was originally identified in both contraction levels, this usually resulted in the identification of additional MUs in the higher contraction levels. This reflects the fact that the identification of low-threshold MUs is facilitated at lower contraction levels where the reduced overall neural drive results in fewer simultaneously active MUs and, consequently, less signal superposition compared to higher contraction levels. Following MU tracking across force levels, the same procedure was applied between contraction directions by applying the MU filters from the abduction to the flexion data and vice versa. The entire process was performed separately for the recruitment protocol and the coherence protocol data.

Motor Unit Recruitment Analysis

For the analysis of RTs and the consistency of RO between trials and contraction directions, data from the recruitment sub-protocol were used. Prior to analysis, intramuscular EMG data was band-pass filtered between 100 and 4400 Hz. Notably, Removal of outlier channels was not necessary, as we did not observe any faulty or noisy signals in the data recorded in this study. EMG amplitude was quantified as the root-mean-square (RMS) calculated in windows of 400 ms with 300 ms overlap, averaged across all 16 channels of the thin-film electrode arrays. Recruitment thresholds were quantified in global EMG amplitude in mV, as this provides a more direct measure of net neural drive than absolute force (N) or relative force (% of maximum force) when considering different contraction directions. The rationale for this choice was twofold: first, motor unit recruitment is driven by increases in neural drive, with force changes representing a downstream consequence; and second, MUs recruited at the same net neural drive in both contraction directions would not return two equal RT values if RT is quantified in absolute or relative force levels, as we described above.

For each recorded trial consisting of six contractions, MU RTs were determined for each individual con-traction ramp. The RT for a given MU and contraction was defined as the mean RMS value of the EMG signal within a 50 ms window centered on the MU’s first firing in that contraction. RT values with a z-score greater than 2 relative to the distribution of thresholds for that MU in that trial were discarded as outliers, as these likely reflected missed or incorrectly identified firings during manual editing, or large deviations in the angle of the performed contraction direction, that were observed in few contractions. Only MUs with identified and accepted recruitment thresholds in at least three of the six ramps were retained for further analysis. The individual contractiońs RTs for each MU were averaged and the resulting values sorted to obtain the recruitment order of all retained MUs for that trial.

The consistency of RTs and ROs were then compared between trials of the same contraction direction, as well as between contraction directions using Pearson correlation and statistical methods described in the Statistical Analysis section. We quantified the amount of relative RT change across trials by normalizing the absolute RT difference between the two trials with respect to trial 1, resulting in the normalized difference in recruitment threshold (ΔRT). For analysis of the relation of MU size and ΔRT, we used the peak-to-peak amplitude of the motor unit action potentials (MUAPs) as an estimate of MU size and grouped all MUs into MUs below or above the median MU size of the respective participants identified MUs. We chose MUAP amplitude over RT to estimate the MU size, as we expect RTs values to change between contraction directions, while MUAP amplitudes would stay consistent. To mitigate the effect of distance of MU muscle fibers to the recording electrodes on the observed MU amplitude, we averaged the largest five channelś peak-to-peak values for each MU.

To quantify the number of MUs that exhibited a change in recruitment between contraction directions, we first quantified the test-retest variability of RT and RO between two trials of the same direction by comput-ing the coefficient of repeatability according to Bland and Altman (1986). The coefficient of repeatability represents the 95% limits of agreement for repeated measurements. Because the RT and RO variability did not systematically depend on RT magnitude, we computed pooled coefficient of repeatabilitys (CRs) across the MUs of each participant. We then considered a MU to have changed its recruitment, if both its nor-malized ΔRT and change in RO position were above the respective participant-specific CRs values. This approach ensures that only changes greater than the expected measurement and physiological variability were interpreted as true alterations in recruitment behavior.

MU DRs were computed by averaging the median instantaneous discharge rate (IDR) during the steady plateau phase across all contractions, separately for abduction and flexion.

Motor Unit Coherence Analysis

To estimate the strength and possible sources of synaptic input to the MU pool in different contraction di-rections, intramuscular coherence (IMC) was computed on the motor unit data obtained from the coherence protocols. Only the steady-force plateau phases of the contractions were considered. As described above, MUs were tracked across contraction levels, resulting in concatenated files that include four contractions (two per contraction level). To correct for inconsistencies of the decomposition algorithm in detecting the exact discharge times of MUs, we realigned the MUs of each participant as described in Ibáñez et al. (2021).

As the strength of common input is affected by the number of MU firings included in the cumulative spike trains (CSTs) (Farina et al. (2014)), which is in turn a consequence of the MU DRs given equal contraction durations, we ensured matching MU DRs across contraction directions by participant and direction-specific selection of contraction intensities (see methods above). We further included only MUs that were consis-tently firing throughout a minimum of two of the four contractions, and MUs with a z-score < 2 with respect to all active MUs of each file.

Coherence was computed on the CSTs of equally sized groups of MUs within each contraction direction. As increasing numbers of MUs per group results in higher coherence values, the number of MUs per group was half the total number of detected MUs in the contraction direction of less identified MUs to ensure comparability between directions. MUs were randomly assigned to the two groups with 100 repetitions per coherence computation. Coherence was calculated using MATLAB 2024b (MathWorks Inc.) using the function sp2a2 R2 mt within the Neurospec 2.11 toolbox (Halliday (2015)) in segments of length 1.6s and multi-tapers. We further transformed the coherence values into z-scores using the Fisher transform, scaled by the number of segments in the coherence computation (Laine et al. (2015)). For comparison of the z-coherence of the MU pools between the contraction direction, the area-under-curve (AuC) was computed for each participant in both contraction direction within the delta (1-5 Hz), alpha (5-13 Hz) and beta band (13-30 Hz). Only z-coherence values greater than the bias (mean z-score value between 250 and 500 Hz) were included in the AuC calculation.

Statistical Analysis

We tested for a significant difference in MU recruitment threshold changes compared between trials of the same contraction direction and between abduction and flexion trials. We computed the median normalized absolute ΔRT value within MUs of each participant, once between trials of the same contraction direction (abduction - abduction, flexion - flexion) and between trials of different contraction directions (abduction - flexion). Normality of the paired differences for each participant’s resulting values was confirmed using Shapiro-Wilk test (p=0.72). Significance was then tested using a paired t-test. Effect sizes were expressed as Cohen’s d for paired samples, and 95% confidence intervals of the mean difference are reported.

To check for the relation between MU RT and DR and potential differences of such relation between con-traction directions, R² values of MU RT and DR were computed for each participant individually in both abduction and flexion. A Wilcoxon Signed-Rank Test was used to determine if the correlation between MU RT and DR (R²) differed significantly between conditions. To examine the relationship between the shift in RT and the adjustment in DR across contraction directions, we performed a one-way ANCOVA to determine the common within-subject association between ΔRT and difference in discharge rate (ΔDR), accounting for the dependence of observations within participants.

To investigate the influence of MU size on RTs, MUs were grouped into MUs with MUAP amplitudes lower or higher relative to the participant-specific median MUAP amplitude. Recruitment thresholds were then analyzed using a linear mixed-effects model (LMM), with MU size group included as a fixed factor and participant modeled as a random intercept (normalized absolute ΔRT ∼ 1 + MU size group (1 | participant)). We further tested if the normalized RT change from abduction and flexion leaned significantly towards an increase or decrease in RT separately for both MU size groups. Linear mixed-effects models were applied on the normalized, not-absolute (as opposed to the tests above, where we looked at absolute changes in RT) ΔRT values of all MUs across all participants, with participant ID as random effect (normalized ΔRT ∼ 1 + (1 | participant)), separately for the two MU size groups.

Differences in the correlation of MU amplitude with RT in abduction and flexion were checked by computing the Pearson correlation coefficient‘s r for each participant in each direction, and performing a paired t-test on the r-values after application of Fisheŕs z-transform. Effect sizes are expressed as Cohen’s d, and 95% confidence intervals of the mean difference are reported.

The difference in synaptic input to FDI MUs between abduction and flexion was checked via the AuC of the z-coherence within the delta band (1-5), alpha band (5-13 Hz) and the beta band (13-30 Hz). Normality of the AuC values was confirmed using Shapiro-Wilk test (p=0.91) and significance checked using a paired t-test. Effect sizes were expressed as Cohen’s d for paired samples, and 95% confidence intervals of the mean difference are reported.

Where applicable, all statistical analyses were two-sided. A significance threshold of p < 0.05 was applied (* < 0.05, ** < 0.01, *** < 0.001). All statistical analyses were performed in MATLAB 2024b (MathWorks Inc.)

Results

Motor Unit Decomposition

We implanted a 16-channel intramuscular thin-film electrode array in the FDI muscle, allowing us to decode and track the activity of larger numbers of MUs in a hand muscle compared to traditional intramuscular fine-wire or needle electromyography (EMG) or state-of-the-art high-density surface-EMG. For the recruitment protocol, we were able to identify per-participant averages of 23.3 ± 5.5 MUs for index finger abduction and 23.0 ± 7.0 MUs in flexion. Out of these MUs, 19.8 ± 7.0 were identified in both contraction directions. In the coherence protocol, similar average MU numbers of 22.8 ± 5.2 in abduction and 22.7 ± 6.9 in flexion were decomposed, with 17.8 ± 6.1 MUs tracked across both contraction directions. As intended and ensured through the calibration process to determine participant-specific target force amplitudes for each contraction direction (see Methods), MU DRs were similar across abduction and flexion for each participant. In the recruitment protocol, average MU DR was 12.86 ± 1.88 spikes/s during index finger abduction and 13.02 ± 2.00 spikes/s in flexion. In the coherence protocol, values were similar with 13.76 ± 1.91 spikes/s and 13.94 ± 2.53 spikes/s in abduction and flexion respectively (Fig. 5B).

We did not find evidence of separate MU pools or single MUs recruited exclusively in one contraction direction. In few cases, we found clear absence of motor unit activity in one direction and clear presence of activity in the other. However, these absences likely reflect contraction direction-specific recruitment thresholds not being reached rather than true direction-specific activation.

Reliable recruitment order within, but not across contraction directions

Figure 2A illustrates the process of RT extraction and the observed recruitment variability within and across contraction directions. Across all six participants, FDI MUs showed highly consistent RTs and consequently ROs between two trials of the same contraction direction. Figure 2B shows one exemplary participant‘s ROs and RTs of two trials of abduction and flexion plotted on the x-axes (trial 1) and y-axes (trial 2) with the respective Pearson correlation coefficients. The narrow spread around the diagonal indicates minimal changes of RT within the same force direction. We observed MUs mostly kept their position within the RO, with rare deviations of 1-2 positions in the RO across the full dataset. Pearson correlation was high across all participants within the abduction (r_RT = 0.90 ± 0.12, r_RO = 0.91 ± 0.08) and flexion contractions (r_RT = 0.96 ± 0.07, r_RO = 0.95 ± 0.05). By contrast, recruitment between abduction and flexion revealed much less consistent RTs and ROs (Fig. 2C) (r_RT = 0.46 ± 0.38, r_RO = 0.47 ± 0.38). While some MUs still showed similar RT values in both contraction directions, others were recruited at drastically different amplitudes of net neural drive, consequently also dramatically altering their position in the RO. This behavior was observed in either direction of RT change, meaning some MUs were recruited earlier in abduction compared to flexion, while others showed the opposite. Figure 2A illustrates the reversal of RO between two exemplary MUs. During index finger abduction, MU#12 was consistently recruited after MU#9, while in flexion MU#9 was recruited after MU#12.

Figure 2:

We quantified the relative change in recruitment threshold across contraction directions by computing the normalized absolute ΔRT. Across participants, absolute ΔRT values were significantly increased from within a contraction direction versus across contraction directions as quantified with the average of all partic-ipant‘s median ΔRTs values (average median ΔRT within direction = 13.5 %, average median ΔRT across directions = 46.7 %. Paired t-test: t(5) = 6.16, p = 0.0016. Cohen’s d = 2.51, 95% CI [0.19, 0.47]) (Fig. 2D). The proportion of MUs that exhibited a clear change in recruitment was determined by employing the coefficient of repeatability of the normalized ΔRT and the change in RO position as criteria for recruitment change (see Methods). Within the same functional task (i.e., identical contraction direction), only 5.0% of MUs across all participants showed a change in recruitment, whereas 45.3% changed their recruitment across contraction directions.

Coordinated Modulation of Recruitment and Discharge Rate

In accordance with the onion-skin principle, we consistently found the MU DRs to be inversely correlated to the specific MUs RT across both contraction directions (Fig. 3A, B), meaning MUs with higher RTs had lower steady-state DRs than lower-threshold MUs. Importantly, the strength of RT-DR correlation did not depend on contraction direction, as we did not find any significantly higher correlation of RT and DR in either contraction direction with average R² values of 0.44 ± 0.18 in abduction and 0.48 ± 0.10 in flexion at the lower force level, and 0.22 ± 0.10 and 0.24 ± 0.24 at the higher force level (Wilcoxon signed-rank

Figure 3:

The preservation of a similar RT-DR correlation, despite the drastic observed changes in RT across contrac-tion directions (see above), was mechanistically explained by corresponding, compensatory adaptations in DR. By correlating the difference in RT to the difference in DR between abduction and flexion, we found that MUs that reduced their RT from abduction to flexion showed increased DR and vice versa. This behav-ior was observed across all participants and both force levels (Fig. 3C, D) with strong statistical significance (mean R² F1= 0.72 ± 0.14, mean R² F2= 0.84 ± 0.11, repeated-measures correlation, F1: p¡0.001, 95% CI [-4.86, 6.93], F2: p¡0.001, 95% CI [-3.91, 5.96]).

Changes in Recruitment do not relate to Motor Unit size

We further examined whether estimated MU size was related to changes in RT across contraction directions. We computed the MUAP peak-to-peak amplitude as measure of MU size. Within each participant, MUs were grouped into those below and above the participant‘s median MU size. A linear mixed model with MU size group as a fixed effect and participant as random effect revealed no statistically significant influence of MU size group on normalized absolute ΔRT across contraction directions (p = 0.22, estimate = 0.15, 95% CI [0.09, 0.38], F(1, 118) = 1.51) (see Fig. 4A).

Figure 4:

Figure 5:

As there was a significant absolute change of ΔRT between contraction directions, we checked if MUs on average significantly increased or decreased their RT from index finger abduction to flexion, and if this might be influenced by MU size. We computed the normalized ΔRT values without taking the absolute of the values to retain the information of increase or decrease of RT. We fitted linear mixed-effect models with the normalized, not-absolute ΔRT values as fixed effect and participant as random effect. For both the small and big MUs, the mean ΔRT was significantly different from 0 as tested with linear mixed models (Small MUs: p = 0.034, estimate = 0.29, 95% CI [0.02, 0.56], F(1, 51) = 4.74. Large MUs: p = 0.002, estimate = 0.19, 95% CI [0.07, 0.32], F(1, 67) = 10.32) with an average increase of RT from abduction to flexion for both MU size groups. However, this reflects only the average change in RT and it is important to note there were many observations of both small and large FDI MUs that decreased their RT from abduction to flexion (Fig. 4B).

Ultimately, we compared the strength of correlation of each participant’s MU amplitudes with their cor-responding RT separately for each contraction direction, as we suspected a potentially higher correlation of MU size with RT during index finger abduction compared to flexion. However, on average participants showed similar correlations across contraction directions (Fig. 4C, mean Pearson r in abduction= 0.59, in flexion= 0.51), and a paired t-test on the Pearson correlation values r in abduction and flexion across partici-pants showed no significant effect of contraction direction on strength of correlation (t(5) = 1.09, p = 0.325. Cohen’s d = 0.45), 95% CI of difference [-0.19, 0.46]).

FDI motor units show stronger beta band coherence in abduction than flexion

We estimated the neural drive to the motor neuron pool of the FDI muscle during index finger abduction and flexion by calculating IMC. For each participant, FDI MUs were randomly divided into two subgroups and combined into CSTs before coherence analysis, with 100 random permutations of group assignment (Fig. 5A).

To avoid confounding effects of differing motor unit DRs across contraction directions, we employed a dedicated calibration protocol to select participant-and direction-specific target force levels that matched global EMG activity between abduction and flexion (see Methods). As a result, MU firing rates were highly similar across directions in all participants (Fig. 5B). This step was crucial, since mismatched firing rates at identical contraction durations would otherwise result in unequal numbers of MU discharges contributing to the coherence estimate, potentially biasing the results (Negro and Farina (2012); Negro et al. (2009); Farina et al. (2014); De la Rocha et al. (2007)).

Neural drive was then compared between contraction directions using z-transformed coherence values (z-coherence). We analyzed IMC separately within the delta (1-5 Hz), alpha (5-13 Hz), and beta band (13–30 Hz). Figure 5C shows the intramuscular z-coherence during index finger abduction and flexion averaged across all participants, with noticeably greater coherence values in the beta band from 13-30 Hz during abduction. Across participants, beta-band z-coherence was significantly lower during flexion than abduction (t(5) = 4.59, p = 0.006, Cohen’s d = 1.87, 95% CI [21.79, 77.24]). By contrast, no significant difference in coherence between contraction directions was found in the delta band (Fig. 5D, t(5) = 0.30, p = 0.776, Cohen’s d = 0.12, 95% CI [10.70, 13.54]) and alpha band (t(5) = 1.64, p = 0.162, Cohen’s d = 0.56, 95% CI [9.85, 44.47]).

Discussion

The primary finding of this study is the significant flexibility in the distribution of neural drive to the MUs of the FDI muscle depending on the functional task. In the MU activity extracted from HD-iEMG in six healthy participants, we found that 45.3% of MUs substantially changed their recruitment thresholds (RT) and re-cruitment order (RO) when switching between abduction and flexion (RT and RO above the participant‘s coefficient of repeatability). This variability in recruitment was far above natural variability in MU discharge timings or recruitment thresholds across multiple contractions, as shown by the highly consistent RTs and ROs that we observed between contractions in the same contraction direction with only 5.0% of MUs above the CRs (average ΔRT across directions = 46.7 %, within directions = 13.5 %) (Figs. 2B,C,D). These findings align with previous observations in the FDI (Desmedt and Godaux (1981); Thomas et al. (1986)) and demonstrate that MU recruitment cannot be fully explained by a rigid application of Henneman’s size principle across different functional tasks. Instead, the CNS appears capable of adjusting the distribution of neural drive independent of the motor neuron size, likely through specific excitatory or inhibitory synaptic currents.

Despite the observed variability in RTs, the fundamental control properties of the identified MU pool re-mained stable. In both index finger abduction and flexion, we found statistically similar negative correla-tions between RT and DR, where later recruited MUs exhibited lower DRs than those recruited earlier (Fig. 3A,B). This hierarchical organization has been described as the onion-skin scheme and has been reported in various conditions at low-to-medium force levels (De Luca et al. (1982); De Luca and Hostage (2010); Hu et al. (2014); Tanji and Kato (1973)). The RT-DR relation was weaker at the higher force level (Fig. 3B), which is in agreement with previous studies and may be caused by the saturation of DRs of low-threshold MUs (De Luca and Hostage (2010); Hu et al. (2014)). Notably, MUs that shifted their RT from abduction to flexion showed a corresponding opposite modulation of their steady-state DR (Fig. 3C,D), preserving the RT-DR relation even though large changes in RT were present. This suggests that while the order of recruitment might be task-dependent, the relative organization of the MU pooĺs firing behavior remains consistent. This indicates the likely reason for the apparent onion-skin organization of MU firing activity to be common synaptic input to the observed MUs, rather than being dictated solely by intrinsic motoneuron properties. Furthermore, given that the inverse RT-DR relationship was preserved across functional tasks that involved contractions in distinct contraction directions, it is plausible that proprioceptive feedback plays a significant role in regulating the organization of the MU pool under these conditions.

Interestingly, Desmedt and Godaux (1981) reported mainly smaller FDI MUs being responsible for the observed reversals in recruitment order as they increase their RTs from abduction to flexion. In the large populations of concurrently active motor units we could identify in this study, we did not find estimated MU size to have any effect on the magnitude of change in RT across contraction directions. When grouped into MUs below and above the participant‘s median MUAP peak-to-peak amplitude, there was no statistically significant effect of MU size group on the normalized absolute ΔRT (Fig. 4A). We did find that small MUs generally have a higher RT in flexion compared to abduction. However, the significance was rather weak (p = 0.034, see results and Fig. 4B) and this observation much less pronounced as suggested. We argue this discrepancy is methodologically rooted in the fact that Desmedt and Godaux (1981) reported RT in absolute forces in Newton. However, for the same absolute force, the net neural drive (EMG amplitude) is generally lower in abduction than in flexion with high inter-subject variability (see Supplementary Fig. S2B). Therefore MUs recruited at same levels of neural drive will artificially appear to have higher RTs in flexion if RT is measured in absolute force. Consequently, across different functional tasks of a muscle, RT of MUs should be reported as a metric measured in EMG amplitude rather than absolute or relative force.

The HD-iEMG electrode arrays used in this study enable the concurrent sampling of large numbers of MUs, allowing for correlation-based measures in the frequency domain (Negro and Farina (2012)), specifically intramuscular coherence to explore the neural origins of the redistribution of net excitatory drive that resulted in the observed recruitment and rate-coding variability. We hypothesized that the FDI might receive a unique dominant input from cortical and spinal pathways when it acts as a prime mover (abduction). On the other hand, when the FDI is functioning as a synergist during index finger flexion, there could be an attenuation of high frequency (cortical) oscillatory activity due to the fact that spinal interneurons might suppress neural input to the FDI MUs related to abduction. Within the delta band, which corresponds to low-frequency correlation from 1-5 Hz and reflects the slow, common synaptic input that drives force fluctuations and is directly related to muscle force production and function (Farmer et al. (1993); De Luca and Erim (1994); Negro et al. (2009)), no significant difference in strength of coherence was found. Similarly, the alpha band between 5 - 13 Hz, related to sensory afferents, muscle spindle feedback coupling and includes the physiological tremor frequency (Christakos et al. (2006); Farina et al. (2017); Laine and Valero-Cuevas (2017)), did not show significant difference.

We did find a significantly lower FDI MU coherence during index finger flexion compared to abduction in the beta band from 13 - 30 Hz (Fig. 5C,D). This difference in beta coherence observed may reflect a reduced corticospinal input to FDI MUs during index finger flexion (Conway et al. (1995); Salenius et al. (1997); Baker et al. (1999); Ibáñez et al. (2021)), where the FDI acts as a synergist muscle as opposed to a prime mover as during abduction. We speculate such a suppressed corticospinal input may be rooted in a shared neural input with the MUs of the long flexor muscles in the forearm and higher involvement of spinal circuits during the synergistic flexion task. However, we acknowledge that differences in intramuscular coherence in the beta range may be caused by other factors than variations in corticospinal drive, such as increased synaptic noise or changes in the synchronization of MUs populations. To further clarify the neural determinants of the observed recruitment variability, in particular within the beta frequency range, future studies should assess corticomuscular coherence between FDI and flexor digitorum motor units using concurrent EMG and EEG recordings with similar experimental protocols that involve distinct functional tasks.

On a side note, we found no evidence of distinct “task-specific” motor unit pools in the FDI which is consistent with previous reports (Thomas et al. (1986)). In the few instances where a MU was identified in only one contraction direction, it is most likely that the direction-specific recruitment threshold was simply not reached in the respective other direction, since we observed large changes in RTs between the directions (mean RT change across contraction directions of 46.7%).

Our experimental setup utilized a 2D dynamometer with real-time feedback to strictly control contraction vectors, eliminating voluntary or involuntary deviation in the applied force vector as a confounding factor (Supplementary Fig. S1B). However, certain limitations exist. First, we acknowledge that MUAP amplitude is not necessarily an accurate measure of MU size as distance of the MUs from the electrodes as well as am-plitude cancellation of action potentials reduces MUAP size (Keenan et al. (2006)). However, as RT changed across tasks, it could not serve as a stable marker of MU size. We further did not observe a significantly higher correlation of MUAP amplitude to RT in either contraction direction (Fig. 4C), which would indicate the RTs in that direction could serve as a reliable measure of MU size. Second, the contractions were per-formed at relatively low forces (10–25% MVC). While this range captures a representative distribution of the FDI motor pool, where the upper limit of recruitment range is reported to be even below 50% maximum abduction force, with most MUs even recruited below 30% (Duchateau and Enoka (2022); De Luca et al. (1982)), we cannot extrapolate these findings to the highest-force conditions.

In summary, the present study provides further insights into the distribution of net neural drive within a pool of MUs across different functional tasks. Using current state-of-the-art HD-iEMG technology, we observed significant changes in recruitment thresholds in the FDI muscle in index finger abduction and flexion with corresponding adaptations in rate coding, and found no influence of MU size on the probability of significant change in recruitment. We further quantified the intramuscular coherence within the observed pool of MUs and found significantly decreased beta coherence during index finger flexion compared to abduction, suggesting a suppression of corticospinal input to the FDI during flexion, where the muscle acts as a synergist as opposed to the prime mover. This work further provides the basis for future investigations in the role of underlying corticospinal mechanisms and their alterations in disease.

Supplementary figures

Figure S1.

Figure S2.

Acknowledgements

M.O. acknowledges the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) - Project 523352235. J.I. acknowledges the European Research Council (ERC) - Project ECHOES - 0107769. M.R. acknowledges the German Federal Ministry of Research, Technology and Space (ACS iIMMUNE, 01EO2105) and Deutsche Forschungsgemeinschaft (DFG, German Research Founda-tion) - Project SFB 1483–442419336

Data availability

The data acquired during this study can not be made public due to german data privacy regulations and the applying ethics vote.

Additional information

Funding

Deutsche Forschungsgemeinschaft (DFG) (523352235)

• Marius Osswald

EC | European Research Council (ERC) (0107769)

• Jaime Ibanez

Bundesministerium für Forschung, Technologie und Raumfahrt (BMBF) (01EO2105)

• Martin Regensburger

Deutsche Forschungsgemeinschaft (DFG) (1483-442419336)

• Martin Regensburger