Common synaptic inputs are not distributed homogeneously among the motor neurons that innervate synergistic muscles

The force generated by the muscles involved in an action is produced by common synaptic inputs received by the engaged motor neurons. The purpose of our study was to identify the low-dimensional latent components, defined hereafter as neural modules, underlying the discharge rates of the motor units from two knee extensors (vastus medialis and lateralis) and two hand muscles (index and thumb muscles) during isometric contractions. The neural modules were extracted by factor analysis from the pooled motor units and no assumptions were made regarding the orthogonality of the modules or the association between the modules and each muscle. Factor analysis identified two independent neural modules that captured most of the covariance in the discharge rates of the motor units in the synergistic muscles. Although the neural modules were strongly correlated with the discharge rates of motor units in each of the synergistic pair of muscles, not all motor units in a muscle were correlated with the neural module for that muscle. The distribution of motor units across the pair of neural modules differed for each muscle: 80% of the motor units in first dorsal interosseous were more strongly correlated with the neural module for that muscle, whereas the proportion was 70%, 60%, and 45% for the thenar, vastus medialis, and vastus lateralis muscles. All other motor units either belonged to both modules or to the module for the other muscle (15% for vastus lateralis). Based on a simulation of 480 integrate-and-fire neurons receiving independent and common inputs, we demonstrate that factor analysis identifies the three neural modules with high levels of accuracy. Our results indicate that the correlated discharge rates of motor units arise from at least two sources of common synaptic input that are not distributed homogeneously among the motor neurons innervating synergistic muscles.

The force generated by the muscles involved in an action is produced by common synaptic inputs 20 received by the engaged motor neurons. The purpose of our study was to identify the low-dimensional 21 latent components, defined hereafter as neural modules, underlying the discharge rates of the motor 22 units from two knee extensors (vastus medialis and lateralis) and two hand muscles (index and thumb 23 muscles) during isometric contractions. The neural modules were extracted by factor analysis from the 24 pooled motor units and no assumptions were made regarding the orthogonality of the modules or the 25 association between the modules and each muscle. Factor analysis identified two independent neural 26 modules that captured most of the covariance in the discharge rates of the motor units in the synergistic 27 muscles. Although the neural modules were strongly correlated with the discharge rates of motor units 28 in each of the synergistic pair of muscles, not all motor units in a muscle were correlated with the neural 29 module for that muscle. The distribution of motor units across the pair of neural modules differed for 30 each muscle: 80% of the motor units in first dorsal interosseous were more strongly correlated with the 31 neural module for that muscle, whereas the proportion was 70%, 60%, and 45% for the thenar, vastus 32 medialis, and vastus lateralis muscles. All other motor units either belonged to both modules or to the 33 module for the other muscle (15% for vastus lateralis). Based on a simulation of 480 integrate-and-fire 34 neurons receiving independent and common inputs, we demonstrate that factor analysis identifies the 35 three neural modules with high levels of accuracy. Our results indicate that the correlated discharge 36 rates of motor units arise from at least two sources of common synaptic input that are not distributed 37 homogeneously among the motor neurons innervating synergistic muscles. 38

