Characterization of Motor Unit Firing and Twitch Properties for Decoding Musculoskeletal Force in the Human Ankle Joint In Vivo

Understanding how motor units (MUs) contribute to skeletal mechanical force is crucial for unraveling the underlying mechanism of human movement. Alterations in MU firing, contractile and force-generating properties emerge in response to physical training, aging or injury. However, how changes in MU firing and twitch properties dictate skeletal muscle force generation in healthy and impaired individuals remains an open question. In this work, we present a MU-specific approach to identify firing and twitch properties of MU samples and employ them to decode musculoskeletal function in vivo. First, MU firing events were decomposed offline from high-density electromyography (HD-EMG) of six lower leg muscles involved in ankle plantar-dorsi flexion. We characterized their twitch responses based on the statistical distributions of their firing properties and employed them to compute MU-specific activation dynamics. Subsequently, we decoded ankle joint moments by linking our framework to a subject-specific musculoskeletal model. We validated our approach at different ankle positions and levels of activation and compared it with traditional EMG-driven models. Our proposed MU-specific formulation achieves higher generalization across conditions than the EMG-driven models, with significantly lower coefficients of variation in torque predictions. Furthermore, our approach shows distinct neural strategies across a large repertoire of contractile conditions in different muscles. Our proposed approach may open new avenues for characterizing the relationship between MU firing and twitch properties and their influence on force capacity. This can facilitate the development of targeted rehabilitation strategies tailored to individuals with specific neuromuscular conditions.


