Interfacing Spinal Motor Units in Non-Human Primates Identifies a Principal Neural Component for Force Control Constrained by the Size Principle

Motor units convert the last neural code of movement into muscle forces. The classic view of motor unit control is that the central nervous system 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 the performance of tasks with some level of behavioural variability. Here we identified and tracked populations of motor units in the brachioradialis muscle of two macaque monkeys during ten sessions spanning over one month during high force isometric contractions with a broad range of rate of force development (1.8 – 38.6 N·m·s-1). During the same sessions we recorded intramuscular EMG signals from 16 arm muscles of both limbs and elicited the full recruitment through neural stimulation of the median and deep radial nerves. 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 observed between motor units with close recruitment thresholds, and only during 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 could be described as a time-shifted version of the first, as it could be predicted from the interplay between the size principle of recruitment and one common input. In conclusion, our results show that motoneurons recruitment is determined by the interplay of the size principle and common input and that this recruitment scheme is not violated over time nor by the speed of the contractions.


INTRODUCTION 57
Theories of motor control are grounded on recording spinal motor unit activity during voluntary force 58 contractions (1-4). Accurate understanding of motor unit function reveals in a direct way the strategies 59 used by the nervous system to control and coordinate muscle forces (4). Generation of force is believed 60 to occur by a combination of recruitment and rate coding of spinal motor neurons. While it is often 61 assumed that recruitment order and rate coding are determined by the size principle (5, 6) and the 62 common inputs that the motor neurons in a pool receive (2), some studies have challenged this view by 63 proposing a more flexible motor unit control (7-9). Although previous evidence supports the size 64 principle during isometric contractions (10, 11), these results have been challenged by the possibility 65 that the motor cortex could provide independent input to spinal motoneurons. Moreover, it is still 66 unclear if the high correlations in motor unit output (2, 12-14) have a functional origin or represent a 67 physiological epiphenomenon. 68 The current lack of definitive evidence for size principle and common input during recruitment with 69 force modulation is due to technical limitations. Accurate measures of the recruitment order and 70 common input necessitate multiple recordings from as many units as possible and the tracking of the 71 same motor units across different days and across rates of muscle force development (4,7,8,10,11,72 15-18). Currently, no studies tracked the same population of motor units in longitudinal experiments 73 during natural tasks in non-human primates. Such tracking of the same population of neurons is crucial 74 to infer functional behaviour. This is even more important when testing intrinsic properties of 75 motoneurons, such as those associated with the size principle. One way to identify motor unit activity 76 during natural tasks is to insert percutaneous wire electrodes into muscles. However, these electrodes 77 may yield limited signal quality and limited number of detected motor units. 78 By tracking the behaviour of the same motor neurons across multiple experimental sessions with a new 79 non-invasive neural interface consisting of high-density grids of electrodes placed on the muscle, we 80 investigated for the first time the variability in motoneuron recruitment and discharge characteristics 81 over a period of one month in two monkeys during natural contractions. The tracking of a relatively 82 large population of spinal motor units during contractions at different rates of force development 83 allowed us to define the neural strategies accomplished by the central nervous system to control muscle 84 force. Moreover, it was possible to investigate the associations between recruitment of motoneurons 85 and estimates of common synaptic inputs. 86 We found a very small day-to-day and trial-to-trial variability in recruitment order and rate coding, 87 suggesting consistent control of the population of motoneuron ensembles. Moreover, with a 88 factorization method we demonstrated that one common input component was sufficient to explain 89 motor unit recruitment. The application of this approach in a primate species with a motor system 90 closely similar to humans opens the future possibility of combining multiple single motor unit 91 measurements with invasive recordings from central pathways. This has the potential to yield 92 substantial new insights into the anatomical source of common drive during different motor tasks. 93

