A graph-based approach to identify motor neuron synergies

Multiple studies have experimentally observed common fluctuations in the discharge rates of spinal motor neurons, which have been classically interpreted as generated by correlated synaptic inputs. However, so far it has not been possible to identify the number of inputs, nor their relative strength, received by each motor neuron. This information would reveal the distribution of inputs and dimensionality of the neural control of movement at the motor neuron level. Here, we propose a method that generates networks of correlation between motor neuron outputs to estimate the number of common inputs to motor neurons and their relative strengths. The method is based on force-directed graphs, the hierarchical clustering of motor neurons in the graphs, and the estimation of input strengths based on the graph structure. To evaluate the accuracy and robustness of the method, we simulated 100 motor neurons driven by a known number of inputs with fixed weights. The simulation results showed that 99.2 ± 0.6%, 94.3 ± 2.2 %, and 95.1 ± 2.7 % of the motor neurons were accurately assigned to the input source with the highest weight for simulations with 2, 3, and 4 inputs, respectively. Moreover, the normalised weigths (range 0 to 1) with which each input was transmitted to individual motor neurons were estimated with a root-mean-squared error of 0.11, 0.18, and 0.28 for simulations with 2, 3, and 4 inputs, respectively. These results were robust to errors introduced in the discharge times (as they may occur due to errors by decomposition algorithms), with up to 5% of missing spikes or false positives. We finally applied this method on various experimental datasets to demonstrate typical case scenario when studying the neural control of movement. Overall, these results show that the proposed graph-based method accurately describes the distribution of inputs across motor neurons. Authors summary An important characteristics for our understanding of the neural control of natural behaviors if the dimensionality in neural control signals to the musculoskeletal system. This dimensionality in turn depends on the number of synaptic inputs transmitted to the elementary units of this control, i.e., the spinal motor neurons, and on their correlation. We propose a graph-based approach applied to the discharge times of motor neurons to estimate the number of inputs and associated strength transmitted to each motor neuron. For this purpose, we first calculated the correlation between motor neuron smoothed discharge rates, assuming that correlated discharge rates result from the reception of a correlated inputs. Then, we derived networks/graphs in which each node represented a motor neuron and where the nodes were positioned close to each or further apart, depending on the level of correlated activities of the corresponding motor neurons. Using simulations of motor neuron behaviour, we showed that the spatial information embedded in the proposed graphs can be used to accurately estimate the number and the relative strengths of the inputs received by each motor neurons. This method allows to reconstruct the distribution of synaptic inputs to motor neurons from the observed motor neuron activity.


56
Even the simplest of human movements involve the coordinated activity of thousands of motor 57 units to generate the desired muscle forces. To achieve multi-muscle coordination, numerous 58 studies have suggested a robust and low-dimensional control of muscles through the 59 combination of motor modules, or muscle synergies (1, 2, 3, 4). Muscle synergies are functional 60 units that generate a motor output by imposing a specific activation pattern on a group of 61 muscles (5, 6). Generally, their identification relies on the combination of electromyographic 62 (EMG) recordings over numerous muscles with factorization algorithms that extract consistent 63 patterns of activity across muscles (7). It has been hypothesized that each synergy would be 64 associated to a single neural command, which would in turn decrease the computational load 65 of movement control (4,8,9). 66 67 The concurrent developments of high-density EMG electrodes and algorithms that decompose 68 EMG signals into individual motor unit spiking activity have opened new perspectives in the 69 study of movement control (10). It is now possible to isolate the activity of tens of individual 70 motor units from multiple muscles and to estimate their correlated activity (11,12,13,14). 71 Several previous works have shown that motor units from a single muscle (15, 16, 17, 18, 19) 72 or multiple synergistic muscles (11, 12, 13, 15) exhibit common fluctuations of their discharge 73 rates, likely due to the presence of correlated inputs. These observations are in agreement with 74 the muscle synergy theory. However, it has also been observed that a pool of motor neurons 75 innervating the same muscle do not always receive the same inputs (14,20,21,22,23,24). 76 Therefore, we recently proposed that the central nervous system may generate a movement by 77 transmitting a few inputs to groups of motor units (25), rather than to entire pools innervating 78 individual muscles. In this view, the study of human movement control should be done at the 79 motor neuron level, using methods that can identify the structure of the common inputs 80 transmitted to motor neurons, rather than at the global muscle level. 81 82 While the projection of multi-muscle EMG envelopes into low-dimensional spaces can be done 83 using several factorization algorithms (7), the problem of the identification of the structure of 84 inputs to motor neurons has not been extensively addressed. Here we propose a graph-based 85 approach to identify the number and relative strength of the inputs received by each motor 86 neuron (14,26,27). To test the accuracy of this approach, we used a realistic model of a 87 population of 100 motor neurons driven by a known number of inputs with fixed weights for 88 each motor neuron (28). The graph-based approach relied on the identification of significant 89 correlations between individual motor neuron smoothed discharge rates to estimate whether 90 they received common inputs (29). Then, an algorithm generated a force-directed graph by 91 positioning motor neurons with correlated outputs close to each other, while motor neurons 92 exhibiting low correlation in their output discharge times were positioned farther apart (30). 93 We then identified the number of common inputs received by the full population of motor 94 neurons by detecting the unconnected sections of the graph and we estimated the relative 95 strength of the inputs for each motor neuron by measuring the distance that separated the motor 96 neurons from the centroid of each unconnected cluster. We tested the robustness of this graph-97 based approach to altered versions of motor neuron spike trains that mimicked potential 98 inaccuracies of the EMG decomposition approach, such as false identification of action 99 potentials, missing spikes, or a limited number of identified motor neurons. Finally, we present 100 applications of our approach on representative experimental data from various muscles. 101 102