I. INTRODUCTION
H UMAN movement is orchestrated by the coordinated interaction of the central nervous system (CNS) and the musculoskeletal system.The smallest incremental unit that the CNS uses to control musculoskeletal force is the motor unit (MU) [1], which consists of an alpha motor neuron (MN) and the muscle fibers it innervates.Despite extensive research on both alpha MNs and the innervated skeletal fibers, there are still gaps in the understanding of how MN function is transduced into skeletal mechanical force [2].An important element hampering progress is the limited understanding of how differences in MU's firing and twitch contractile properties impact musculoskeletal force generation in intact humans in vivo.Such differences are demarcated by the continuum spectrum of MU phenotypes that directly dictate MU firing rate, contractile speed, and rate of force generation [3].In this context, the impact of alterations in MU properties on skeletal muscle force-generation capacity remains unclear.Transitions in MU phenotypes, i.e., MU remodeling, can arise from different factors such as aging [4], injury [5], and exercise training [4].For example, transitions towards smaller, slower MUs have been reported due to the loss of larger MUs in elderly humans [4], resulting in a decline in functional capacity.Skeletal muscle training can help attenuate these age-related alterations by facilitating reinnervation [4], increasing MU size (hypertrophy) [4], decreasing recruitment thresholds [6] and increasing discharge rates [6].Thus, understanding how different MU types and their associated properties influence motor function is central for designing personalized training programs and neuro-rehabilitation treatments.
An established way to non-invasively study MUs in humans in vivo is based on the use of high-density electromyograms (HD-EMGs) in combination with blind-source separation [7], [8], which can capture the concurrent activity of multiple MUs.Interfacing with MUs in vivo provides valuable insights into healthy [9] and pathological [10], [11] neural patterns under different conditions (e.g.nerve or spinal electrical stimulation [12], [13]).Moreover, this enables the longitudinal tracking of MU adaptations resulting from exercise training (e.g., endurance [14] and strength [6]).Nonetheless, relying solely on HD-EMG decomposition is insufficient to establish causal associations between MU properties and musculoskeletal force.
In contrast, mechanistic, numerical models of the composite neuro-musculoskeletal system could potentially reveal causalities between measured neural activity (e.g., spinal MN firing patterns) and resulting mechanical force generated in innervated skeletal muscle-tendon units, as well as moments in biological joints.This is critical for studying and testing specific hypotheses underlying the neuro-mechanics of skeletal muscle contraction and multi-joint motor control [15], [16].
In particular, EMG-driven musculoskeletal modeling frameworks [18], [20] are becoming widely-used computational methods for simulating the mechanics of the human musculoskeletal system as controlled by neural surrogates, i.e.EMG-derived neural activations.However, global EMGs can be constrained by intrinsic limitations, including action potential amplitude cancellation [21] and the volume conductor effect [22], which limit their association with the neural drive to muscles [23] and subsequent force output [24].Moreover, established models driven by global EMGs cannot explain the contribution of individual MUs.Therefore, these approaches are unable to explain how the CNS controls muscle force by modulating the number of MUs (i.e., recruitment coding) and the rate at which each of them discharges an activation signal (i.e., rate coding) in a person-and task-specific way.This limits our ability to understand how individual MUs, which produce millinewton-level forces, collectively generate thousands of Newtons in a fraction of a second.
In this context, we propose a novel method to estimate muscle force and joint moments by explicitly considering both the firing and twitch properties of individual MUs experimentally decoded in vivo.Although previous studies have presented musculoskeletal models driven by spike trains, these models were either derived from animal data or dissociated from in vivo human MU firing activity [15].More recently, high-density EMG-driven modeling formulations were proposed that utilized MU spike trains (i.e., firing function) as the main drive of the model but disregarded the twitch dynamics of each MU [20], [25], [26].These formulations would not enable investigating how alteration in individual MU firing and contractile behavior (e.g., in response to aging, injury, or training) affects skeletal muscle force generation, thereby hampering the development of personalized training procedures.Thus, by ignoring the mechanical differences between MUs and assuming they have the same force production capabilities, these models may overlook important aspects of neuro-motor control, thereby limiting our ability to accurately predict and understand movement and pathological patterns.
In this study, we address the limitations of state-of-theart methodologies in dictating how MU firing and twitch properties affect force generation in a person-and multimuscle-specific manner.First, we employ MU statistical distributions of firing properties to estimate individual MU twitch responses.Second, we use the twitch responses to generate MU-specific activation dynamics.Third, we integrate our proposed framework with person-specific musculoskeletal models [18], [19] to predict joint moments.By these means, our methodology may open new directions to study longitudinal changes in MU phenotype due to aging, injury, or training.
To determine hip, knee, and ankle joint centers of rotation, we recorded motion capture data and ground reaction forces synchronously using a seven-camera system (Qualisys, Göteborg, Sweden, 256 Hz) and a force plate (Bertec Co., Columbus, OH, USA, 2048 Hz).We placed 18 retro-reflective markers on each participant's pelvis and right lower extremity.Data were recorded during one static anatomical pose and a set of functional trials and low-pass filtered with a fourth-order Butterworth filter (cut-off frequency: 6 Hz).
We recorded ankle angular moments and positions using a dynamometer (M3, Biodex Medical Systems Inc., Shirley, NY, USA), synchronized with HD-EMGs from a 256-channel amplifier (EMG-USB2, OT Bioelettronica, Torino, Italy, sampled at 2048 Hz) during isometric plantar-dorsi flexion tasks.Dynamometer data were filtered using a fourth-order lowpass Butterworth filter with a cut-off frequency of 2 Hz.HD-EMG data were measured from the right lower leg muscles (Fig. 1.a): tibialis anterior (TA), soleus (SOL), medial and lateral gastrocnemius (GM and GL), peroneus longus (PL) and peroneus tertius (PT).Prior to electrode placement, the skin was shaved and lightly abraded.Subsequently, the electrode grids were placed on the skin using a 1-mm thick double-adhesive foam and a conductive paste.For the TA, SOL, and GM, 64-channel grids (10 mm inter-electrode distance, 13 rows, and 5 columns) were used.For the GL and peroneus group, 32-channel grids (10 mm inter-electrode distance, 8 rows, and 4 columns) were used.The EMGs of the peroneus group were split into PL (during plantar flexion) and PT (during dorsiflexion).The reference electrode was placed on the right malleolus.The HD-EMG signals were processed to be decomposed into MN firing events (see Section II-D) and for validation against linear envelopes.For decomposition, the HD-EMGs were bandpass filtered with a second-order, zero-lag Butterworth filter (cut-off frequency: 10-500 Hz).To generate the envelopes, HD-EMGs were averaged across all channels per muscle.Then, the resulting signals were high-pass filtered (cut-off frequency: 30 Hz), fully rectified and low-pass filtered (cut-off frequency: 20 HZ) with secondorder, zero-lag Butterworth filters.For each muscle and each subject, the obtained envelopes were normalized with respect to the maximum value obtained from the entire set of trials.

