Introduction

Human muscles are versatile effectors producing forces that span several orders of magnitude. They allow humans to perform extremely diverse motor tasks with the same limbs, like a surgeon closing an incision, and a climber who grasps supports on a cliff. This versatility relies on a unique structure that converts neural inputs into muscle force; that is, the motor unit, which comprises an alpha motoneuron and the muscle fibres it innervates (1, 2). The nervous system controls muscle force by modulating the number of active motor units (recruitment) and their firing rates (rate coding) (3). This control strategy involves projecting a substantial level of common synaptic inputs to motoneurons (4-6) that recruits them in a fixed order following the size principle (7).

While we have a clear understanding of the general mechanisms of force control, we still lack a full picture of the detailed organisation of the firing activities of motor units. This is challenged by our inability, so far, to decode the concurrent firing activity of the full spectrum of motor units spanning the entire range of recruitment thresholds in human muscles (8). Indeed, much of our knowledge on the modulation of motor unit activity still derives from animal preparations (9-11) or from serial recordings of motor units in humans, mainly decoded over narrow force ranges (12-14).

Animal studies have shown two stages in the rate coding of individual motor units to increase muscle force (amplification and rate limiting) in response to linear increases in the net excitatory synaptic inputs (11, 15-17). These stages correspond to different responses of motoneurons to ionotropic inputs due to modulation of their intrinsic properties by metabotropic inputs acting on the somato-dendritic surfaces of motoneurons through several ion channels (2, 16, 17). Similar characteristics in firing activity have been observed in humans (12, 13, 18), although these observations have been limited to small samples of motor units (typically < 40) with firing activities decoded over a narrow range of submaximal forces (typically <30% maximal voluntary contraction (MVC) force) (12, 13, 18). It is therefore difficult to infer from these studies a control scheme that is generalisable to entire pool and to all contraction levels.

The combination of arrays of electromyographic (EMG) electrodes with modern source-separation algorithms set the path for the decoding of populations of active motor units - and their motoneurons - in humans (8, 19). Recent studies have identified the concurrent firing activity of up to 60 motor units per contraction and per muscle by increasing the density of electrodes (20, 21). Here, we further extended this approach by tracking motor units across contractions at target forces that ranged from 10 to 80 %MVC, based on the unique spatial distribution of motor unit action potentials across the array of surface electrodes (22, 23). We were able to identify up to ∼200 unique active motor units per muscle and per participant in two human muscles in vivo, yielding extensive samples of motor units that are representative of the entire motoneuron pools (24). With this approach, we described the non-linear transformation of the net excitatory synaptic inputs into firing activity for individual motoneurons (2, 17). We specifically focused on the non-linearities between the changes in firing rate and force: amplification, saturation, and hysteresis (9, 17, 25). The results provided a framework picture of the rate coding of motor units during large increases and decreases in an applied force.

Having access to the rate coding of entire motoneuron pools allowed us to understand how ionotropic and neuromodulatory inputs combine to modulate force, confirming some details of hypothetical force-control schemes in humans. For example, experimental (26, 27) and computational (25) studies have proposed a gain-control mechanism driven by neuromodulatory inputs to motoneurons. One common conclusion of these studies is that pools of motor units involved in fine motor tasks have a low gain, whereas additional motor units activated during more forceful tasks exhibit a progressively increase in gain (17, 25-27). However, this conclusion is based on results from motor units in animals (16) and humans (26, 28, 29) while experimentally modulating inputs, or by fitting a non-linear model predicting muscle activation from excitatory synaptic inputs at multiple contraction levels (27). In contrast, we addressed this question by decoding the concurrent activity of hundreds of motor units spanning the entire range of recruitment thresholds in vivo in humans.

Results

Identification and tracking of individual motor units

