Interfacing Motor Units in Nonhuman Primates Identifies a Principal Neural Component for Force Control Constrained by the Size Principle

Abstract

Motor units convert the last neural code of movement into muscle forces. The classic view of motor unit control is that the CNS sends common synaptic inputs to motoneuron pools and that motoneurons respond in an orderly fashion dictated by the size principle. This view, however, is in contrast with the large number of dimensions observed in motor cortex, which may allow individual and flexible control of motor units. Evidence for flexible control of motor units may be obtained by tracking motor units longitudinally during tasks with some level of behavioral variability. Here we identified and tracked populations of motor units in the brachioradialis muscle of two macaque monkeys during 10 sessions spanning >1 month with a broad range of rate of force development (1.8-38.6 N · m · s^-1). We found a very stable recruitment order and discharge characteristics of the motor units over sessions and contraction trials. The small deviations from orderly recruitment were fully predicted by the motor unit recruitment intervals, so that small shifts in recruitment thresholds happened only during contractions at a high rate of force development. Moreover, we also found that one component explained more than ∼50% of the motor unit discharge rate variance, and that the remaining components represented a time-shifted version of the first. In conclusion, our results show that the recruitment of motoneurons is determined by the interplay of the size principle and common input and that this recruitment scheme is not violated over time or by the speed of the contractions.SIGNIFICANCE STATEMENT With a new noninvasive high-density electromyographic framework, we show the activity of motor unit ensembles in macaques during voluntary contractions. The discharge characteristics of brachioradialis motor units revealed a relatively fixed recruitment order and discharge characteristics across days and rate of force developments. These results were further confirmed through invasive axonal stimulation and recordings of intramuscular electromyographic activity from 16 arm muscles. The study shows for the first time the feasibility of longitudinal noninvasive motor unit interfacing and tracking of the same motor units in nonhuman primates.

Keywords: EMG; common synaptic input; motor units; peripheral stimulation; size principle.

Figure 1.

Motor unit decomposition in awake behaving macaques, experimental framework, and analysis. A, From left to right, 64 monopolar EMG signals during three individual contractions. Each contraction lasted ∼2 s. The monopolar EMG signals were spatially filtered with a double-differential derivation. After this process, blind source separation identified the spike trains belonging to individual motor units. The spike trains for each motor unit were used to spike trigger the average 2D motor unit waveform. The 2D motor unit waveforms were used for the longitudinal tracking, through a 2D cross-correlation function. B, Monkey 1 (MI) individual motoneuron spike trains across the 10 d (color coded). Note that during the different days, we identified a relatively similar number of motor units. The center of the figure shows the experimental setup and an individual voluntary contraction (force signal in red) extracted from day 3. STA, Spike-triggered average. R represents the two-dimensional cross-correlation value.

Figure 2.

Motor unit action potentials and total numbers of identified and tracked motor units across the 10 d (color coded). A, Two-dimensional motor unit action potential propagating under the high-density EMG electrode array. The highlighted yellow inset shows the respective column and row of the high-density EMG matrix during the experiment on the brachioradialis muscle. B, Raster plot of 12 identified motor units (color coded) for seven representative contractions. C, Three-dimensional representation of the motor unit action potential in a specific time instant (highlighted with a red dot in A). Note that each action potential has a unique 3D signature that allows the independent component analysis to converge to the time series of discharge timings of the motor unit. D, Shimmer plots for two action potential waveforms. Each action potential was averaged across an individual contraction and then superimposed across all contractions for a specific day. Note the high similarity across channels for two representative motor units. The left side of the figure shows four EMG channels corresponding to the largest amplitude across the full 64 EMG electrodes, which are shown in d and a. On the right side, E shows a full representative column of the grid (13 monopolar EMG signals). F, H, The total number of identified motor units across the 10 d (black-edged circles) and tracked motor units (open circles with vertical line depicting the SD) for monkeys MI (F) and MA (H). G, I, Bar plot of the number of motor units that were successfully tracked across the 10 d (color coded). Note that the black-edged bar plot corresponds to the number of motor units that were identified at the respective day and used for tracking those motor units in the other days.

Figure 3.

Motor unit discharge characteristics for the tracked motor units. A, The average instantaneous motor unit discharge rate was plotted for all tracked motor units at any given day. Note that some motor units may show different discharge rates because of changes in synaptic input. B, The day-to-day variability was very low (<6%), and this low variability is demonstrated by very high correlation values for the tracked motor units. C, The absolute variability in the discharge rate of the tracked motor units (i.e., the average motor unit discharge rate at day 1 minus the discharge rate of the same motor unit in the other days). Note that this correlation can be significant only if the motor units are tracked successfully, since the motor unit discharge rate shows high variability across the different units (see the figures below). D–F, The same plots as in A–C for monkey (MA). *p < 0.01, **p < 0.001.

Figure 4.