Introduction 44
The motor unit is the final common pathway by which an activation signal is transmitted to muscle and 45 transformed into contractile activity (1). As such, all voluntary actions are accomplished by varying the 46 amount of motor unit activity. Despite early claims to the contrary (2, 3), it is not possible to control 47 the activation of individual motor units (4). Instead, synaptic inputs are distributed broadly among the 48 neurons that comprise a motor nucleus and the motor units that are activated in response to these inputs 49 depends on their relative excitability (5-7). As a consequence of this scheme, the order in which motor 50 units are recruited during a voluntary action is relatively fixed (8-10). 51 It is the shared synaptic inputs received by the motor neurons that innervate a muscle and not the activity 52 of individual motor units that is responsible for the force it generates (11,12). In general, the shared 53 inputs can arise from three sources (cortical, brain stem, spinal, and afferent pathways) with varying 54 distributions across different motor nuclei (5, 13-15). One advantage of this scheme is that the shared 55 inputs can engage the motor nuclei of the various muscles involved in an action and thereby facilitate 56 control of the net muscle torque. 57 It has been hypothesised that the control of multiple muscles is achieved by the activation of sets of 58 motor neurons, that have been referred to as "neural modules" or "motor primitives" (16-20). Neural 59 modules emerge from common synaptic inputs, or "neural manifolds" (21), that synergistically activate 60 a group of muscles to perform a specific action. For example, evidence from animal studies indicates 61 that the electrical stimulation of spinal interneurons produces coordinated movements that depend on 62 the location of stimulation (22-24). The modularity of neural control in humans has been estimated by 63 measuring the covariation in muscle activation patterns (EMG signals). The modules extracted by 64 factorization analysis have been termed muscle synergies (19,25) and are assumed to emerge from 65 synaptic inputs that are common to the motor neurons involved in the action. 66 If the synaptic input is shared among the motor neurons that innervate synergistic muscles, it should 67 generate at least one latent neural module based on the covariance in the discharge times of the activated 68 motor units (21). Previous work has addressed this issue by factorizing EMG signals from different 69 muscles (19, 20, 26-29), which assumes that the motor neurons innervating the synergistic muscles 70 receive similar proportions of common synaptic input from one or more sources. 71 The neural modules determining coordinated control of multiple muscles can be investigated by 72 pairwise spike train correlations, an approach that gives access to the full statistical operating principles 73 of a neural network (30-32). The purpose of our study was to identify the low-dimensional latent 74 components, defined hereafter as neural modules, underlying the discharge rates of the motor units from 75 two knee extensors (vastus medialis and lateralis) and two hand muscles (index and thumb muscles) 76 during isometric contractions (19-21, 33). We hypothesized that the discharge rates of the motor 77 neurons innervating each muscle would be explained by more than one neural module. 78 We found that the discharge rates of motor units in individual quadriceps and hand muscles could be 79 characterized by two independent muscle-specific neural modules. The discharge rates of most motor 80 units were associated with the neural module for the muscle in which they resided, but others were 81 correlated with either the neural module for the synergistic muscle or both neural modules. We then 82 simulated the delivery of two independent common synaptic currents into integrate-and-fire motor 83 neurons and to validate our approach to identifying latent components. Our findings provide a greater 84 level of detail about the distribution of common synaptic input within and across the motor nuclei that 85 innervate synergistic muscles. 86 Our approach involved extending the classic method for muscle synergy analysis (17,19,20,25,28,89 34) to motor unit recordings. Instead of treating muscles as individual elements, the discharge times of 90 motor units from different muscles were grouped together. We used a factor analysis that maximizes 91 the correlation between each motor unit and a set of unknown factors, referred to as neural modules. 92 We demonstrate that the factor analysis outperforms other factorization approaches (see Methods), such 93 as principal component analysis and non-negative matrix factorization (19, 35-37), in maximizing the 94 correlation between individual motor unit discharge rates and the latent low-dimensional modules. 95 Theoretical and indirect experimental observations suggest that a common motor command is 96 distributed to sets of motor nuclei (25, 26, 33, 38-41), which results in the discharge rates of motor 97 units across muscles being strongly correlated during voluntary contractions in humans (42-44) and 98 non-human primates (45). In our approach, we did not assume that there is only one latent signal for 99 each muscle (37, 46, 47). Instead, we decoded populations of motor units from surface EMG signals 100 into series of motor unit discharge times during two tasks that involved the synergistic activation of 101 pairs of muscles. 102 The experimental setup and correlation analysis for the two vastii muscles is shown in Figure 1. The 103 motor unit discharge times were decomposed with a blind source separation procedure, which identifies 104 each event with no a-priori knowledge on the physiological information conveyed by the individual 105 motor units (48-50). We identified on average across participants 6.9 ± 4.3 and 4.37 ± 2.34 motor units 106 for VL and VM, respectively during isometric contractions at 10% of maximum. As described in the 107 methods, we subsequently smoothed the motor unit discharge times with a Hann window, which 108 retained all the frequencies responsible for muscle force (<5 Hz (11)). As the applied force has a cut-109 off frequency of ≤20 Hz, the low-pass filtered discharge times (the time series of zeros and ones, Fig.  110 1C) are strongly correlated with the variance in force during steady contractions (11, 51-53). 111 Consequently, we focused on finding the latent components (i.e., the neural modules) for the low-pass 112 filtered signals (Fig. 1D).

119
The motor unit discharge times (series of zeros and ones) were convolved with a 400 ms Hann window, which

126
After converting the discharge times to rates and smoothing the signal, we computed pairwise 127 correlations between each motor unit within the same muscle ( Fig. 1E-F). We consistently found 128 correlated and uncorrelated discharge rates of some motor units from the same vastus muscle, which 129 indirectly indicates that not all motor neurons received the same common input (30, 31). Because most 130 previous studies report high correlations among motor units within a motor nucleus (26, 42, 54), it was 131 necessary to assess the level of cross-talk between muscles. Based on a recently developed method (55, 132 56), we found that the identified motor units had action-potential amplitudes that were statistically 133 similar to the other units in the homonymous muscle and, therefore, were not the result of cross-talk 134 (see Methods).