16 participants performed either isometric dorsiflexion (n = 8) or knee extension tasks (n = 8) while we recorded the EMG activity of the tibialis anterior (TA - dorsiflexion) or the vastus lateralis (VL – knee extension), with four arrays of 64 surface electrodes (256 electrodes per muscle). The experimental tasks consisted of isometric contractions involving a ramp-up, a plateau, and a ramp-down phase. The ramp-up and ramp-down phases were relatively slow compared with contractile speeds observed during activities of daily living. They were performed at a constant pace of 5 %MVC.s^-1 and the force plateau was maintained for either 10-s (70-80 %MVC), 15-s (50-60 %MVC) or 20-s (10-40 %MVC) (Figure 1A). The level of force during the plateau ranged between 10 and 80 %MVC, with increments of 10 %MVC (randomised order). A source-separation algorithm (8, 19) was applied to the EMG signals to extract motor unit pulse trains, from which we automatically identified the series of discharge times. All identified motor unit firings were visually inspected and manually edited when necessary (30). Because the signal was stationary during the plateau, we were able to estimate reliable separation vectors for large samples of motor units with source-separation algorithms (8, 19) (Figure 1A). The average number of motor units identified in each contraction per participant was 42 ± 24 (25^th–75^th percentile: 24 – 53, up to 95) motor units from the TA and 33 ± 15 (25^th–75^th percentile: 23 – 47, up to 71) motor units from the VL. The datasets from all contraction levels were merged to obtain a sample of unique motor units per muscle and participant spanning the full range of recruitment thresholds (1^st–99^th percentile: 0.9-73.4% MVC; see below).

Figure 1.

To estimate the proportion of the EMG signal represented by the identified motor units, a synthetic EMG signal was reconstructed from the firing activity. To do so, the discharge times were used as triggers to segment the differentiated EMG signals over a window of 25 ms and estimate the averaged action potential waveforms for each motor unit (Figure 1A). The action potentials were then convolved with the discharge times to obtain trains of action potentials, and all the trains of the identified motor units were summed to reconstruct the synthetic EMG signal. Finally, the ratio between the powers of synthetic and original EMG signals was calculated (Figure 1B). This ratio reached 69.3 ± 17.3% (25^th–75^th percentile: 59.3 – 83.6%, up to 94.2%) for TA and 55.2 ± 19.5% (25^th–75^th percentile: 50.0 – 71.9%, up to 86.5%) for VL (Figure 1C), meaning that most of the recorded EMG signals were successfully decomposed into motor unit activity.

To investigate motor unit activity across large force ranges, motor units were tracked between contractions using their unique spatial distribution of action potentials (23) (Figure 1D). This method was validated beforehand using both simulated EMG signals and two-source validation with simultaneous recordings of firing activity from intramuscular and surface EMG signals (Figure S1). On average, we tracked 67.1 ± 10.0% (25^th–75^th percentile: 53.9 – 80.1%) of the motor units between consecutive contraction levels (10% increments, e.g., between 10% and 20% MVC) for TA and 57.2 ± 5.1% (25^th–75^th percentile: 46.6 – 68.3%) of the motor units for VL (Figure S2). There are two explanations for the inability to track all motor units across consecutive contraction levels. First, some motor units are recruited at higher targets only. Second, it is challenging to track small motor units beyond a few contraction levels due to signal cancellation (31, 32). In total, we identified 129 ± 44 motor units (25^th–75^th percentile: 100 – 164, up to 194) per participant in TA and 130 ± 63 motor units (25^th–75^th percentile: 103 – 173, up to 199) per participant in VL.

The accuracy of the tracking method was further tested by confirming the uniqueness of each motor unit within the identified sample, assuming that each motor unit must have a unique representation of their action potentials across the array of surface electrodes (22). This was accomplished by calculating the root mean-square error (RMSE) between their action potentials across the electrodes and the action potentials of the rest of the motor units (Figure 1E). The RMSE between the action potentials of the same motor unit tracked across contractions was calculated as a reference and compared with the 5^th percentile of the distribution of RMSE between motor units (Figure 1F). This reference value typically less than the 5^th percentile of the RMSE calculated with the action potentials of the other motor units (TA: 5.6 ± 2.4 vs. 8.0 ± 1.5%; p<0.001; VL: 6.2 ± 2.8% vs. 8.5 ± 1.6%; p<0.001). Inspection of data for each motor unit revealed that 92.1% (TA) and 87.0% (VL) of the motor units had a lower RMSE between contractions than with the other motor units (<5^th percentile). Of note, motor units with the highest reference values (>95^th percentile) were considered as outliers due to tracking errors and were excluded from the subsequent analyses (Figure 1F). As a result of this approach, a total of 34 motor units were excluded from TA and 28 from VL.

Non-linear rate coding during the linear increase in force – amplification and rate limiting