103
A graph-based approach to identify motor neuron synergies

104
The general idea of this approach is to consider motor neurons as nodes connected with springs, 105 such that they attract each other if they exhibit correlated activity (30). We propose here to take 106 advantage of the spatial properties of this network of motor neurons to infer the number of 107 inputs they receive and their relative strength for each motor neuron. 108 109

Overview of the simulation paradigm 110
We simulated 100 S-type motor neurons receiving correlated and independent synaptic inputs 111 using a motor neuron model previously proposed (28). This model has two compartments for 112 the soma and the dendritic tree, respectively, Na + , fast and slow K + conductance that generate 113 action potentials and afterhyperpolarization, with the addition of a L-type Ca ++ channel into 114 the dendritic compartment to consider the non-linearity between motor neuron inputs and 115 outputs due to the activation of persistent inward currents (31). Geometric and electrical 116 parameters of the 100 motor neurons were linearly interpolated from the range of values 117 reported in Table 1 of (28). We defined the total input for each motor neuron as: 118 where µ is a constant offset that represents the mean of the synaptic inputs, N is the number of 120 inputs, ci is the temporal profile of the i th input, wi,j is the relative strength of the i th input for 121 the motor neuron j, and nj is a synaptic input independent from the other motor neurons. We 122 simulated c and n as Gaussian noises with equal variances. The inputs c had a bandwidth of 0-123 2.5Hz while the independent inputs had a bandwidth of 0-50Hz. Mean and variance of the 124 synaptic current were tuned to generate physiologically realistic Poisson distributions of the 125 interspike intervals (18). 126

127
We performed three simulations with N = 2, 3, and 4 inputs and a simulation where motor 128 neurons only received independent inputs (Figure 1). For the simulation with N = 2, 40 motor 129 neurons had either w1,j or w2,j set to 1, i.e., 20 motor neurons received only the first source of 130 common input and 20 other motor neurons received only the second source of common input. 131 For the 60 remaining motor neurons, w1,j was linearly interpolated between 0.99 and 0.01, and 132 w2,j was equaled to 1 -w1,j. For the simulations with 3 or 4 correlated inputs, groups of 10 133 motor neurons received a single source of common inputs. The remaining motor neurons were 134 grouped by 10 with the same set of w values. These values were randomly distributed across 135 sources of correlated inputs to reach a total of 1. The duration of each simulation was set to 60 136 s.

Correlation between motor neurons to identify common inputs 153
The output of each simulation was a matrix of 100 motor neuron spike trains organized row-154 wise and computed as binary vectors, where 'ones' were discharge instances. We computed 155 the smoothed discharge rate of each motor neuron by convolving its spike train with a 400-ms 156 Hanning window. The duration of the Hanning window was chosen to consider the fluctuations 157 of motor neuron discharge rates effective for force modulation (16, 33). The smoothed 158 discharge rates were then high-pass filtered with a cut-off frequency of 0.75 Hz to remove 159 offsets and trends (Figure 1, (16)). 160 161 We assumed that pairs of motor neurons received common inputs when their smoothed 162 discharge rates were significantly correlated (29). To identify these correlated activities, we 163 performed cross-correlation on the smoothed discharge rates with a maximal time lag of 100 164 ms (16). We considered the statistical significance threshold as the 99 th percentile of the cross-165 correlation coefficient distribution generated with resampled versions of all motor neuron spike 166 trains. To this end, we generated random spike trains for each motor neuron by bootstrapping 167 the interspike intervals (random sampling with replacement). This random spike train had the 168 same number of spikes and the same discharge rate (mean and standard deviation) as the 169 original motor neuron spike train. We repeated this step four times per pair of motor neurons. 170 At the end, all the correlation coefficients below the significance threshold were removed from 171 the 100 by 100 matrix of correlation, which was used to generate the force-directed graph 172