147
The average number of identified motor units for the hand muscles was 12.2 ± 3.0 and 4.3 ± 1.2 for the 148 first dorsal interosseous and thenar muscles, respectively, across participants. In contrast to the vastii 149 muscles, Figure 2 shows that the discharge rates of all motor units in each hand muscle were strongly 150 correlated (>0.9), and there were few cases of low correlations (see cluster analysis below). Because of 151 differences in the strength of the correlations between motor units in the vastii and hand muscles, we 152 then examined the latent components between motor units across participants and muscles with a 153 factorization approach. 154

Factor analysis reveals a distinct organization of common synaptic inputs 155
Factorization analysis identifies the latent components that covary among sets of variables. This method 156 enables the identification of the potential 'neural modules'. Figure 3 shows the results obtained from 157 the factor analysis for the vastii muscles of two participants. The factor analysis was applied to all motor 158 units from both muscles; therefore, the extracted neural modules did not have any a-priori muscle-159 specific constraint. The latent neural modules are superimposed on each muscle (grey lines indicate 160 individual motor unit discharge rates from that muscle). The first two modules explained most of the 161 signal (>80%); therefore, we only used these two factors in the subsequent analysis. 162 motor unit discharge rates (grey lines, with the mean = 0 spikes/s) and two neural modules derived from 165 the factor analysis (green for the vastus lateralis and violet for vastus medialis). Note the high correlation 166 between the two factors and the discharge rates of some, but not all, motor units for the two subjects. 3. Also note that some motor units are only correlated with one module. 173 We then determined the level of correlation between the discharge rate of each motor unit and the two 174 neural modules (Fig. 3B). This analysis shows, for example, that motor unit number 2 in vastus lateralis 175 for Subject 1 had a stronger correlation with the first neural module, whereas the two other motor units 176 were more correlated with the second factor (Fig. 3D). However, the two modules were not correlated 177 (Fig. 3C). Projections of the two modules ( Fig. 3D) indicated that one motor unit in vastus lateralis was 178 located in the module of the vastus medialis motor units. Subject 3 exhibited more intermingling of the 179 motor unit data in the space of the two neural modules ( Fig. 3

lower right graph). 180
We then looked at the overall distribution of the identified motor units within each neural module across 181 participants for the vastii (n=8) and hand muscles (n=8) (Fig. 4). We found that although many motor 182 units from each muscle shared the same module ( Fig. 4A-D), the discharge rates of some motor units 183 were correlated with both neural modules. There were also some motor units that showed negative 184 correlations with one of the modules. These negative correlations were more common for the hand 185 muscles. for all muscles was the group belonging to the homonymous muscle; that is, a motor nucleus-specific 189 cluster. Interestingly, the proportion of motor units belonging to the shared cluster was greater for vastus 190 lateralis than vastus medialis (Fig. 4B). The motor nucleus-specific cluster was stronger for the hand 191 muscles, with few motor units present in the shared cluster (<20% for both first dorsal interosseous and 192 thenar muscles). Moreover, there were some motor units for both sets of muscles that diverged from 193 the homonymous control and were more correlated with the other synergistic muscle. This was more 194 evident for the vastii (>10% of motor units) than hand muscles (<3% of motor units). 195 length) and sign of the correlation between the discharge rate of one motor unit and the neural module. 199 Note that some motor units shared the same module space (indicated in grey dots), whereas others 200 diverged from synergistic control (blue and red) and a few invaded the territory of the other muscle. B.

201
A cluster analysis identified three main clusters. Note the grey cluster that indicates the percentage of 202 motor units that shared both neural modules. C-D. The same analysis as in A-B but for the hand muscles. 203 Note the smaller proportion of motor units belonging to the shared (grey) cluster in comparison with 204 the vastii muscles. E-F An example of two motor units that occupy different module space. The black 205 ellipse is a visual guidance of the territory occupied by the motor unit 1 (top panel in F), which shows 206 the firing rate of that unit correlated to both module 1 and module 2. In contrast, the lower panel in F 207 shows a motor unit that is only correlated with module 1. 208 We then removed the motor units that shared both neural modules and performed coherence analysis 209 between the motor pools. We found approximately a two-fold decrease in the coherence value without 210 the common motor units. For some subjects, the coherence in any of the physiological bandwidths (0-211 50 Hz), did not differ than for frequencies >50 Hz, which means that there was no coherence between 212 motor units that did not share the same neural module. Conversely, the coherence for the motor units 213 that shared both modules was similar to what previously reported (26, 44). This finding strongly 214 indicates that previous coherence found between muscles from thigh and from the hand is due to the 215 shared inputs that generate the significant coherence value (26, 56). 216