B. Experimental Protocol
Each isometric plantar-dorsi flexion task comprised four repetitions of symmetric triangular moment profiles targeting a percentage of maximum voluntary contraction (%MVC at a known ankle position.The ankle positions were the same for all subjects: anatomical, 10 • dorsi-, and 10 • plantarflexed.Regarding the MVC condition, three subjects followed moment profiles from rest to 30, 50, 70, and 90 %, and back  [7], [8].The firing properties are calculated and employed to map MU twitch properties.(c) The estimated twitch properties are used to design MU-specific twitch responses defined as an impulse of a second-order response system [17].The convolution between each j th spike train (u(t)) and its respective twitch response provides the MU-specific activation dynamics (a(t)).(d) Together with the ankle joint angles, the MU-specific activation serves as input to a person-specific musculoskeletal framework [18], [19].RT and DR: recruitment threshold and normalized discharge rate, A p and T c twitch peak amplitude and contraction time, ch: channel, CE: contractile element, SEE: series elastic element, PEE: parallel elastic element, L M o : optimal fiber length, L T s : tendon slack length, F max : max.isometric force.to rest within 4 seconds (i.e., 2 seconds per ramp).These tasks were designed in a way that the slopes of the ramps were variable across MVC targets, i.e., the higher the moment target, the higher the slope.Conversely, the fourth participant Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I CONTRACTION TIMES FOUND ON HUMANS
followed moments with fixed slopes of 10% MVC/s across different MVC conditions: 20, 30, 40, 50, and 60% MVC.Namely, the duration of the ramps varied for this subject to keep the slopes constant across conditions (ankle positions and %MVC).

C. High-Density Electromyography Decomposition
The band-pass filtered HD-EMG signals were visually inspected and poor-quality channels (i.e., exhibiting low signal-to-noise ratio) were excluded from further analysis.We then decomposed the HD-EMGs into MU spike trains using a de-convolutive blind source separation technique [7].Only MU spike trains with mean discharge rates (see section II-D) falling between 3 and 40 Hz were included.We assessed the quality of the decomposition using the pulse-to-noise ratio (PNR) metric [31].For the present study, only MUs with PNR>26 dB were considered.The identified MUs were further inspected by an experienced operator who retained MUs with consistent discharge patterns.Non-physiological discharges, i.e., with unusual inter-spike intervals (<25 ms or >250 ms), were manually edited.

D. From Firing Behavior to Twitch Characteristics
1) Decoding Firing Properties: we estimated the recruitment threshold and normalized mean discharge rate for each MU.Recruitment threshold (RT ) was defined as the mean percentage of MVC around the first firing event (1): where w is a time window of 300 samples (∼150 ms) around the first firing event (t 1 ), and m(t) is the normalized joint moment.The mean discharge rate was defined as the mean inverse difference of the time intervals between consecutive firing events (2): where t is the time of every i th discharge and N is the total number of spikes in a single spike train.The mean discharge rate was normalized to 40 Hz, as we did not observe discharge rates above this value, and previously reported peak rates for human MUs are typically below 50 Hz [2].
2) Decoding Twitch Characteristics: to estimate twitch profiles, we first employed the firing features for discriminating a diverse continuum of MU types (Fig. 1.b).For this purpose, we computed a linear combination of normalized firing properties (mean-centered) through principal component analysis (using singular value decomposition).We obtained a compound feature by projecting the data onto the first principal component.As we sampled MUs in a broad range of force levels (20-90 %MVC), this new compound feature enabled the identification of a spectrum of MU types, ranging from slow (low-threshold) to fast (high-threshold) MUs.
Subsequently, we mapped the compound feature into contractile twitch characteristics based on Hennenman's principle of orderly recruitment [32] and extensively reported matching correlations between MU firing and twitch properties [33], [34], [35].Namely, we employed the compound feature to assign twitch properties from slow (long contraction times and low peak amplitudes) to fast twitch responses (shorter contraction times and higher peak amplitudes).Although highly flexible MU twitch models with multiple twitch parameters have been proposed [16], [36], we previously found that contraction time (or time-to-peak) and peak amplitude of a MU twitch are sufficient to describe the activation dynamics accurately, without compromising computational efficiency [37].We computed a linear regression between minimum and maximum values of the compound feature and the contraction times intervals found in experiments on humans [27], [28], [29], [30] (see Table I).In such manner, we assigned longer contraction times to slow MUs and shorter contraction times to fast MUs in a linearly continuous fashion.Similarly, we computed a linear regression between the limits of the compound feature and the normalized peak amplitudes (from 0 to 1), thereby assigning lower amplitudes to slow MUs and higher amplitudes to fast MUs.
Our proposed transformation from firing to twitch characteristics is summarized by the following vectorial equation (3): where y is either of the twitch properties with its corresponding PCA coefficients (c 1 and c 2 ) and center points (dr 0 and rt 0 ), and its interpolation coefficients ( p 1 and p 2 ) that take into account the ranges found in literature (see Table I).