219
we only kept the groups of motor neurons that were fully independant to each other. We then calculated the 220 centroid of each of these independant clusters (C), and calculated the distance of each motor neuron to these 221 centroids (D). We considered the number of independant clusters as the number of inputs, and the distance 222 between each motor neuron and the centroids as the relative strength of these inputs for each motor neuron.

224
To estimate the number of clusters in the graph (i.e., inputs), we identified fully independent 225 groups of motor neurons, i.e., groups that did not share any connections or inputs with each 226 other. While the algorithm used to draw the graph pushes these groups to the edges of the two-227 dimensional space, motor neurons that receive several correlated inputs tend to be positioned 228 close to the center of this plane. Therefore, we performed an iterative procedure starting from 229 the closest motor neuron from the centroid of the graph to the farthest. The algorithm deleted 230 one motor neuron per iteration if it was connected to at least one motor neuron belonging to a 231 different cluster. This procedure stopped when it converged to fully independent groups of 232 motor neurons. Then, we estimated the weight for each motor neuron based on the Euclidean 233 distance between the motor neuron and the centroid of each refined cluster. Indeed, as the 234 position of each motor neuron within the two-dimensional space depends on the number of 235 connections it shares with the others, the closer the motor neurons from each centroid, the 236 greater the number of significant correlations they have with the motor neurons from this 237 cluster. We performed this analysis for each cluster to obtain a vector of weights for each input. 238 Finally, the Euclidean distances were normalized between 0 and 1 and subtracted from 1, given 239 that the shorter the distance, the higher the number of significant correlations. 240 241 Robustness and accuracy of the graph-based approach 242 We first verified the assumption that motor neurons with a significant correlation between their 243 smoothed discharge rates received correlated inputs using the simulation with 2 inputs. To this 244 end, we displayed the distributions of differences between the weights of input #1 for pairs 245 with significant and non significant correlations. A difference equal to 1 means that one motor 246 neuron receive 100% of input #1 while the other motor neuron receives 0% of input #1. On the 247 opposite, a difference of 0 means that both motor neurons receive the same level of input #1. 248 We repeated this analysis with correlations calculated over segments of 10s, 20s, 30s, 40s, 50s, 249 and 60s to test the robustness of our approach across different durations of recordings. 250