Integrate and fire model: factor analysis accurately reflects the interplay of two common inputs 217
We performed computer simulations to generate a dataset of motor neuron discharge times to determine 218 the optimal convergence and accuracy of factorization analysis on the extraction of neural modules 219 from motor unit data. The aim was to assess the influence of known distributions of two synaptic inputs 220 (Isyn1,2) and their average Isyn3 = (Isyn1 + Isyn2)/2 + independent inputs on the number of identifiable neural 221 modules. The common and independent inputs, as well as the spike times, were approximated by tuning 222 the parameters of an integrate-and-fire model (32, 57). 223 Because we have no information on the dimensions of the latent components, which reflect common 224 inhibitory and excitatory synaptic inputs, we can model these inputs realistically with an integrate-and-225 fire model (32, 57) and study the outputs with pairwise correlations and factor analysis. Moreover, the 226 model allowed us to test the influence of time (total duration of spike times), net synaptic currents, and 227 the strengths of the common and independent inputs. 228 at three time points during the simulation. The absolute difference between the modules corresponds to 249 the accuracy of factor analysis in converging in that specific module. The values for the shared module 250 (green) were close to 0, which indicates perfect separation from the two modules. We injected low 251 percentages of common input (0% indicates that the common and independent input are the same). The 252 dashed vertical lines indicate the range of values that reflect in-vivo motor unit correlation values. There 253 was a strong influence of time, so that 10 s of data were insufficient to obtain reliable estimates of the 254 proportion of common input, whereas there were no differences for the data at 50 s and 80 s. I. Three 255 representative neural modules extracted by factor analysis at three time points (10, 50, and 80 s) when 256 the common input was twice as much as the independent input. 257 We simulated 480 integrate-and-fire neurons that were activated by applying an independent input and 258 a common input. Two-thirds of the population of neurons received the uncorrelated inputs, Isyn1 and 259 Isyn2, and a one-third received Isyn3, which represented the average of the two other inputs (Fig. 5A). We 260 then looked at the correlations between the inputs and outputs (smoothed motor unit discharge rates) as 261 well as the average and standard deviation of the neural modules extracted by the factorization method 262 (Fig. 5). Due to the influence of low-pass filtering of the discharge rates, there was a strong influence 263 of trial duration with the 10-s data being unable to distinguish between the unique and shared inputs.

264
With longer time signals, we were able to retrieve the full dimensions of the three common synaptic 265 input signals by increasing the simulation to 50 s or 80 s (Fig 5. H-I), independently of common and 266 independent input strength. 267