E. Motor Unit-Specific Activation Dynamics
We employed the contraction times and peak amplitudes to design twitch models for each MU.Each twitch response was computed as an impulse response of a second-order system (4) [17]: where a j is the activation of the j th MU, A c the peak amplitude, T c the contraction time and u the excitation input, i.e., a MU spike train.We discretized the system (as shown by [15]) to filter each spike train with its respective twitch model in an online fashion, thereby obtaining the MU-specific activation dynamics (Fig. 1.c).The sum of every MU contribution within a muscle was defined as muscle activation ( 5): where N R is the total number of recruited MUs at a certain contraction force level.However, we could only access a limited subset of MUs via surface HD-EMG decomposition (i.e., N R → N D , where N D is the number of decomposed MUs).Thus, we accounted for the contribution of the unidentified MUs using a parametric approach, i.e., by multiplying the activation by the ratio between the total number of recruited and decomposed MUs (N R /N D ) as follows ( 6): This step is particularly important to match the increase in activation and force amplitude, despite decomposing approximately the same amount of MUs across MVC levels.We estimated N R as the proportion of recruited MUs (rec) multiplied by the total number of MUs (N T ) in the pool (i.e., N R = N T rec).To describe the proportion of recruited MUs (rec) as a function of % of MVC (i.e., % of force f ), we adapted the exponential equation derived from experimental data of the TA (7) [38]: where m is a shape factor that dictates how fast the pool is recruited.In this study, we assumed the same shape factor (m = 0.045, derived from [38]) for all muscles, as comparable recruitment curves have been found in different muscles (e.g.TA and first dorsal interosseus in [38]).Thus, (6) can be expressed as ( 8): As we integrated this formulation into a subject-specific musculoskeletal framework [18], [19] that takes a normalized activation (from 0 to 1) as input to compute angular moments (Section II-F), we normalized (8) by dividing it by N T (9): We finally further normalized the total activation dynamics by the maximum activation found per subject.

F. Neuro-Musculoskeletal Modeling
We estimated ankle joint plantar-dorsi flexion angular moments using two approaches: MU-specific (Fig. 2) and conventional EMG-driven musculoskeletal models [18], [19], [20].The difference between these models lies in the way activation dynamics were estimated.In the MU-specific framework, muscle activation is computed directly from decoded spike trains (Section II-E), while in the EMG-driven models, it is computed as a linear envelope (Section II-A).The musculotendon kinematics (i.e., musculotendon unit length and moment arms) was computed via a family of B-spline functions dependent on knee and ankle joint angles [19].Musculotendon kinematics and activations were employed as inputs to a Hill-type muscle model (contraction dynamics in Fig. 2), which estimated muscle-tendon force (Fig. 1.d).Net joint moments were computed by multiplying muscule-tendon Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE II DESCRIPTIVE STATISTICS OF FIRING AND TWITCH PROPERTIES
force by moment arms.The calibration tool adjusted musculoskeletal parameters using a simulated annealing algorithm to minimize the mean squared error sum between the predicted and reference moments.The calibrated parameters were tendon slack length, optimal fiber length, and a scaling coefficient of the maximum isometric force.The tendon slack length and optimal fiber length were adjusted within ± 5% of their initial value for fine-tuning muscle-tendon forcelength-velocity relationships.The initial values were linearly scaled and pre-optimized [39].The scaling coefficient of the maximum isometric force varied within ( 1 2 ,2).We calibrated these parameters using a single repetition from the 30% and 50% MVC conditions for the first three subjects, and from the 20% and 40% of MVC conditions for the fourth subject across all ankle positions.

G. Statistical Analysis
Normalized root mean squared error (RMSE) values were calculated to compare the performance in moment prediction between the MU-specific and EMG-driven models.The normalization was done by the reference moment to have comparable RMSE values across different levels of activation.
To check the variability of the predicted moments, the coefficient of variation (CoV) was computed for each model and for the actual moment.We defined CoV as the standard deviation of the detrended moments divided by the mean.We tested statistical differences in RMSE and CoV between MU-specific and EMG-driven predictions using non-parametric Wilcoxon signed-rank tests.Median and interquartile ranges (IQR) were calculated for the firing and twitch properties.Confidence intervals (95%) of the medians were estimated using the percentile bootstrap method (10,000 bootstrap samples).