251
To assess the robustness of the graph-based approach to identify the number of inputs and their 252 relative strength for each motor neuron, we simulated perturbations of motor neuron spike 253 trains that mimicked inaccuracies in the spike identification algorithms (EMG decomposition). 254 First, we generated altered versions of each motor neuron spike trains by randomly removing 255 or adding discharge instances. Specifically, we removed 1% to 9% of the total number of 256 discharge instances or added false discharge instances equal to 1% to 9% of the initial number 257 of discharge instances per motor neuron. We chose this range of values to get slightly lower 258 sensitivity values that those generally observed after automatic EMG decomposition without 259 manual edition (35). We then applied the graph-based approach as described above to identify 260 the number of inputs and their relative weights for each motor neuron. shaved and cleansed with an abrasive pad. Electrode-to-skin contact was maintained with a 296 biadhesive perforated foam layer filled with conductive paste. Signals were recorded in 297 monopolar mode with a sampling frequency of 2,048 Hz, amplified (x150), band-pass filtered 298 (10-500 Hz), and digitised using a 400 channels acquisition system with a 16-bit resolution 299 (EMG-Quattrocento; OT Bioelettronica, Italy). 300 301 Once recorded, the data were decomposed using convolutive blind-source separation (37). In 302 short, a fixed-point algorithm that maximized the sparsity was applied to identify the sources 303 embedded in the electromyographic signals, i.e., the motor unit spike trains. Motor unit spike 304 trains can be considered as sparse sources with most samples being 0 (i.e., absence of spikes) 305 and a few samples being 1 (i.e., spikes). In this algorithm, a contrast function was iteratively 306 applied to the EMG signals to estimate the level of sparsity of the identified source, and the 307 convergence was reached once the level of sparsity did not vary when compared to the previous 308 iteration, with a tolerance fixed at 10 -4 [see (37)

Estimation of inputs to motor neurons 321
Our preliminary analyses aimed at determining whether the significant correlations between 322 smoothed discharge rates accurately reflected the presence of correlated input between motor 323 neurons. To this end, we used the simulations with N = 2 inputs. We report these results for 324 durations ranging from 10s to 60s, as trial durations can vary during experiments and 325 potentially impact the correlation values. We first calculated the significance thresholds 326 defined as the 99 th percentile of the cross-correlation coefficient distribution generated with 327 resampled versions of all motor neuron spike trains. These thresholds were set at 0.37, 0.26, 328 0.21, 0.18, 0.17, and 0.15 for 10s, 20s, 30s, 40s, 50s, and 60s, respectively. Figure 3A depicts 329 the distribution of motor neuron pairs with significant and non-significant correlations as a 330 function of their difference in weight of input #1. On average, pairs of motor neurons which 331 did not exhibit correlated activity had a large difference in weights of 0.87 ± 0.13, 0.89 ± 0.11, 332 0.88 ± 0.12, 0.90 ± 0.10, 0.91 ± 0.09, and 0.92 ± 0.09 for 10s, 20s, 30s, 40s, 50s, and 60s, 333 respectively; 1 being the maximal difference. This means that they received a small proportion 334 of common inputs. Importantly, motor neurons that received only uncorrelated inputs, i.e., with 335 a weight difference equal to 1, never exhibited significant correlation between their discharge 336 rates, as expected theoretically. Moreover, we observed that the trial duration barely impacted 337 the discrimination of pairs of motor neurons receiving common inputs, i.e., 91,9% of the pairs 338 of motor neurons were either always or never correlated across the six durations. Consequently, 339 the shape of force-directed graphs constructed using the matrix of significant correlations was 340 similar across durations ( Figure 3B). We performed all the following analyses using a trial 341 duration of 60s. 342

Estimation of the number of inputs 353
We generated force-directed graphs using the algorithm developed by Fruchterman and 354 Reingold (30). One property of this algorithm is that nodes, i.e., motor neurons, are positioned 355 evenly in a two-dimensional space. In other words, groups of motor neurons loosely connected 356 due to the presence of minimal, or no, correlated inputs tend to be pushed toward the edges of 357 the space while keeping an equal distance between them. Thus, the shape of the graph gives a 358 sense of the number of inputs received by the pool of motor neurons ( Figure 4A). We applied 359 a hierarchical clustering method to identify these groups of motor neurons loosely connected, 360 and to estimate the number of inputs ( Figure 4B). We identified 2, 3, and 4 clusters over the 361 population of 100 motor neurons for simulation of 2, 3, and 4 inputs. In addition, when applied 362 to the simulation where motor neurons received fully independent inputs, the force-directed 363 graph approach did not identify any clear structure, with 11 clusters, and 17 motor neurons not 364 connected to the graph ( Figure 4A and B).

374
After clustering the motor neurons, we used an iterative algorithm that identifies groups of 375 motor neurons fully independent from the other clusters ( Figure 5A; condition "0%"). On 376 average, these refined clusters grouped motor neurons with the relative weight of one the 377 sources of common input reaching 0.98 ± 0.04, 1.0 ± 0.0, and 0.99 ± 0.05 for graphs with 2, 3, 378 and 4 sources of inputs, respectively. Given that a value of 1 means that the motor neuron 379 received only a single source of correlated inputs, each cluster grouped motor neurons sharing 380 almost 100% of the same source of common input. We then assessed the robustness of this 381 method to estimate the number of inputs by perturbating motor neurons spike trains. 382 Specifically, we removed 1% to 9% of the total number of firing instances or added false firings 383 instances equal to 1% to 9% of the initial number of firing instances per motor neuron. The 384 number of inputs was accurately identified for all the perturbed conditions for 2 inputs, for all 385 the perturbed conditions except -8% and +8% for 3 inputs, and from -6% to +9% for 4 inputs 386 ( Figure 5B

Estimation of the relative strengths of the inputs and motor neuron classification. 398
For this analysis we only considered the conditions where the right number of inputs was 399 estimated for all the simulations, i.e., from -5% to +5%. The classification of each motor neuron 400 according to the cluster with the highest weight was barely affected by the perturbation of 401 motor neuron spike trains for the three simulation. The accuracies were 99.2 ± 0.6%, 94.3 402 ± 2.2 %, and 95.1 ± 2.7 %, for 2, 3, and 4 inputs, respectively ( Figure 6C). When considering 403 the condition without any perturbation of the spike trains (0%), the RMSE of the estimated 404 weights reached 0.11, 0.18, and 0.28 for 2, 3, and 4 inputs, respectively ( Figure 6B 0.7, 0.1, 0.1, and 0.1 for inputs 1, 2, 3, and 4, respectively, it must me classified with the input 1.

Impact of the proportion of identified motor neurons. 426
We tested the impact of the proportion of identified motor neurons on the estimation of the 427 number of inputs and their relative weights (Figure 7). To this end, we randomly selected motor 428 neurons out of the initial group of 100 motor neurons by steps of 10 motor neurons and repeated 429 this procedure 50 times per size of graph ( Figure 7A). For the sake of clarity, we only present 430 the results for the simulation with N = 2 common inputs. Our approach always identified two 431 inputs, independently of the number of motor neurons considered in the analysis. RMSE 432 reached 0.11 ± 0.0. Decreasing the number of motor neurons barely impacted the RMSE, with 433 values of 0.10 ± 0.01 and 0.11 ± 0.0 for graphs of 10 and 100 motor neurons, respectively 434 ( Figure 7B). When considering the classification of motor neurons into a cluster, an accuracy 435 of 97.4 ± 2.4% was observed across conditions. Specifically, the accuracy ranged from 96.8 ± 436 5.2 % to 97.4 ± 5.4 % for graphs of 10 and 100 motor neurons, respectively ( Figure 7C

Application to experimental data 448
We finally applied our approach to experimental data previously published and obtained from 449 various lower leg muscles and various isometric motor tasks (14, 36) (Figure 8). For each 450 dataset, only motor neurons with continuous firing activity during the isometric contraction 451 were considered for the analysis. We first computed the cross-correlation coefficient between 452 the smoothed discharge rates for each pair of motor neurons, and thresholded the adjacency 453 matrix with a significance threshold calculated for each dataset (matrices in Figure 8). We then 454 generated the force-directed graphs, applied hierarchical clustering, and used our approach to 455 identify groups of motor neurons fully independent from the other clusters, as a way to estimate 456 the number of common inputs to motor neurons. 457 458 For the ankle plantarflexion task, we identified two clusters of motor neurons innervating the 459 gastrocnemius lateralis (GL) and gastrocnemius medialis (GM) muscles, with one cluster 460 mostly grouping the motor neurons innervating the GL and the other cluster mostly grouping 461 the motor neurons innervating the GM ( Figure 8A). For the ankle dorsiflexion task, the 462 smoothed discharge rates of all the motor neurons innervating the tibialis anterior were highly 463 correlated, leading to the identification of only one cluster grouping all the identified motor 464 neurons ( Figure 8B). For the knee extension task, we observed that motor neurons innervating 465 the vastus medialis and vastus lateralis muscles were intermingled within the graph ( Figure 8C   neurons, and considered that they received correlated inputs if the correlation in their outputs 514 was significant (14, 16, 29, 41). This approach largely attenuates the problem of the non-515 linearity between motor neuron inputs and outputs (29, 33, 42, 43, 44). Indeed, while the 516 strength of correlation between inputs cannot be estimated from the correlation between 517 outputs (42), a significant correlation between inputs very likely corresponds to a certain level 518 of significant correlation between outputs. Thus, we limited the analysis of the output spike 519 trains to the detection of the presence, not the strength, of common inputs. The validity of this 520 assumption was directly tested by proving i) that pairs of motor neurons with a significant 521 correlation always received most of the same inputs (Figure 3), and ii) that no clear graph 522 structures emerged from motor neurons only receiving independent inputs (Figure 4). 523 Interestingly, the duration of the signals used to compute the correlations did not impact the 524 distribution of significant correlations across pairs of motor neurons nor the shape of the graphs 525 ( Figure 3). 526 We also found that our refinement of clusters, with the identification of groups of motor 527 neurons that do not share any connections with the other clusters, can discriminate motor 528 neurons that exclusively receive the same single source of common input ( Figure 5, 'Condition 529 0%). This is particurlary interesting to quantify the sources of common inputs within the 530 population of identified motor neurons, i.e., the dimensionality of the neural control of motor 531 neurons. Indeed, whether this 'common drivers' are consistent across human behaviors, or 532 change depending on the mechanical contrains of the task, would determine their functional 533 relevance to build human behaviors (1, 45). Additionally, a number of sources of common 534 input significantly lower than the number of identified motor neurons would support the theory 535 of synergies (8, 14, 22, 23). Conversely, having a number of 'common drivers' in the same 536 order of dimension as the number of identified motor neurons woud support the hypothesis of 537 a flexible control of single motor neurons (24, 46, 47). 538 Overall, our results made us confident that our approach can accurately map the distribution of 539 correlated inputs across large datasets of motor neurons receiving several sources of common 540 inputs. Moreover, the proposed graph-based method can also be interpreted as a way of 541 visualizing high-dimensional spaces into two-dimensions. Indeed, conventional plots cannot 542 display datasets with more than three dimensions. It comes as an alternative of several methods 543 previously proposed in the field of computational neuroscience (48), biomechanics (49), or 544 machine learning (50) to solve this issue, even if the output plots of the lattest methods can be 545 abstract and difficult to interpret without complementary information. 546 547 The second step consisted in clustering motor neurons densely connected to each other while 548 separating groups of motor neurons loosely connected. Many approaches exist to solve this 549 problem (51), such as Isomap using the method of random walks (52), or the Louvain algorithm 550 that optimizes the modularity (i.e., a metric that quantifies the strength of the clustering by 551 comparing the number of connections within and between clusters (53)). They all have 552 limitations due to stochastic processes or arbitrary choices of parameters (e.g., the resolution 553 of the clusters; (54, 55)). This means that different clusters can emerge after each iteration, or 554 that the optimization process stops at different sizes of clusters. One way to bypass these issues 555 is to cluster the graph at multiple resolutions, going from one partition including all the motor 556 neurons to smaller clusters (32, 56). The addition of consensus clustering, which identifies the 557 partition with the clusters having the highest probability to be found across all the partitions 558 generated by multiple iterations or multiple resolutions, also improves the robustness of the 559 clustering (32, 34). Therefore, we used this method to accurately identify the number of 560 simulated correlated inputs to motor neurons (Figures 4 and 5). Once refined, these clusters 561 only grouped motor neurons that received a single source of the correlated inputs ( Figure 5). 562 The application of the proposed clustering approach to representative experimental data 563 highlighted consistent results with previous studies that used coherence analyses or 564 factorization algorithms. For example, we reported either one common input to both VL and 565 VM muscles (11, 13), or two common inputs to intermingled motor neurons innervating each 566 muscle (23). However, the scruture of the networks can sometime be noisier with a large 567 number of clusters grouping a few motor neurons, as observed for some participants (see 568 Supplementary Figure 1). Complementarity metrics that caracterize clusters as 'strong' or 569 'weak', such as the conductance (57), may help researchers to remove the weakest clusters, 570 and thus select the right level of resolution within the multiple levels detected by hierarchical 571 clustering (32). 572

573
Besides the estimation of the number of inputs, we also proposed a distance-based metric to 574 estimate the relative strength of these inputs for each motor neuron (Figures 6 and 7). 575 Dimensionality reduction on large scale neural recordings is usually achieved using approaches 576 such as principal component analysis (22, 23), non-negative matrix factorization (58), or 577 factorization algorithms (23). These algorithms usually identify neural latents, with time 578 courses and weights, that can explain most of the variance of the recorded motor neuron spiking 579 activity (59). One caveat of such approaches is the necessity of identifying the number of neural 580 latents a priori, using an arbitrary threshold on the level of variance explained by the model (7, 581 59). Alternatively, a lengthy and computationally demanding approach consists of identifying 582 the number of latents that maximizes the cross-validated data likelihood with probabilistic 583 methods (60). Conversely, the proposed approach first estimates the number of inputs and then 584 uses a simple metric based on the Euclidean distance between each node and the centroid of 585 each refined cluster to estimate the relative strength of each input to each motor neuron ( Figures  586   6 and 7). Importantly, our method was robust to the perturbation of the spiking activity of 587 individual motor neurons that mimic the errors of automatic electromyographic decomposition 588 (35, 38). Additonally, the root-mean squared error in the estimate of input gains did not vary 589 when decreasing the number of motor neurons used to build the graph (Figure 7). This is 590 particularly important when considering that currently it is only possible to experimentally 591 identify a subset of the active motor neurons during natural tasks. As such, this method could 592 add a layer of information for reasearchers trying to quantify and describe the dimensionality 593 of the neural control of motor neurons during and across human movements. 594