The input-output function for each motoneuron was characterised as the relation between its instantaneous firing rate and the muscle force (Figure 2A). The linear increase in isometric force was assumed to reflect a proportional increase in the net synaptic excitatory inputs received by the motoneurons, as proposed previously (13, 33). Thus, any deviation from a linear increase in firing rate with respect to a linear increase in force should reflect the influence of neuromodulatory inputs on motoneuron gain (13, 14, 18), or a saturation of the motoneuron firing rate (12, 13). Therefore, the association between the firing rate and the force was characterised by comparing three curve fits: i) a linear fit, which would indicate an increase in firing rate proportional to the increase in the net synaptic excitatory input, ii) a rising exponential fit, which would correspond to an initial acceleration of firing rates followed by full saturation (12, 13, 34), and iii) a natural logarithm fit, which would denote an initial acceleration of firing rate followed by a slower constant increase in firing rate that reflects a rate limiting effect (11, 15, 35) (Figure 2A). Of note, this analysis was performed on each unique motor unit with the data pooled across all contractions (Figure 2A). The result was that 69.5% of the motor units in TA and 72.3% of those in VL had an input-output function better fitted by a natural logarithm (Figure 2B, inset panels).

Figure 2.

The input-output functions of the motor units were then compared based on the distribution of recruitment thresholds (1^st–99^th percentile: 0.9-73.4% MVC). We clustered motor units into three groups with equal ranges of recruitment thresholds: low- (0 to 25% MVC), medium- (25 to 50% MVC), and high- (50 to 75% MVC) threshold motor units (Figure 2B). A significant percentage of input-output functions of low threshold motor units was better fitted with a rising exponential function (12.5% for the TA, 30.9% for the VL), highlighting a steep acceleration of firing rate followed by a full saturation. This pattern was not observed in medium- and high-threshold motor units (0% for TA and VL). In contrast, a significant percentage of input-output functions of high threshold motoneurons was better fitted with a linear function for TA (39.6%), although the input-output functions of high threshold motor units from VL were still better fitted with natural logarithms (87.9%) than linear functions (12.1%).

We further described the first stage of rate coding by estimating the acceleration of firing rate (Figure 2C; Figure S3). The initial acceleration of the firing rate, likely due to the amplification of synaptic inputs by persistent inward currents (11), was calculated as the value of the first derivative of the natural logarithm function at the recruitment threshold (Figure 2C). These values were compared between motor units using a linear mixed effect model, with Muscle (TA; VL) and Threshold (low; medium; high) as fixed effects and Participant as a random effect. There was a significant effect of Muscle (F=30.4; p<0.001), Threshold (F=24.7; p<0.001), and a significant interaction Muscle × Threshold (F=4.1; p=0.017). The initial acceleration of firing rate in TA was greater for the low- (1.0±0.7 pps²; p<0.001) and the high- (1.0±0.3 pps²; p=0.032) threshold motor units than for medium-threshold motor units (0.6±0.3 pps²), with no difference between the low- and high-threshold motor units (p=0.999) (Figure 2C). Similarly, the initial acceleration of firing rate was greater for low-threshold motor units (0.6±0.6 pps²) than for medium- (0.4±0.2 pps²; p<0.001) and the high-threshold motor units (0.4±0.2 pps²; p=0.013) in VL, with no difference between the medium- and high-threshold motor units (p=0.990) (Figure 2C). The initial acceleration of firing rate was also greater for the low- and high-threshold motor units from TA comparing with VL (p<0.001 for both).

Finally, we tested whether the distribution of initial accelerations of firing rate observed across the entire motor unit pool was generalisable to all the participants. Participants with either less than 20 motor units or with no motor units recruited below 5%MVC were excluded from this analysis (TA: n=4; VL: n=3). The function representing the relation between recruitment thresholds and initial accelerations were accurately fitted using the non-linear least-square fitting method (TA; adjusted R-squared: 0.77; VL; adjusted R-squared: 0.80; green lines on Figure 2C). We then fitted the same non-linear functions for each participant; the adjusted R-squared were increased to 0.86 ± 0.04 for the TA and 0.84 ± 0.11 for the VL (inset panels in Figure 2C).