A. Motor Unit Statistical Distributions
Fig. 4.a and Fig. 5.a show the MU distributions of the firing properties (discharge rates and recruitment thresholds) and the twitch characteristics (contraction times and peak amplitude)  Table III reports the number of identified MUs per muscle levels of activation.Fig. 4.b and 4.c show the discharge rates and recruitment thresholds across % levels of activation for the TA and the SOL, respectively.For the TA, both firing properties exhibit a shift towards higher frequencies and recruitment thresholds as the % level of activation increases.For the SOL, the recruitment thresholds also exhibit a shift as the % level of activation increases, however, a moderate modulation in discharge rate across activations.It is noteworthy to mention that in this study, we observed no modulation of the firing properties across ankle positions..cshow the contraction times and peak amplitudes across % levels of activation for the TA and for the SOL, respectively.For both muscles, the MU distributions exhibit a shift towards lower contraction times and higher peak amplitudes as the level of activation increases.Similar to the firing distributions, we observed no modulation of the twitch properties across ankle positions.
Table IV shows the PCA coefficients, center points (dr 0 and r t 0 interpolation coefficients ( p 1 and p 2 ) of contraction times (T ) for each muscle, i.e., the transformation coefficients from firing to twitch properties.As we normalized the twitch amplitudes, the interpolation coefficients for the peak amplitude (A c ) for all muscles were p 1 = 0.86 and p 2 = 0.44.

B. Moment Prediction
We assessed the ability of our myoelectrically-driven models to predict accurate moments across joint angles (fiber lengths) and multiple sub-maximal contractions.Fig. 6 depicts the mean moment profiles across all conditions for each myoelectrically-driven model and for the measured moment.For each repetition, we evaluated the performance of our framework with normalized RMSEs and CoV.Fig. 7 shows these performance metrics in three ways: a jitter plot, a boxplot, and probability density functions.
In terms of RMSE (Fig. 7.a), the performance of the MU-specific model did not differ significantly from the EMG-driven model ( p > 0.05).Although the median of the MU-specific framework is slightly higher than the EMG-driven approach, the prediction errors of the EMG-driven model displayed higher errors, ranging beyond 100% of normalized RMSE.For our proposed MU-specific approach, a deeper analysis across activations shows that the medians of the low activation trials (≤ 50% MVC) were significantly lower than the higher activation trials.On the other hand, the medians of the EMG-driven models did not display a specific trend.The lowest RMSEs were within 40-50% MVC, followed by the 20-30% and 60-70% MVC with no statistical difference between them, and lastly the 70-90 % MVC group with the highest RMSEs.
Regarding the variability of the predicted moments (Fig. 7.b), the CoVs of the actual moments were significantly smaller than the ones of both models ( p < 0.001).Likewise, the CoVs of the MU-specific model were significantly smaller than the ones of the EMG-driven predictions ( p < 0.001).Across activations, the measured moments displayed higher CoVs in the lowest range of activation (20-30% MVC) and similar CoVs for higher % levels of activation (≥ 40% MVC).The CoVs of the MU-specific model predictions were also the highest for the lowest range of activation (20-30% MVC), however, followed by the 90% MVC group, and then by the 40-50% and 60-70% MVC ranges with no statistical difference between them.The CoVs of the EMG-driven model predictions were the highest for the lowest range of activation (20-30% MVC), followed by the 40-50% MVC range, and then by the 60-70% and 90% MVC groups with no statistical difference between them.

IV. DISCUSSION
We presented a signal-driven model-based framework to simultaneously identify MU firing and twitch properties as well as muscle-tendon force-generating parameters.With this, we reconstructed net moments at the ankle joint by taking into account 6 lower leg muscles in 4 healthy individuals performing a large repertoire of plantar-dorsi flexion contractions.Our approach differs from previous state-of-the-art methodologies, which disregard the twitch dynamics of each specific MU [18], [19], [20], [25], [26].Instead, we formulated a MU-specific model that captures the continuum diversity of MU firing and twitch properties.This approach can open new avenues in the future for tracking firing and twitch adaptations in MUs with subsequent impact on force-generation capacity, due to neuromuscular injuries [5] (e.g., stroke, spinal cord injury), degenerative or autoimmune disorders [40] (e.g., Parkinson's disease, multiple sclerosis), or aging [4] (e.g., sarcopenia).
Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.Fig. 6.Mean moment profiles (dashed lines) and standard deviation (shaded area) across ankle positions (rows), and levels of activation (columns) for each model (red and blue for the high-density electromyography (EMG) and the motor unit (MU)-specific models) and measured moment (yellow).
Importantly, future iterations of this technology can help track training-induced neuromuscular remodeling for personalized interventions.
Our methodology comprised four main elements (Fig. 1).The first was an interface to decode the firing events of MUs across multiple muscles (Fig. 1.a).This allowed the calculation of statistical distributions of the firing properties of the sampled MUs.Second, we established a link between MU firing and twitch properties (Fig. 1.b) to create MU-specific activation dynamics for each muscle (Fig. 1.c).Lastly, we coupled our framework with person-specific neuromusculoskeletal modeling [18], [19] (Fig. 1.d and Fig. 2).
Our proposed analysis of the resulting statistical distributions revealed distinct firing and twitch characteristics across the major muscles of the ankle (Fig. 4.a and Fig. 5.a) and across different levels of sub-maximal contractions (Fig. 4.b-c and Fig. 5.b-c), with consistent behavior across ankle positions.Results showed that the SOL presented a tight concentration towards low discharge rates and high contraction times (Fig. 4.a, Fig. 5.a, and Table II).This effectively reflects the MU composition of the SOL, as it is mainly constituted by slow MU types [41].Conversely, the rest of the muscles spread distributions (Fig. 5).This is well in line with the more heterogeneous compositions of the TA, GM, GL, and PL [41].
The major analyses of our work focused on TA and SOL due to their greater sample size of decoded MUs, and thus the better representation of the MU pool.This could be due to the chosen motor task (plantar-dorsi flexion), as TA and SOL are the primary dorsi and plantar flexors.Although the gastrocnemii can greatly contribute to plantar flexion when the knee is fully extended [42], [43], [44], [45], the selected knee angle (120 • ) decreased their engagement.Future work should test whether a muscle-specific experimental design enables replicating this methodology on the rest of the ankle muscles.
In this study, the SOL and the TA showed different strategies of neuromuscular control across levels of activation (% of MVC).Across levels of activation, distributions showed more discharge rate modulation for TA than for SOL (Fig. 4.b).In addition to the aforementioned difference in MU composition, this suggests dominance of recruitment coding for SOL.As slow-twitch fibers can be tetanized with relatively slow discharge rates [46], the SOL needs to progressively recruit more MUs to achieve contraction strengths close to MVCs [47].Notably, MUs of the SOL exhibited larger modulation of contraction times than discharge rates (Fig. 5.b).This may be due to the influence of the recruitment threshold on the hereby proposed mapping to twitch properties.This highlights the relevance of considering both firing features to characterize MU phenotypes and contractile properties.
Across ankle positions (i.e., muscle lengths), decoded MU populations showed similar distributions in discharge rates.This is in line with recent studies [9], [48] where tracking individual MUs revealed similar discharge rates across different muscle lengths.Likewise, contractile twitch properties of recruited MUs remained unaltered across joint angles for both, TA and SOL.Although this differs from other studies [48], the reason for such conflict may be because our proposed twitch dynamics is modeled at the activation level and serves as input to a forward musculoskeletal model, where changes in muscle length are explicitly accounted for.This can be observed in our predictions where changes in ankle position are reflected in changes in moment amplitudes (Fig. 6).Another reason may be the limited range of ankle positions included in this study.As we studied multiple muscles (including both plantar and dorsi flexors), we chose angles around the anatomical position (−10 • to 10 • ).On the other hand, studies of individual muscles, such as the one by Cudicio et al. [48] for the TA, observed changes at angles around the optimal fiber length (0 • to 40 • ).
Our proposed MU-specific musculoskeletal model enabled accurate moment predictions for a wide range of submaximal voluntary contractions (ranging from 20 to 90 % of MVC) and joint angles (from 10 • of dorsiflexion to 10 • of plantar flexion).While the results of the MU-specific model were comparable to those obtained by EMG-driven models (Fig. 7.a), this approach offered the additional benefit of elucidating how the CNS orderly recruits different types of MUs for force generation, a crucial step towards understanding neural control of force [30], [49].Moreover, the MU-specific model achieved higher generalization than traditional EMG-driven models (i.e., lower RMSE variability across conditions), despite only calibrating it with a single repetition of two % levels of activation.Contrarily, the EMG-driven model's performance varied for each condition (e.g. higher for 40-60% MVC and lower for 30 % MVC, Fig. 6).Notably, the distributions of twitch properties are not integrated into the calibration tool yet, i.e., they may not be optimal.This indicates a greater degree of adaptability of our methodology across different motor tasks, even without optimization of twitch parameters.Future work will focus on optimizing the twitch parameters by adjusting the means and standard deviations of their resulting fitted Gaussian curves.This will improve moment estimation and provide a more tailored set of contractile parameters.
In terms of variability, the MU-specific framework showed significantly lower coefficients of variation (Fig. 7.b) than the EMG-driven models, as observed in the shaded areas of moment profiles depicted in Fig. 6.This is due to an improvement in the activation dynamics by assigning physiological twitch properties to each MU and the elimination of noise components via MU decomposition and selection.Although the model's predictions do not reach the same CoVs as the experimental moments, it represents a significant improvement over the state-of-the-art and may lead to better control schemes of neuro-rehabilitation technologies.
We successfully achieved matching the increase in activation with the increase in force magnitude by accounting for non-identified MUs (i.e., normalization in Section II-D).However, this methodology inherits the intrinsic limitations of MU decomposition, including silent periods in the spike trains.The high-moment profiles in Fig. 6 (see columns >50%MVC) show that, even though our twitch models matched the reference amplitude, they produced a late moment response as we were not able to decompose firing events from the onset of the ramp.This late response is due to the challenge of decomposing smaller MUs that are recruited earlier, in the presence of bigger MUs (with greater action potentials).In these cases, the standard EMG-driven models provided less error in the predictions.Future work will focus on improving the decomposition in high-moment trials and coalescing the present methodology with a subject-specific MN pool modeling framework capable of generating complete in silico MN pools based on the neuronal characteristics of decomposed MUs [50], [51].
Another limitation of our methodology is the assumption of constant MU twitch responses across time.While this assumption ensures computational efficiency and served well for predicting joint moments accurately in our experiments (involving short-duration and slow contractions), it does not fully account for the non-linear nature of tetanic force development [36].In particular, MU potentiation [52] and fatigue [53] play significant roles in dictating progressive changes in twitch responses.Future iterations will incorporate state-of-art formulations (e.g., for potentiation [54] and fatigue [55]) to fine-tune these time-dependent variations and enhance the physiological realism of our models.This will further improve the accuracy of muscle force predictions, especially during long, fast, and repetitive contractions.
Although the present study was conducted offline, our methodology was specifically designed to interface with real-time musculoskeletal models [56], [57] that operate in short time frames (<14 ms), well below the muscle electromechanical delay.With this intention, we integrated computationally efficient formulations to obtain firing properties (Eq.1-2), twitch characteristics (Eq.3), and MU-specific activation Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
dynamics (via discretized twitch models [15]).Future work will further incorporate real-time decomposition techniques (e.g., [58], [59]) for establishing human-machine interfaces that concurrently decode MN firing patterns, contractile properties, muscle forces, and joint moments in real-time.This has multiple implications for clinical implementation since, once calibrated, our model-based approach can operate as a function of MU firing events and joint angles, both of which can be directly derived or measured from wearable sensors.Nonetheless, a broader and more diverse participant cohort is needed before clinical translation, since MU function may vary depending on age [4], pathology [5], [40], motor tasks [60], [61], and level of training [6], [62].Future work will expand the generalizability of our methodology to encompass individuals of varying ages and impairments as well as a wider range of motor tasks.
Notably, we performed our predictions in an open-loop fashion with no corrective adjustments to compensate for errors in our estimated twitch parameters.Although we validated our framework by predicting accurate moments, further validation is needed at the twitch level.Future work will employ imaging techniques, such as ultra-fast ultrasound imaging [63] or microendoscopy [64], to corroborate our estimated properties against identified twitch dynamics of individual MUs.
In conclusion, our proposed MU-specific framework demonstrated robustness in accurately predicting joint moments across different ankle positions and levels of activation.Moreover, this approach allowed a more detailed understanding of the interaction between MUs and the innervated muscles at multiple levels, including neural, muscular, and joint levels.This is valuable information to help track condition-specific neural and mechanical adaptations in the MU pool composition.Consequently, MU-specific models could greatly enhance the design of effective interventions to target specific MU types, thereby optimizing recovery of neuromuscular function following impairment.

Fig. 1 .
Fig. 1.Study overview.(a) Experimental set-up for ankle dorsi-plantar flexion tasks.Ankle joint moments and high-density electromyography (HD-EMG) from the lower leg muscles are recorded.(b) HD-EMGs are decomposed into motor unit (MU) spike trains using a convolutive blind-source separation (BSS) technique[7],[8].The firing properties are calculated and employed to map MU twitch properties.(c) The estimated twitch properties are used to design MU-specific twitch responses defined as an impulse of a second-order response system[17].The convolution between each j th spike train (u(t)) and its respective twitch response provides the MU-specific activation dynamics (a(t)).(d) Together with the ankle joint angles, the MU-specific activation serves as input to a person-specific musculoskeletal framework[18],[19].RT and DR: recruitment threshold and normalized discharge rate, A p and T c twitch peak amplitude and contraction time, ch: channel, CE: contractile element, SEE: series elastic element, PEE: parallel elastic element, L M o : optimal fiber length, L T s : tendon slack length, F max : max.isometric force.

Fig. 2 .
Fig. 2. MU-specific neuro-musculoskeletal modeling framework: each muscle contains N decoded spike trains, each assigned to a twitch response (contraction time, T c , and peak amplitude, A p ) depending on their firing properties.The activation dynamics of a MU is the tetanic summation of its respective individual twitches of a motor unit (MU).The musculotendon kinematics (musculotendon unit length (L MT ) and moment arms (r)) and activation dynamics (a(t)) are inputs to a Hill-type muscle model, which estimates muscle-tendon force (F T ) accounting for the force-length (F M -L), force-velocity (F M -V), and tendon force-strain (F Tϵ), relationships.The joint moment for each muscle m(t) is the multiplication between muscle-tendon force and moment arms.The total joint moment M(t) is the summation of the moments produced by each contributing muscle.The diagonal arrows indicate the parameter tuning for each subject.

Fig. 3 .
Fig. 3. Translation from firing behavior to twitch contraction time (tibialis anterior example).(a) Representation of projection of normalized data onto the first principal component (non-centered for visual purposes).Each point represents a sample of a motor unit (MU).(b) Histogram of data projection normalized by probability density (PD) function estimate (blue line).The red dashed lines correspond to different MU populations (from slow to fast MUs).(c) Firing properties are mapped by linear interpolation with corresponding twitch contraction times (T c ) found in the literature (Table I).(d) Transformed distribution into twitch contraction times (T c ). DR: discharge rate, a.u.: arbitrary unit).

Fig. 4 .
Fig. 4. Firing Properties across muscles (a) and across levels of activation for the tibialis anterior (TA, b) and the soleus (SOL, c).The discharge rates are displayed on the left and the recruitment threshold on the right.GL, GM: lateral and medial gastrocnemius, PT, PL: peroneus tertius and longus, MVC: maximum voluntary contraction, and a.u.: arbitrary unit.

Fig. 5 .
Fig. 5. Twitch Properties across muscles (a) and across levels of activation for the tibialis anterior (TA, b) and the soleus (SOL, c).Contraction times (reversed x-axis) are on the left and peak amplitudes on the right.GL, GM: lateral and medial gastrocnemius, PT, PL: peroneus tertius and longus, MVC: max.voluntary contraction, and a.u.: arbitrary unit.

Fig. 5 .
Fig.5.b and 5.c show the contraction times and peak amplitudes across % levels of activation for the TA and for the SOL, respectively.For both muscles, the MU distributions exhibit a shift towards lower contraction times and higher peak amplitudes as the level of activation increases.Similar to the firing distributions, we observed no modulation of the twitch properties across ankle positions.TableIVshows the PCA coefficients, center points (dr 0 and r t 0 interpolation coefficients ( p 1 and p 2 ) of contraction times (T ) for each muscle, i.e., the transformation coefficients from firing to twitch properties.As we normalized the twitch amplitudes, the interpolation coefficients for the peak amplitude (A c ) for all muscles were p 1 = 0.86 and p 2 = 0.44.

Fig. 7 .
Fig. 7. Normalized root mean squared error (RMSE) values (a) and coefficient of variation (CoV) values (b) for the motor unit (MU)-Twitch models (blue), EMG-driven models (red), and measured torque (yellow).Each point of the jitter plots represents a single repetition, i.e., a dorsi-plantar flexion contraction.The black triangles of the boxplots represent notches, i.e., 95% confidence interval of the median values (black dot).On the right, the probability density functions are displayed.(***) indicates significant differences (p<0.001).

TABLE III NUMBER
OF IDENTIFIED MUS PER MUSCLE ACROSS %MVC across muscles.TableIIsummarizes the descriptive statistics of these properties (median, its 95% CI, and IQRs).

TABLE IV PCA
COEFFICIENTS, PCA CENTER POINTS, AND INTERPOLATION COEFFICIENTS OF CONTRACTION TIMES