Motor unit recruitment thresholds and intervals across different contractions, motor units, and days. A, For each motor unit, we calculated the shifts in recruitment order with respect to the average recruitment threshold of that unit. The motor units in the two contractions in A are color coded with respect to the recruitment threshold in the first contraction. For example, it is possible to observe the shift in the recruitment order of 2 to 4 in the second contraction. However, these changes only happen for motor units with very similar thresholds. For example, the motor unit (1 in left panel; red), and the highest threshold motor unit (12; green) shows a consistent recruitment order. This can be well appreciated in the following figures when showing the motor unit recruitment order with respect to the average across the specific day. B, Swarm plots of the recruitment order across all motor units. We first computed the recruitment threshold as the first spike of the motor unit during a specific contraction. We then averaged the recruitment threshold across all contractions for the specific motor unit that was tracked across all contractions (each dot in the swarm plot represents the recruitment threshold of a motor unit in an individual contraction). The average recruitment threshold was then used to sort the recruitment interval of all motor units. Note that each motor unit shows a stable behavior across all contractions. C, Three-dimensional swarm plot for all the motor units across the 10 d. For both monkeys, the relationship between recruitment order and motor unit number was linear across the 10 experimental sessions spaced over a month (monkey MI, R = 0.88 ± 0.04; monkey MA, R = 0.88 ± 0.04; p < 0.00001). D, The variability in recruitment order across days and contractions was highly correlated with the recruitment speed of motoneurons. The recruitment speed of motoneurons is an estimate of supraspinal drive and corresponds to the time derivative of the first discharge timings of all motor units during an individual contraction. Each regression line in D shows the variability across contractions for a specific day. Note the high variability in recruitment speed, which indicates the variance in rate of force development across the contractions for a specific day. E, F, Recruitment threshold versus peak-to-peak amplitude of the motor units for monkey MI (E) and monkey MA (F). Each color represents 1 d, and the peak-to-peak amplitude is calculated after spike trigger averaging the high-density EMG signal from the firings of the individual units. The final peak-to-peak value corresponds to the average of the column with the highest activity in the electrode grid. The correlation coefficient is calculated after pooling all days and motor units.

Figure 5.

Encoding of muscle force by motor units. We aimed at decoding and encoding the temporal motor unit information into components by non-negative matrix factorization. A, Raster plot of 12 motor units during a subset of macaque voluntary isometric contractions (gray lines indicate the torque signal). Note the variability in peak forces and the rate of force developments. B–D, The first three contractions in A. C, The motor unit spike trains in B were convoluted with a 2.5 Hz Hanning window. Note the high correlation between the motor unit smoothed discharge rates and muscle force. D, We applied the reduction dimensionality technique non-negative matrix factorization. We constrained the model to learn the components in the motor unit discharge rates up to 10 factors. In this example, the two modules that together explained ∼80% of the variance are shown. Note that these two modules are highly correlated, and time shifted. The inset in D shows the reconstruction accuracy (variance percentage) of the neural modules with respect to the original signal (smoothed motor unit discharge rates). E, We applied cross-correlation analysis between the modules and muscle force. This example shows the correlation between the first module and the second module as well with voluntary force. F, The same method was then applied for the three modules in both monkeys. Note the high correlation across all days and for both monkeys. Moreover, there was always one module with a dominant component (the lag between the different modules was never zero). G, This indicates that there is only one component constrained by the size principle, since the motor unit recruitment thresholds are highly preserved across all contractions. The reconstruction accuracy (variance percentage) is explained across the 10 d for both monkeys.

Figure 6.

Motor unit synchronization across the different contractions for monkeys MI and MA (color coded). A–C, Pipeline for the estimate of motor unit synchronization for an individual contraction. A, Raster plot of 12 motor units (gray indicates the torque with the scale shown in B). B, The discharge timings of the motor units were filtered with a Hanning window of 200 ms. C, The synchronization value was obtained by performing the cross-correlation function between two groups of randomly permutated groups of motor units (number of permutations, 100). Note that the synchronization value was relatively high and comparable to what was observed in humans during rapid force contractions. D, E, Histogram of the synchronization value across the individual contractions for both monkeys. F, The synchronization value was stable across the 10 d (average and SD values for each day are shown). For monkey 2, the first 2 d resulted in a lower synchronization value because of a lower number of identified motor units, as shown previously. Note that the small variability in synchronization value in D and E was fully explained by the instantaneous discharge rate of the motor units, as previously shown (de la Rocha et al., 2007; Del Vecchio et al., 2019b).

Figure 7.

Neuromuscular implants in macaques. A, Both monkeys were implanted bilaterally with a nerve cuff around the median and radial nerves. Implanted intramuscular EMG signals recorded the gross myoelectric activity of 16 muscles bilaterally (8 muscles per side). B, During each experiment, the nerve cuff delivered stimulation pulses at supramaximal intensity (M-waves) and ramped down in small decrements of 0.1 µA. The left side of B (dark green lines) shows the iEMG recording sessions from supramaximal intensity to the smallest intensity (light green). On the right side of the panel, 12 M-waves obtained during the different days (color coded). C, The iEMG signals from the voluntary contractions during one experimental session. Individual contractions as well as the average (black line) are shown. Note the high intertrial variability in gross EMG responses. D, The average iEMG traces across days (color coded), for monkeys MI and MA. E, Non-negative matrix factorization analysis applied to the gross iEMG signals. The neural module that explained most of the variance is shown for each monkey. F, The reconstruction accuracy (variance percentage) of the components extracted by NNMF. Note that one component explained >80% of the variance. G, The cross-correlation of the first two modules for the respective muscles.