Discussion 268
We analysed the correlation between the discharge times of motor units from different synergistic 269 muscles during isometric contractions with the knee extensors and index finger and thumb muscles. We 270 found two neural modules for the motor units of the vastii muscles, which contrasts with previous 271 findings of only one dominant common input to individual (42) or synergistic muscles (26). As shown 272 in Figure 6, large groups of motor units innervating the VL and VM muscles were associated with 273 specific neural module for each of these muscles, but some motor units were associated with both neural 274 modules. In contrast, fewer motor units innervating the hand muscles (<20%) were associated with both 275 neural modules. Moreover, the discharge rates of some motor units were not correlated with the neural 276 module for the muscle in which they reside, but instead were correlated with the neural module for the 277 synergistic muscle. 278 based on global EMG signals have found strong coherence between the two muscles, indicating a 282 unique common input to the synergistic muscles. After we removed those motor units that shared both 283 neural modules, the correlation between the two pools of motor units was significantly reduced, as 284 indicated in the coherence graph. This indicates that the coherence found in previous studies is mainly 285 attributable to those motor units that shared two distinct sources of common synaptic. B. Visual 286 representation of our current findings. With factorization dimensionality techniques, we found that there 287 are at least two sources of common synaptic input to motor neurons that innervate the vastii and hand 288 muscles. Most motor neurons, but not all of them, innervating each vastus muscle receive common 289 input from a unique source (green and violet lines), but some motor neurons receive inputs from the 290 source directed to the other muscle (dashed green and violet lines; upper graph) and some receive inputs 291 from both sources (lower graph). 292 The correlation between each motor unit and its neural module reveals the potential nature of shared 293 synaptic inputs by the motor neuron pools that is inevitably obscured in the global EMG signal. With 294 this analysis we show that the motor neurons from two hand muscles during an independent task can 295 by fully retrieved by the module they carry, but the motor units for each knee extensor muscle receives 296 common input from two unique sources. 297 Previous experiments reported a single dominant common input governing coordination of the vastus 298 medialis and lateralis muscles (26). Similarly, previous studies on one motor unit pool have identified 299 a single common synaptic input (42, 44). The identification of a single dominant common input in 300 previous studies is based on a pooled-coherence approach to estimate neural connectivity. This analysis 301 averages the correlation between motor unit spike trains with several permutations, therefore, the 302 averaging process inevitably generates significant correlations because ~50% of the homonymous 303 motor unit pool receives a similar input. Because we found that the motor units innervating an individual 304 muscle can receive more than one common input, our results demonstrate that pooled coherence is not 305 an appropriate approach to assess neural connectivity. 306 It is important to note that our experimental task involved isometric contractions in contrast to the 307 dynamic actions typically used to identify 'muscle synergies'. Perhaps, the sources of common input, 308 such as the type and intensity of feedback from sensory receptors (58, 59), differ during isometric and 309 dynamic contractions and the common input received by the motor neurons innervating the synergistic 310 muscles is more homogeneous. For example, we found that in macaque monkeys during rapid high 311 force contractions most motor units share the same neural module (60). 312 Even for isometric contractions, however, the sources of common input may differ with the details of 313 the task being performed. Based on the interpretation that the fluctuations in force during steady 314 isometric contractions are attributable to the variance in the common synaptic input (11, 12, 53), 315 differences in the coefficient of variation for force during a specific action suggest adjustments in the 316 common input across tasks. For example, the coefficient of variation for force during index finger 317 abduction, which is mainly due to the activity of the first dorsal interosseus muscle, was 2x greater 318 when participants performed index finger abduction and wrist extension at the same time even though 319 the abduction force was the same in both tasks (61). Based on the finding of an increase in the coefficient 320 of variation for force during the double-action task (index finger abduction + wrist extension) in the 321 study by Almuklass et al. (2016), it seems reasonable to predict that the neural modules for the two 322 hand muscles in the current study would differ from those observed during the independent actions. 323 Consistent with this possibility, Desmedt and Godaux suggested that the synaptic inputs delivered to 324 the motor neurons that innervate the first dorsal interossei muscle differed when the direction of the 325 force applied by the index finger changed from abduction to flexion (62). The basis for this conclusion 326 was the finding that the recruitment order for some pairs of motor units (~8%) consistently reversed 327 recruitment order when the task was changed from abduction to flexion. They hypothesized that this 328 effect, although relatively modest, was attributable to differences in the distribution of the motor 329 command for each task. 330 Despite the limited scope of the tasks examined in our current study, the findings indicate that the 331 derivation of muscle synergies is based on the common synaptic input that is shared by the motor 332 neurons involved in the action but that this common input is not shared among most of the motor 333 neurons within a given motor nucleus. Moreover, we found that the modulation of discharge rate for all 334 motor units could be classified into three clusters distributed across two neural modules. These results 335 indicate that synergistic motor neuron pools receive common synaptic inputs from at least two different 336 sources during submaximal isometric contractions. 337

Methods 338
Participants 339 Eight subjects were recruited for each experiment (hand and knee extensor). All procedures were 340 approved by the local ethical committees at the University of Rome Foro Italico (approval number 341 44680, knee extension experiments) and Imperial College London (approval number 18IC4685, hand 342 experiments) and conformed to the standards set by the Declaration of Helsinki. The subjects signed an 343 informed consent before participating in the study. Some results from these datasets have been 344 published previously (56, 63). 345 As described subsequently, high-density EMG recordings (Quattrocento, OTBioelettronica, Turin, 346 Italy) and decomposition of the acquired signals (64) were performed in both experiments. 347