Motor unit decomposition and tracking 95
We describe the strategies of control of macaque motor units and evaluate the performance of a new 96 non-invasive neural interface framework to monitor the changes in the number and properties of 97 longitudinally tracked units over 10 experimental days (gathered over one month) in two animals. 98 We decomposed spike trains of individual motor units from high-density EMG signals using blind 99 source separation techniques ( Figure 1A; see details in Methods). After this process, the spike trains 100 belonging to each decomposed motor unit were used to estimate the average 2D waveform of the 101 corresponding action potentials ( Figure 1A shows one column of the recording grid). The motor unit 102 waveforms were used to track the same motor unit with a 2D cross-correlation function (19,20). Figure  103 1B shows the raster plot of all motor units across the ten days for Monkey MI. The y-axis in Figure 1B  104 shows the total number of identified motoneurons across days (color-coded). The central panel of Figure  105 1B shows an example of force signal and raster plot of the motor units during a contraction. 106

117
On average, each recording session (one per day, ten days in total) lasted 9.8 ± 2.5 min (Monkey MI) 118 and 8.6 ± 2.8 min (Monkey MA). During these sessions, the monkeys performed on average 118.0 ± 119 30.1 (MI) and 103.5 ± 33.9 (MA) contractions, that were used for the subsequent EMG analyses. The 120 monkeys were instructed to reach a target without a specific training on the rate of force development. 121 Therefore, we obtained a relatively large variance in rate of force development and motor unit 122 recruitment speeds across contractions. During these contractions the rate of force development ranged 123 widely, with an average and standard deviation of 6.44 ± 4.00 N•m•s -1 (range 1.86 -38.66 N•m•s -1 ). 124 Moreover, the peak force obtained across days also showed high variability, spanning two-fold 125 maximum EMG amplitudes. 126 We identified a total of 389 motor units (192 MI and 197 for MA) in the individual recordings. Of these, 127 only a subset ( Figure 2) could be tracked and reliably matched with a unit from one or more different 128 days on the basis of a two-dimensional correlation coefficients R>0.7 (see details in Methods). The 129 average number of identified motor units for each experimental session was 19.2 ± 2.97 and 19.7 ± 2.4 130 (mean and standard deviation), for MI and MA respectively. We were able to track on average 9.07 ± 131 1.06 and 8.13 ± 2.08 motor units across all 10 days. Figure 2 shows the total number of identified motor 132 units at each day and the number of tracked motor units across sessions, for the two monkeys. The upper 133 panel of Figure 2 shows examples of 2D and 3D motor unit waveforms as well as the total number of 134 motor units across contractions and days (bottom panels F-I). Figure 2F     The motor unit action potential similarity across sessions was assessed with the two-dimensional (2D) 162 cross-correlation function (see details in Methods). Because the motor unit action potential waveform 163 and motor unit discharge characteristics are independent, we first computed quality measures of 164 decomposition based on the action potential waveform, and successively we computed correlation 165 measures between the tracked motor units firing characteristics (discharge rate and recruitment 166 threshold across days). 167 The consistency of each motor unit action potential that was accepted to belong to the same cluster, was 168 very high (Silhouette measure averaged across all the identified motor units and the 10 days, 0.91 ± 169 0.01 and 0.92 ± 0.01, for MI and MA respectively). Silhouette measures above 0.9 have been associated 170 with highly accurate decomposition with respect to intramuscular EMG signals (21). Moreover, the 171 tracked units across sessions exhibited very high 2D correlation coefficients of the motor unit waveform 172 (>0.7 for the tracked units) and similar discharge rates across the different days. Figure 2D shows the 173 action potentials that were spike-trigger averaged across the individual contractions (all the action 174 potentials for a representative contraction were used to generate the motor unit action potential 175 waveform, across all 64 channels). The variability in the action potential waveforms across contractions 176 for the same day were minimal, with action potentials 2D correlation values always above 0.9. This 177 indicates very high reliability in identifying the same motor unit across contractions. 178 Physiological characteristics of macaque motor units 179 Figure 3 shows the discharge characteristics of the tracked motor units across and between days. The 180 inter-day motor unit discharge rate variability was very low, at 3.51 and 5.41 % for MI and MA 181 respectively. For Monkey MI, the bivariate Pearson correlation coefficients between the average 182 discharge rate across the different days were significant in all cases (P<0.001 after Bonferroni 183 correction, Fig.3A-B). Indeed, the absolute differences in discharge rate across the units over the 184 different days (Fig. 3C) was very low (0.14 ± 3.45 spikes/s). For Monkey MA, the results were similar, 185 although with a smaller number of reliably decomposed motor units during the first two days (Fig. 2B), 186 that resulted in poorer tracking performance during those two days ( implies that for a given synaptic input, motoneurons are recruited according to intrinsic properties (2). 204 However, some current and previous studies suggests a flexible control of spinal motor units in the 205 mammalian nervous system (7, 8), so that a strict recruitment order is seen as a special case of a flexible 206 control. According to this view, it is conceivable that variability in recruitment may occur over multiple 207 experimental sessions where the monkeys are instructed to reach a target force level according to a 208 broad range of contraction speeds. Contrary to this idea, we found a consistent recruitment order of 209 motor units that was maintained across contractions and days ( Figure 5). The recruitment order across 210 the 10 experimental sessions was occasionally violated for motor units with very close recruitment 211 thresholds ( Fig. 5A-D). In these cases, the occasional reversals of recruitment order were highly 212 correlated with the speed of recruitment (and therefore with the rate of force development) (Fig. 5C).

213
With very fast recruitment, the difference in threshold between motor units with close recruitment 214 threshold compresses to very small values so that the variability in synaptic input may likely explain 215 the occasional reversals (that happened in a small range). 216

238
The present results are in accordance with previous human and in-vitro experiments indicating that 239 motor units are recruited in a specific order. We therefore wanted to understand if there are specific 240 patterns in the motor unit discharge timings that control the recruitment and muscle force. We applied 241 a non-negative matrix factorization analysis (22) to the motor unit discharge timings. Because of the 242 large amount of motor unit data, we were able to discern the exact patterns common to all and to sub-243 groups of motor units.

262
The non-negative matrix factorization revealed a principal component that explained ~50% of the 263 variance. This component was present in the activity of virtually all low-thresholds motor units. There 264 was a significant second factor that explained ~25% of the variance. Interestingly, this originated mainly 265 from high threshold motor units and was an undistorted, time-shifted version of the first component. 266 We then performed correlation analysis between all the components (10 in total, see Methods) and 267 looked at the specific weight distributions across the individual motor unit recruitment thresholds. We 268 found that these components were consistently time-shifted and with very high correlation values 269 between each other ( Figure 5F). Moreover, the second component was consistently present only in the 270 high-threshold motor units. These results indicate that motor unit discharge rates during natural tasks 271 in macaque monkeys are driven by one dominant command, which manifests in time-shifted form 272 because of the progressive recruitment imposed by the size principle. Because the motoneuron is a non-273 linear system, the ensemble activity strongly indicates that these common fluctuations must originate 274 from common input from cortical, afferents, or brainstem pathways. We provide strong evidence that a 275 main component drives a pool of macaque brachioradialis motor units that is mediated by the 276 recruitment order of the motor units. 277

Motor unit synchronization 278
It has been reported that the discharge timings of spinal motor units show very high synchronization 279 values (2), which are associated with the generation of muscle force (23). Accordingly, we found high 280 values of motor unit synchronization similar to what is typically observed in humans (24). We analysed 281 synchronization in two frequency bandwidths; one which retains most of the information of the 282 corticospinal pathways, 0-40 Hz (25), and a narrowed one (0-5 Hz), which retains the information that 283 is correlated to force generation (corresponding to the muscle low-pass filtering bandwidth, <5Hz (26)). 284 The cross-correlation value for the low pass filtered signals (5 Hz) at lag 0 was 0.78 ± 0.01 and 0.72 ± 285 0.10 for MI and MA, respectively. The values across the different bandwidths were consistently very 286 high. These values also showed very small deviations across the contractions (1.55% and 2.30% for MI 287 and MA; see Figure 6). Interestingly, the value of synchronization was in the highest portion of the 288 range observed in humans (R = 0.5 -0.8).

303
The high correlation further indicates that the motoneurons likely received a strong common excitatory 304 synaptic input and that this input was stable across days (Fig. 6F).  The potentials evoked by electrical stimulation showed high reliability across days, with negligible 317 deviations around the mean (Figure 7). This demonstrated stability of the recordings over days. On the 318 other hand, the voluntary EMG amplitudes showed very high variability, with some muscles (including 319 the brachioradialis) showing a 2-fold difference in maximal amplitude. This indicated relatively large 320 variability in the way contractions were executed. 321 We then applied the same method for the identification of motor unit components (Fig. 5) to identify 322 the neural modules within the muscles, as classically referred to as muscle synergies (28-30). We found 323 one invariant neural component that explained more than 90% of the variance. This component was 324 present either in the iEMG signals only from the trained limb, or in the combined iEMG signals from 325 both limbs. This result further supports the role of a common input that is distributed between and within 326 motor nuclei that is processed by the size principle and spinal cord circuitries, despite the large 327 variability in the muscle activities.  We found relatively high motor unit discharge rates in macaque monkeys (41.7 ± 1.4, 38.4 ± 2.0 361 spikes/s for MI and MA respectively across the ten days). These discharge rates were higher than 362 those observed in isometric contractions at low and moderate forces in humans (<50 % of maximal 363 voluntary force, <30 spikes/s) (33). Conversely, when related to fast human isometric contractions of 364 the tibialis anterior muscle, the observed rates are similar (40.09 and 42.85 spikes/s, for the non-365 human and human motor units, respectively (31)). 366 The discharge timings of the motor units represent the neural code that generates muscle force. 367 Recordings of motor unit activity during voluntary force contractions allow us to test the recruitment 368 of motor units by the central nervous system in a detailed way, clarifying current debates in motor 369 control. It has been debated for decades whether the common motoneuron fluctuations observed at the 370 motor unit level are an epiphenomenon or have a functional origin. Similarly, the Henneman size 371 principle has been constantly under investigation, due to the lack of in-vivo evidence with contractions 372 at different rates of force development (6-9, 34, 35). These problems arise because of the lack of 373 adequate methods. 374 Here we showed that the neural drive to the muscle is highly structured in a hierarchical fashion. We 375 found strong associations between hierarchy and behaviour, so that for a given common input signal, 376 the motoneurons behave synchronously once they reach their threshold to discharge, likely dictated by 377 the intrinsic motoneuron properties. Our results are in strong accordance with simulations suggesting 378 that the spinal cord decodes inputs from descending pathways by modulating the recruitment and 379 derecruitment of motoneurons (36). The factorization analysis applied to individual motor unit 380 discharge timings and gross intramuscular EMG signals from the trained and untrained limb, revealed 381 that one component explained more than 80% of the variance. The motor unit findings revealed that 382 this component is filtered by size principle. Our results demonstrate the interplay between common 383 synaptic input and size principle. 384 In conclusion we presented a new non-invasive framework to decode populations of single spinal neural 385 cells in macaque monkeys, which allows us to move from simple measures of behaviour (force) to the 386 inputs that determine that behaviour. In addition to being non-invasive, this framework identifies the 387 same motor units across months over the full force range. This is critical since inferring the patterns of 388 motor behaviour by random sampling small population of active units may be inadequate (7, 11, 16, 17,  389 37). We anticipate that this approach may find further utility when combined with invasive recordings 390 of central motor circuits, which can provide direct access to the various putative sources of common 391 drive (38-40). 392

Animals 394
Recordings were performed from two adult female awake behaving monkeys (M. mulatta; monkeys MI 395 and MA, age 6, weight 6.2 and 6.7 kg respectively). All animal procedures were performed under 396 appropriate licences issued by the UK Home Office in accordance with the Animals (Scientific 397 Procedures) Act (1986) and were approved by the Animal Welfare and Ethical Review Board of 398 Newcastle University. 399

Behavioural Task 400
The monkeys were trained to perform an isometric elbow flexion task with their right arm. Monkey 401 MA was also trained to perform this task with her left arm. The forearm was placed into a rigid plastic 402 cast. This was 3D printed from a digital model of the forearm made using a laser scanner (Go!Scan, 403 Creaform 3D, Levis, Quebec, Canada), ensuring a close but comfortable fit. A further support held the 404 upper arm; the supports were attached to the training cage to fix the elbow in 90° flexion, and the 405 forearm in semi-pronation so that the radius and ulnar were oriented in a vertical plane. torque had to be held in this window for 1 s before releasing to obtain a food reward. Auditory cues 415 were used to indicate to the monkey that the exerted force was within the required window, or else it 416 was too high. Auditory feedback was also given to mark the end of the hold period. Recordings were 417 collected from 10 sessions spanning 30 and 24 days for monkey MI and MA, respectively. 418

423
Methylprednisolone was infused to reduce oedema (5.4mg·kg −1 ·hr -1 IV). Blood-oxygen saturation, 424 heart rate, arterial blood pressure (using a non-invasive blood pressure cuff on the leg), core and 425 peripheral temperature and end-tidal CO2 were monitored throughout; ventilation was supported with a 426 positive pressure ventilator. Hartmann's solution was infused to prevent dehydration (total infusion rate 427 including drug solutions 5-10 ml·kg −1 ·h −1 ). Body temperature was maintained at 37°C using a 428 thermostatically controlled heating blanket and also a source of warmed air. Intraoperative prophylactic 429 antibiotics (cefotaxime 20mg·kg -1 IV) and analgesia (carprofen 5 mg·kg -1 SC) were given. 430 In monkey MI, nerve cuff electrodes (Microprobe, Gaithersburg, MD, USA) were implanted around 431 the median and deep radial nerves bilaterally and secured with the integral sutures. Each cuff contained 432 eight contacts, arranged as two sets of four wires placed radially around the inner circumference. A 433 plastic headpiece (TECAPEEK MT CF30, Ensinger, Nufringen, Germany) was manufactured based on 434 an MRI scan to fit the skull and fixed using ceramic bone screws (Thomas Recording Inc, Giessen, 435 Germany) and dental acrylic. Intramuscular electrodes comprising Teflon-insulated stainless-steel 436 wires were implanted in eight arm and forearm muscles bilaterally for gross electromyography (EMG) 437 recording. Specifically, the muscles that were implanted with intramuscular electrodes corresponded 438 to: deltoids, extensor carpi radialis (ECR), extensor digitorum communis (EDC), biceps brachii, 439 brachialis, brachioradialis, flexor carpi radialis, and the flexor digitorum superficialis muscle (FD). The 440 EMG and nerve cuff wires were tunnelled subcutaneously to connectors fixed to the headpiece. Nine 441 weeks after monkey MI's first implant surgery, several wires connected to the deep radial nerve cuffs 442 bilaterally were found to be broken, and stimulation through these cuffs was no longer possible.

443
Replacement cuffs (with three contacts each, organised radially around the inner circumference) were 444 then implanted bilaterally on the radial nerve below the spiral groove in a further brief surgery, again 445 with wires tunnelled subcutaneously to the head. Monkey MA underwent the implant surgery at a later 446 stage to monkey MI and so was implanted with the same three contact cuffs around the median and 447 radial nerves, along with EMG electrodes in the same muscles and fitted headpiece. All recordings were 448 subsequently collected using the three contact nerve cuffs. 449 Post-operative care included a full programme of antibiotic (co-amoxiclav, dose as above) and 450 analgesics (meloxicam, 0.2mg kg -1 oral plus a single dose of buprenorphine 0.02mg kg -1 IM). 451

Nerve cuff stimulation and recording 452
Biopolar current pulses (0.2ms per phase) were delivered through the first and third contacts of the 453 three contact radial cuffs with a bi-phasic constant current isolated stimulator (Model DS4, Digitimer, 454 Hertfordshire, UK). Stimulus current was delivered at supramaximal intensity (0.45mA for monkey MI 455 and 0.4mA for monkey MA) and ramped down in decrements of 0.1µA to threshold intensity. Left and 456 right arms were stimulated in different sessions, following recordings of the motor task. 457

Electrophysiological Recordings 458
Recordings were made from the brachioradialis muscle using a high-density surface EMG grids 459 (GR04MMI305, OT Bioelettronica, Turin, Italy) with 64 electrodes (spacing 4mm). A bi-adhesive foam 460 strip with holes aligned to the matrix was placed on the grid, and the holes filled with conductive paste 461 (CC1, OT Bioelettronica, Turin, Italy). This assembly was then stuck to the skin over the muscle. To 462 ensure good skin contact the forearm was shaved and cleansed with alcohol wipes. The location of the 463 grid on the skin was marked each day with permanent marker pen to ensure reproducible placement 464 from session to session. Standard surface adhesive electrodes (Neuroline 720; Ambu A/S, Ballerup, 465 Denmark) were placed over the flexor and extensor tendons at the wrist to act as reference and ground; 466 in the implanted animal (monkey MI), one of the unused nerve cuff electrodes was used as the ground. 467 The surface grid electrode was connected to a custom printed circuit board containing a 64-channel 468 amplifier (gain 192; bandwidth 30Hz -2 kHz) and an analogue-to-digital convertor (RHD2164; Intan 469 Technologies LLC, Los Angeles, CA, USA). Digitized signals were sent over a serial peripheral 470 interface (SPI) cable to an RHD USB interface board (also Intan Technologies). This allowed data to 471 be captured to a computer hard disc (5 kSamples/s) along with the elbow torque signal and digital 472 markers signalling the phases of task performance and stimulus timing. Voluntary brachioradialis 473 activity was recorded from the grid electrode during performance of the behavioural task (typically 100 474 successful trials per session). Involuntary contractions were recorded by the intramuscular electrodes 475 and the grid electrode during the radial nerve stimulation protocol. 476

Motor unit decomposition and analysis 477
The high-density EMG recordings were offline digitally filtered with a 20-500 Hz Butterworth filter. 478 Semi-automated MATLAB software extracted the area under the power spectrum and amplitude of 479 each of the 64 channels and highlighted the channels with poor signal to noise ratio for visual inspection 480 and exclusion from subsequent analysis. After this procedure, the monopolar signals were used for the 481 decomposition. Identification of the individual motor unit firings was accomplished through a 482 previously proposed algorithm (21), modified for these large datasets to use a graphical processing unit 483 (GPU) running CUDA software (Nvidia Inc, Santa Clara, California, USA). 484 Briefly, this algorithm takes advantage of the unique two-dimensional spatiotemporal features of 485 individual motor unit action potentials, to converge on an estimate of the motor unit spike trains. The 486 decomposition blindly identifies the motor unit firings; only motor units with high silhouette-measure 487 (>0.92 SIL) are initially maintained. SIL represents a qualitative measure of decomposition accuracy 488 which is comparable to the pulse to noise ratio, ranging from 0 to 1, where 1 indicates perfect clustering 489 of the motor unit action potential. The blind source separation procedure leverages the high spatial and 490 temporal dimensionality of motor unit action potentials. This information is used to converge in an 491 iterative way in the unique time-series representation of the firing times of the alpha motoneurons. We 492 briefly describe here the general steps of decomposition. For a more detailed look into the details of 493 high density EMG decomposition, the technical and physiological details have been described 494 previously (41, 42) 495 The EMG signal corresponds to the filtering of the motoneuron action potential by the muscle tissue 496 with some added noise. Therefore, it is possible to represent in a mathematical form the signal that is 497 carried by each channel of a multidimensional arrays of EMG signals. The EMG signal can be described 498 as a convolution of the motoneuron discharge timings (sources) by the muscle tissue (muscle unit action 499 potentials). The sources (s) are the motoneuron axonal action potentials when reaching the muscle fibres 500 and can be written as Dirac delta function. 501 where represent the spike times of the jth motor unit. We can then write the EMG signal in a matrix 503 form (e.g., when recorded with multidimensional arrays such as the high-density EMG grids used in 504 this study) as: 505 where s (k) = [s1 (k), s2 (k),…, sn (k)] T represent the n motor unit discharge times that generate the EMG 507 signal (x) and n is the noise to for each electrode. The matrix H (l) in eq. 2 contains the spatial 508 information of the motor unit action potential and has size m x l with lth sample of motor unit action 509 potentials for the n motor units and m channels (two-dimensional format, hereafter referred to 2D motor 510 unit waveform). The high spatial sampling given by the 64 electrodes further enhanced by extending 511 the observation numbers (41) allows the recovery of the sources in an iterative blind way with a function 512 that maximizes the sparsity between each motor unit action potential (Fig. 1A). This process is obtained 513 in a fully automatic and blind way; therefore, we can inspect the validity of decomposition by spike-514 triggered averaging. With spike-trigger averaging it is also possible to retrieve by correlation analysis 515 the information that is carried by the action potential (H) in different days, in a fully automatic way. By 516 using 2D correlation analysis it indeed possible track motor unit waveform across weeks (32) and even 517 months (24). The motor unit tracking uses the information carried in H to compare across sessions the 518 two dimensional cross-correlations across all possible combinations of motor unit action potentials. The 519 two-dimensional cross-correlation (2D correlation hereafter) is comparable to a one-dimensional cross-520 correlation, but with a weighted average across the time-space features of the motor unit waveforms 521 (see Figure 1). The output of the two-dimensional cross-correlation ranges from 0 to 1, where 1 indicate 522 maximal similarity. For example, two randomly-selected motor units have a two-dimensional cross-523 correlation lower than 0.3 (20). 524

Motor unit characteristics 525
We first displayed all motor units 2D correlation values with R > 0.55 and with a total number of 526 discharge timings (impulses) >100 and visually inspected the waveforms for potential errors. The 527 unique combinations of motor unit waveform that were preserved after this visual inspection stage had 528 R > 0.70. From all the retained motor units, we computed the instantaneous discharge rate (inverse of 529 the inter-spike interval) averaged across the hold period for all trials of the task on a given day. 530 Synchronization of the motor unit pool was also assessed, as the magnitude of the cross-correlation 531 between two equally sized groups of motor unit spike trains. The number of motor units in each group 532 was randomly assigned for a total of 100 permutations. For each iteration two random unique subsets 533 of units were selected for each group (each group being half of the total number of the identified units 534 during a specific contraction). The spike trains (binary signals) for each motor unit group were then 535 summed and smoothed using a Hann window with two corner frequencies of 40 Hz and 2.5 Hz. We 536 chose two Hann window because this cut-off retains most of the oscillatory activity of the motoneuron 537 pool (40 Hz) and the low frequency is mainly associated to the neural drive that is responsible for force 538 production (i.e., the correlation between a force signal and the low pass filtered motor unit discharge 539 timings is minimally distorted by the musculotendinous unit). 540 For the motor unit recruitment threshold estimates, we first looked at the recruitment order (in seconds) 541 of the motor unit during the individual contractions. This was estimated by taking the time point when 542 that unit was active for the first time. We then calculated the average recruitment threshold (in seconds) 543 for all units across all contractions. Afterwards, we labelled each unit from 1 to the maximum number 544 of identified units in a specific contraction (i.e., a motor unit takes the value of 1 if it is the first 545 recruited). Then we plotted the recruitment thresholds for each specific unit across all contractions. 546 Because the labelling is not dependent on the average, if there is a correlation between the average 547 recruitment threshold and the binarized recruitment threshold across all contractions, this relationship 548 indicates the amount of flexibility in recruitment order obtained by the nervous system in a direct way. 549 We then computed the derivative of the recruitment of the motor units. After calculating the recruitment 550 thresholds (in seconds) of all the motor units in each contraction, the recruitment threshold was sorted 551 from the smallest to the largest and we computed the derivative of this vector. The derivative of this 552 vector corresponds to the number of motor units recruited per seconds, which is an estimate of the 553 efferent drive received by the population of motor units, i.e., a faster recruitment speed of motoneurons 554 results in a faster rate of force development (31, 43) . We then associated for each contraction the 555 variability in recruitment order, that was calculated as the standard deviation of the binarized 556 recruitment thresholds versus the motor unit recruitment speed (the first derivative of the recruitment 557 thresholds). If there would be an association between these two variables it would indicate that a faster 558 recruitment (which could be due to higher synaptic input) is associated with a violation in the 559 recruitment order. 560

Factorization of motor unit activities 561
We factorized the motor unit discharge timings with a non-negative matrix factorization method 562 (NNMF, 29). This method can learn specific features in 2D images such human face characteristics or 563 sematic properties of a written text with the use of linear algebra. In the context of neural signals, we 564 constrained this method to learn the unique components in the motor unit discharge rates that are 565 responsible for force production. Figure 6 shows the overall architecture for this analysis. 566 The force level developed by a muscle is driven by the number of motor unit activation signals, which 567 can be represented as time sequences of M dimensional vectors, that correspond to the activation of the 568 motoneurons m(t) in response to common and independent synaptic inputs arising from afferent and 569 efferent volleys. Therefore, we can express the motoneuron behaviour as combinations of N varying 570 synaptic inputs which construct a specific motor unit firing characteristic, or neural module, expressed 571 as { ( )} =1,… 572 where is a non-negative scaling coefficient of the i-th neural module. We are interested in finding the 574 vectors within the low-frequency motor unit discharge rates. Because motor unit firing rates are non-575 negative, we can utilize NNMF (22) to constrain to be non-negative. This procedure maximizes the 576 interpretability of the data since the representation of the neural motor unit ensemble only includes 577 additive and not subtractive combinations, therefore having an output module with the same scale as 578 the input signal. NNMF iteratively finds the non-negative factors W and H with an interactive procedure 579 that minimize the residuals between D (the sources) and W*H. So that W*H is a lower-rank and then low pass filtered at 2.5 Hz (Figure 6). NNMF is an iterative algorithm that starts with random 584 initial value of W and H. Because the root mean square of D can have local minima, we performed up 585 to 1000 iterations to converge to a representative reconstruction of D = W*H. 586 We then evaluated the output of NNMF with different decoding-encoding functions. First, we 587 constrained the number of factors number equal to the number of identified motor units across a specific 588 day. After this initial procedure, we consistently found that >10 factor explained 99% of the variance. 589 The reconstruction accuracy (residual variance or variance explained) was calculated by computing the 590 residuals (D -W*H) and then computing the deviation from the mean (R 2 ). Second, we evaluated the 591 decomposition by looking at the decoding-encoding of the individual neurons with the respect to the 592 matrix W. This analysis was computed by performing the cross-correlation between the low-pass 593 filtered motor unit discharge rates (D) and the individual neural modules (W) extracted by NNMF. The 594 same method was applied on the gross EMG signals from the intramuscular electrode. After 595 rectification and averaging, the average EMG signals for each day were processed by NNMF and the 596 residual variance was calculated in the same way for the motor units ( Figure 7 shows the results and 597 analysis of the intramuscular EMG signals). All of the analyses were performed in MATLAB. 598

Funding 599
Supported by a grant from The William Leech Charity. JI received the support from "la Caixa" 600 Foundation (ID 100010434; fellowship code LCF/BQ/PI21/11830018). 601