The second stage of rate coding involved a linear increase in firing rate after the initial acceleration phase. This phase was analysed using the average firing rate during the successive force plateaus for each tracked motor unit (Figure 3). Specifically, the increase in firing rate between consecutive contraction levels (i.e., change in force equal to 10% MVC) was compared with a linear mixed effect model with Muscle (TA; VL) and Threshold (low; medium; high) as fixed effects and Participant as a random effect. There were significant effects for Muscle (F=39.7; p<0.001) and Threshold (F=145.3; p<0.001), and a significant interaction Muscle × Threshold (F=37.7; p<0.001). The rate of increase was significantly greater for TA (3.5 ± 1.7 pps; 25^th-75^th percentile: 2.4 – 4.1 pps) than for VL (2.0 ± 1.1 pps; 25^th-75^th percentile: 1.3 – 2.6 pps. p<0.001). The rate of increase was significantly greater in TA for high-threshold than for low- (p<0.001) and medium- (p<0.001) threshold motor units and for medium-than for low-threshold motor units (p<0.001). In contrast, the rate of increase in VL was greater for high- (p<0.001) and medium- (p<0.001) threshold motor units than for low-threshold motor units, with no difference between high- and medium-threshold motor units in VL (p=0.993).

Figure 3.

The data for each participant was characterised with a linear function (Figure 3C). There was a significant positive association between the rate of increase in firing rate and the recruitment threshold for all the participants (Pearson’s r: TA: 0.71±0.14 pps and VL:0.57±0.22 pps; all p<0.022), except in one participant for each muscle (TA: r=0.17; p=0.272; VL: r= 0.12; p=0.411). Most motor units (92.0% in TA and 80.9% in VL) increased firing rates by more than one pulse per second between contractions despite reaching the maximal force tested or until tracking was not possible.

Differences between ascending and descending levels of force - hysteresis

Neuromodulatory inputs can prolong excitatory synaptic inputs and produce a hysteresis between the recruitment and the derecruitment thresholds. The left columns in Figure 4 show the relation between recruitment and derecruitment thresholds within the motor unit pools (Figure 4). The relation for most participants was better characterised with a non-linear than a linear regression (TA: 6 out of 8 participants, adjusted R-squared >0.92); VL: 7 out of 8 participants, adjusted R-squared >0.83), indicating potential differences in the level of hysteresis between motor units of different size. The difference between the recruitement and the derecruitment threshold was compared with a linear mixed effect model with Muscle (TA; VL) and Threshold (low; medium; high) as fixed effects and Participant as a random effect. There was a significant effect of Threshold (F=83.0; p<0.001) and a significant interaction Muscle × Threshold (F=30.1; p<0.001), but no effect of Muscle (F=2.5; p<0.137). We further tested for each group of motor units whether the difference between the two thresholds was significantly different from 0 (absence of hysteresis). Thus, only medium-threshold motor units had a significant positive hysteresis in both muscles (TA: −3.8 ± 7.2 pps; p = 0.013; VL: −3.1 ± 6.0 pps; p = 0.026), whereas low-threshold motor units from both muscles (TA: −1.5 ± 4.8; VL: −1.0 ± 3.2 pps; p>0.221) did not exhibit significant hysteresis. Although high-threshold motor units from TA followed the same trend (−1.7 ± 6.0 pps; p<0.437), high-threshold motor units from VL exhibited a negative hysteresis (VL: +3.2 ± 6.7 pps; p=0.039), that indicated that derecruitment threshold was greater than the recruitment threshold.

Figure 4.

We also compared the relation between motor unit firing rate and force between the ramp-up and ramp-down phases of the isometric contraction. Most motor units (TA: 60.2%; VL: 76.7%) kept the same relation during the two phases. Although the percentage of associations better fitted with linear functions increased for the ramp-down phase for the TA (from 25.4% to 41.9% of motor units), there was no change between ramps for the VL (from 12.8% to 19.4% of motor units). Thus, most motor units kept a non-linear relation during the ramp-down phase.

Discussion

No previous studies have decoded to date the firing activity of >100 motor units spanning the full spectrum of recruitment thresholds, whether in animals or humans. According to previous estimates, the number of motor units in human would be 200±61 for TA and 146±29 for VL (36). We have therefore identified in this study most of the active motor units for each muscle (Figure 1), which has allowed us to infer the rate coding of entire pools of motor units. Overall, we found that motor units within a pool exhibit distinct rate coding with changes in force level (Figure 2 and 3), which contrasts with the long-held belief that rate coding is similar across motor units from the same pool.

The fast initial acceleration phase reflects the amplification of synaptic inputs through the activation of persistent inward currents via voltage-gated sodium and calcium channels (2, 11, 15, 17). Although noradrenergic and serotonergic synapses that generate persistent inward currents are widespread within motor unit pools (37), the influence of the neuromodulatory inputs is greater during this phase for low-than for high-threshold motor units. This is presumably due to a larger ratio between the amplitude of inward currents and the input conductance in low threshold motoneurons (38).