Experiment 1 (knee extensors) 348
Participants visited the laboratory on two occasions. In the first visit, they were familiarized with the 349 experimental procedures by performing a series of maximal and submaximal isometric contractions 350 with the knee extensors. In the second visit, which occurred 24 hours after the familiarization session, 351 simultaneous recordings of the force generated by the knee extensors during maximal and submaximal 352 voluntary contractions and HDsEMG signals were recorded from vastus lateralis and vastus medialis. 353 After standardized warm-up contractions, participants were verbally encouraged to push 'as hard as 354 possible' for ∼3-5 s to achieve peak maximal voluntary force (MVC comprised a ramp-up phase (5% MVIF s-1), a plateau (10 s of constant force at target), and a ramp-359 down phase (-5% MVIF s-1). Three minutes of rest was provided between all submaximal contractions. 360 In this study we only used the submaximal steady state contractions at 10% of maximum. 361 All measurements were performed with both legs with the order determined randomly. Participants 362 were asked to avoid exercise and caffeine intake for 48 hours before testing. Participants were 363 comfortably seated and secured in a Kin-Com dynamometer (KinCom, Denver, CO, USA) by means 364 of three Velcro straps (thigh, chest, pelvis), with the knee joint fixed at 45° of flexion (full knee 365 extension at 0°). HDsEMG signals were acquired from the vastii muscles with two grids of 64 electrodes 366 each (5 columns × 13 rows; gold-coated; 1 mm diameter; 8 mm inter-electrode distance; OT 367 Biolettronica, Turin, Italy) (Fig. 1A). Placement of the electrode grids was based on existing guidelines 368 (Barbero et al. 2012) and adjusted as necessary. After shaving and cleaning the skin (70% ethanol), 369 both electrode grids were attached to muscle surfaces using two layers of disposable double-sided foam 370 (SpesMedica, Battapaglia, Italy). Skin-electrode contact was ensured by filling the holes of the foam 371 layer with conductive paste (SpesMedica). A ground electrode was placed on the contralateral wrist, 372 whereas the reference electrodes for both vastus lateralis and vastus medialis grids were attached to the 373 skin over the ipsilateral patella and medial malleolus, respectively. Monopolar HDsEMG signals were 374 recorded using a multichannel amplifier (EMG-Quattrocento, A/D converted to 16 bits; bandwidth 10-375 500 Hz; OT Bioelettronica). 376

Experiment 2 (hand muscles) 377
The experimental setup involved a chair, table, and computer monitor. Participants were comfortably 378 seated with both arms resting on the table. A custom-made apparatus that was secured to the table 379 supported the dominant hand (self-reported) in a position midway between pronation and supination 380 and the forearm and wrist were immobilized. The index finger was aligned with the longitudinal axis 381 of the forearm, and the thumb was held in a resting position at the same height as the index finger. The 382 applied force was displayed on a monitor that was positioned 60 cm in front of the subject. The visual 383 gain was fixed at 66 pixels per percentages of MVC force for each muscle (axis). The forces exerted by 384 the index finger and thumb were measured with a three-axis force transducer (Nano25, ATI Industrial 385 Automation, Apex, NC, USA), digitized at 2048 Hz (USB-6225, National Instruments, Austin, TX, 386 USA), and low-pass filtered with a cut-off frequency of 15 Hz. HDsEMG signals were recorded with a 387 multichannel amplifier (OT Bioelettronica Quattrocento, Turin, Italy; bandwidth: 10-500 Hz; 388 resolution: 16 bits) at a sampling rate of 2048 Hz. Two flexible grids of high-density EMG electrodes 389 (13 × 5 pins, 4 mm interelectrode spacing) were placed on the skin over the FDI and thenar muscles 390 (flexor pollicis brevis and abductor pollicis brevis). 391 Participants performed force-matching tasks (10% MVC force) involving concurrent abduction of the 392 index finger and flexion of the thumb (Fig. 2A). Subjects performed two sustained index finger 393 abduction and thumb flexion contractions for 60 s. Visual feedback was provided as a moving dot cursor 394 in which the x-axis and y-axis corresponded to the thumb and index finger forces, respectively. Subjects 395 had to maintain the force signal within 10% of the target. 396 The experiments began with MVCs (as described in Experiment 1). After the MVCs were determined, 397 the required target forces were displayed on a monitor. Participants performed two 60-s trials with a 30 398 s of rest between trials. As noted in the introduction, we designed our tasks to determine the extent to 399 which distinct motor neuron pools would receive common inputs. To achieve this goal, subjects were 400 instructed to exert forces in the same sagittal plane, which required ~10 minutes of practice. 401

Data analysis 402
The 64 monopolar HDsEMG signals were filtered offline with a zero-lag, high-pass (10 Hz) and low-403 pass filter (500 Hz). The force signals were corrected for the influence of gravity and normalized to 404 MVC force. HDsEMG channels with poor signal-to-noise ratios were inspected with a semi-automated 405 function that identified spurious EMG signals based on the power spectrum. Those channels with a poor 406 signal-to-noise ratio (defined as 3 standard deviations from the mean, power spectrum averaged across 407 all signals in the band 10 -500 Hz) were visually inspected and removed from the analysis. The number 408 of EMG channels containing noise was low; > 95% of the channels had good signal-to-noise ratios. Before the beginning the blind source separation procedure, the spatial sparsity of the matrix x was 430 enhanced by extending the observation numbers. This procedure improves the decomposition as the 431 gradient descent update rule maximises the diversity of the motor unit waveform to converge the 432 discharge times of each motor unit (the sources, s). Because this process is blind, it is possible to inspect 433 the shapes of the motor unit action potentials obtained by spike triggered averaging and to perform 434 visual inspections of the 2D and 3D waveforms (see 50, 60). 435

Factorization analysis 436
The neural control of muscles by motor neurons can be described and predicted analytically. If the 437 discharge times of the motor units are known, it is possible to predict modulation of muscle force with 438 near perfect correlations (52, 65). By recording of large samples of motor units, it is possible to 439 reconstruct modulation of muscle force (11) due the low-pass filtering properties of the muscle to a 440 given neural drive (51, 52). When motor unit discharge rates are filtered in the low-frequency range 441 (muscle bandwidth <20 Hz), it is possible to predict oscillations in the applied force close to ~1% MVC 442 (51). Consequently, the factorization analysis used in the current study focussed on the low-pass filtered 443 motor unit discharge rates. The discharge rates were filtered by convolving with a Hann window of 400 444 ms (2.5 Hz). The motor unit discharge times were factorized with three methods: factorization analysis 445 (66, 67), principal component analysis, and non-negative matrix factorization (see Figure 1 in 446 supplementary materials). 447 These factorization methods were applied on the matrix containing the smoothed discharge rates with 448 rows equal to the number of motor units identified for both muscles and columns equal to the smoothed 449 discharge times. Figures 1 and 2 shows examples of this procedure for the vastii and hand muscles.
where is a non-negative scaling coefficient of the i-th neural module. We were interested in finding 458 the matrix without making any assumptions about the relations between muscles or motor neurons. 459 We found that factor analysis was the best method in terms of correlations of the neural modules to the 460 discharge times of individual motor units. Moreover, we demonstrate with an integrate-and-fire model 461 (see below) that factor analysis can separate the neural modules with high levels of accuracy. We also 462 examined the performance of non-negative matrix factorization and principal component analysis by 463 using previous approaches to identify muscle synergies (i.e., >100 iterations and reconstruction of the 464 original signal, (19, 20, 25)). 465 The factor analysis models the associations between variables into fewer latent variables (factors). It 466 assumes that for a collection of observed variables ( ) there are a set factors (f) that explain most of the 467 total variance, which is the common variance. The function factoran (in Matlab) computes the 468 maximum likelihood estimate of the factor loadings matrix (Λ) 469 where is the vector of independent specific factors. Alternatively, the model can be specified as 471 Where ΛΛ T is the common variance matrix and Ψ = ( ) is the diagonal matrix of specific variances. 473 The unique variance in the model with no a priori assumption of orthogonality between factors (when 474 allowing for factor rotations such as promax) makes the factor analysis an appropriate choice to extract 475 the latent discharge rate of the synergistic motor nucleus. It is supposed that the model mimics the 476 common and independent inputs impinging into the two motor nuclei. 477

Crosstalk and realigning 478
Motor unit action potentials from the first dorsal interosseous into the thenar muscles (and vice versa) 479 and from the vastus medialis into vastus lateralis can experience cross-talk up to 95%; that is, 95% of 480 motor units from one muscle being conducted with minimal shape distortion to the neighbouring muscle 481 (55). Consequently, we examined the level of cross-talk with a validated method (55, 56). Briefly, this 482 procedure takes advantage of the distance from the activated muscle fibers (muscle unit) to the 483 electrode, which is less for motor units in a targeted muscle. Motor units from a targeted muscle are 484 expected to show greater action potential amplitudes with minimal shape distortion (action potential 485 conduction velocity in the grid, see Germer et al. 2021 for the detailed assessment of statistically 486 significant crosstalk of individual motor units). 487 Another step was to realign motor unit action potentials with respect to the averaged motor unit action 488 potential that was obtained after spike-triggered averaging. Because the action potential at individual 489 time instants shows some variability due to surface EMG stochasticity, we convoluted the average 490 action potential to retrieve the time instants of activation of the motor units. Although this procedure 491 was critical for assessing accurate brain-spinal transmission latencies (68), it did not influence our 492 results because corticospinal latencies are so small (<50 ms). In our study, the firing rate was smoothed 493 in the frequency of force production (2.5 Hz, 400 ms). The effects of action potential onset timing, 494 therefore, are removed when low-pass filtering the motor unit discharge timings at <2.5 Hz. 495

Computer simulations 496
We simulated 480 integrate-and-fire motor neurons, each of which received a computer-generated 497 current input (set at 20 nA). The synaptic currents that were shared among all neurons to represent the 498 common synaptic input, but the neurons also received some independent synaptic inputs. Because motor 499 neurons can exhibit synchronous discharge times, the common input currents were close to maximal 500 values in the low frequency range <2.5 Hz (see below). The resting membrane potential and reset 501 voltage was set at -70 mV, the spike threshold was set at -50 mV, and a membrane time constant was 502 20 ms. The timestep duration was set at 0.1 ms. 503 Our model comprised randomly distributed gaussian noises at each time step to represent the common 504 and independent synaptic inputs. Figure 6 shows the overall architecture of the model. Two random 505 uncorrelated gaussian input currents were created at each time step to represent the neural modules that 506 were identified experimentally. One-third of the neurons (160 neurons) received Isyn1 as a unique 507 common input. Isyn2 received the same common input strength but orthogonal to Isyn1. Isyn3 was the 508 average of Isyn1 and Isyn2 plus its unique independent noise (see eq. 5). The input current for each neuron 509 i and population j (j = 1,2,3) can be summarized by the following equation: 510 , = µ + ( ( ) + ( )) (4) 511 where µi is the temporal average of the current that was set at 20 nA and σi sets the global network state 512 by taking into account the unique independent inputs for each cell ( ) and the gaussian-distributed 513 random common inputs ( ) The tuning of these parameters was matched to those observed in vivo.

514
The values of µ, , and were adjusted to reflect physiological values for the variability in motor unit 515 interspike intervals and common input. Interspike interval variability was examined with histogram 516 distributions, as found in the current study and by others (69). The common and independent inputs 517 were tuned based on the cross-correlogram function derived from previously reported motor unit data 518 (44). Therefore, each combination of three randomly assigned groups of 160 neurons from the total 519 pool (n=480) received two independent synaptic currents ( 1,2 equal to equation 4) and a third 520 randomized subpopulation (j = 3) of motor neuron received the average of the two inputs: 521 , = 1 + 2 2 + ( ) eq. 5 522 The bivariate correlations between the synaptic currents are shown in Figure 5. The total duration of 523 the simulation was set at 10, 50, and 80 s. We removed the first and last 2 s of spiking activity for all 524 simulations to minimize the influence of spike-frequency adaptation. The spike trains emitted by the 525 Ith neuron after generating the spike times were stored as a binary time series, which was equal to 1 526 when the neuron reached voltage threshold. The successive analysis followed the same steps as the 527 experimental data. Briefly, the binary spike trains were low-pass filtered with a 2.5 Hann window and 528 the factorization analysis was then performed on the low-pass filtered signals. Because the distribution 529 of inputs to each neuron were known (i.e., 1−3 ), it was possible to retrieve the performance accuracy 530 of the factorization analysis. Moreover, it was possible to investigate the relation between motor neuron 531 responses to increased synaptic currents with changes in common and independent inputs. As reported 532 in the results, we found that when a large number of spikes was included in the analysis (simulation of 533 80 s), the factorization analysis provided a perfect prediction of the three sources of common synaptic 534 inputs. Our model confirmed that the three clusters observed experimentally in the neural modules arise 535 due from two distinct oscillatory common synaptic inputs and that the third component (the shared 536 neural modules) is an average of the two other modules. 537

Statistical analysis 538
We performed linear regression analysis between the smoothed motor unit discharge rates within and 539 between muscles. The significant level was extracted from bivariate Pearson's correlations and 540 Bonferroni method was applied with multiple comparisons. The same procedure was used to find the 541 modules carried by each neuron (decoding function). Each neural module extracted by the factorization 542 algorithm (66, 67) was compared to the firing rate of the individual motor units. Significance was 543 accepted for P values < 0.05. 544