The second phase of rate coding involves a slower linear increase in firing rate, as characterised by the second part of a logarithmic function (Figures 2-3). This phase corresponds to mechanisms previously referred to as ‘rate limiting’ (2, 25, 39) or ‘saturation’ (12, 34), and was more pronounced for low-than high-threshold motor units (Figure 3). This difference may be explained by smaller excitatory synaptic inputs onto low-than high-threshold motoneurons (2, 16), lower synaptic driving potential of the dendritic membrane (12, 40, 41), and longer and larger afterhyperpolarisation phase in low-than high-threshold motoneurons (42-45).

Taken together, these results show how ionotropic and neuromodulatory inputs to motoneurons uniquely combine to generate distinct rate coding across the pool, following a similar pattern across the pool, but with a pace that depends on motor unit size. While the size of the motor unit determines its recruitment threshold, neuromodulatory inputs determine its gain. Thus, in a similar fashion as size imposes a fixed recruitment order (7), neuromodulatory inputs determines a spectrum of variations in motor unit firing rates without the need to differentiate the ionotropic inputs to each motoneuron. (Figure 2C-3-S3) (25, 46, 47). Because high-threshold motor units exhibit a higher gain than low-threshold motor units (Figure 3), their firing rate could eventually reach similar or even greater values than that of low-threshold motor units during strong contractions (47-50). This indicates that the onion skin principle (34, 51) may not hold in all muscles and for all contraction forces (47-50).

It is also worth noting that the profile of rate coding mostly followed the same natural logarithmic trend during the ramp-up and ramp-down phases of the contraction, but with prolonged firing activity only evident for medium-threshold motor units (Figure 4). This hysteresis is possible because of the bistability of the membrane potential mediated by inward currents from calcium channels, the two stable membrane states being the resting potential and a depolarised state that enables self-sustained firing (9, 52). The absence of statistically significant hysteresis in low threshold motoneurons may be simply explained by their early recruitment during the first percentages of force, mathematically lowering the potential amplitude of their hysteresis. Alternatively, prominent outward currents may attenuate the impact of persistent inward currents on the prolongation of their firing activity (53). Similarly, the absence of hysteresis (TA) and even a negative hysteresis (VL) in high-threshold motoneurons can be explained by a briefer membrane bistability due to a faster decay of inward currents from calcium channels, and a narrower range of membrane depolarisation over which these inward currents are activated (52).

The increase in firing rate was also significantly greater for TA motor units than for those in VL. This difference may reflect a varying balance between excitatory/inhibitory synaptic inputs and neuromodulation (14, 25, 39, 46, 47, 54). Specifically, the strength of recurrent and reciprocal inhibitory inputs to motoneurons innervating VL and TA, and their proportional or inverse covariation with excitatory inputs, respectively, may explain the differences in rate limiting and maximal firing rates (14, 25, 39, 46, 47, 54). It is worth noting that similar differences in rate coding were previously observed between proximal and distal muscles of the upper limb (55).

Our results on rate coding characteristics of the entire motor unit pool provide insight on force control strategies. The accurate control of force depends on the firing activity of the population of active motor units, which effectively filters out the synaptic noise of individual motor units and primarily transforms common synaptic inputs into muscle force (4, 5, 56). Because the bandwidth of muscle force is <10 Hz during isometric contractions (57), the activation of inward currents to low-threshold motoneurons at recruitment may enable them to promptly discharge at a rate (>8 pps) that transmits the effective common synaptic inputs (<8-10 Hz) without phase distortion (4, 5, 56). This mechanism facilitates the accuracy of force control.

Moreover, force generation can be described as the filtering of the cumulative firing activity of active motor units with the average twitch force of their muscle units (5, 58). From this perspective, at low forces, recruitment of new motor units has a higher impact on force modulation than rate coding since each recruitment impacts both the cumulative firing activity of the pool and the average twitch force. Rate coding only modulates the cumulative firing activity. Thus, the amplification of the firing rate of low-threshold motor units near their recruitment threshold instantly favours rate coding over the recruitment of additional motor units, which likely allows for smoother force control. Similarly, the progressive recruitment of motor units with a higher gain promotes faster changes in firing rates, which promotes force control accuracy across the full force range. Overall, our results may help to design future non-linear decoders that aim to predict muscle activation or muscle force from descending inputs recorded in the motor cortex, with the will to generalise their performance across movements (27).

Methods

Detailed materials and methods are provided in supplementary files. Data and code are provided at https://github.com/simonavrillon/MUdictionnary.