Estimation of the firing behaviour of a complete motoneuron pool by combining electromyography signal decomposition and realistic motoneuron modelling

Our understanding of the firing behaviour of motoneuron (MN) pools during human voluntary muscle contractions is currently limited to electrophysiological findings from animal experiments extrapolated to humans, mathematical models of MN pools not validated for human data, and experimental results obtained from decomposition of electromyographical (EMG) signals. These approaches are limited in accuracy or provide information on only small partitions of the MN population. Here, we propose a method based on the combination of high-density EMG (HDEMG) data and realistic modelling for predicting the behaviour of entire pools of motoneurons in humans. The method builds on a physiologically realistic model of a MN pool which predicts, from the experimental spike trains of a smaller number of individual MNs identified from decomposed HDEMG signals, the unknown recruitment and firing activity of the remaining unidentified MNs in the complete MN pool. The MN pool model is described as a cohort of single-compartment leaky fire-and-integrate (LIF) models of MNs scaled by a physiologically realistic distribution of MN electrophysiological properties and driven by a spinal synaptic input, both derived from decomposed HDEMG data. The MN spike trains and effective neural drive to muscle, predicted with this method, have been successfully validated experimentally. A representative application of the method in MN-driven neuromuscular modelling is also presented. The proposed approach provides a validated tool for neuroscientists, experimentalists, and modelers to infer the firing activity of MNs that cannot be observed experimentally, investigate the neuromechanics of human MN pools, support future experimental investigations, and advance neuromuscular modelling for investigating the neural strategies controlling human voluntary contractions.


Introduction
During voluntary muscle contractions, pools of spinal alpha-motoneurons (MNs) convert the synaptic input they receive into a neural command that drives the contractile activity of the innervated muscle fibres, determining limb motion. Identifying the recruitment and firing dynamics of MNs is fundamental for understanding the neural strategies controlling human voluntary motion, with applications in sport sciences [1][2][3], and neurological and musculoskeletal rehabilitation [4][5][6][7][8]. Determining the MN-specific contributions to the MN population activity also allows more realistic control of neuromuscular models [9][10][11][12], investigation of muscle neuromechanics [13,14], prediction of limb motion from MN-specific behaviour [15], or improvement in human-machine interfacing and neuroprosthetics [16,17].
Our understanding of MN pool firing behaviour during human voluntary tasks is however currently limited. While the MN membrane afterhyperpolarization and axonal conduction velocity can be inferred from indirect specialized techniques [18,19], most of the other electrochemical MN membrane properties and mechanisms that define the MN recruitment and discharge behaviour cannot be directly observed in humans in vivo. Analysis of commonly adopted bipolar surface EMG recordings, which often lump the motor unit (MU) trains of action potentials into a single signal assimilated as the neural drive to muscle, cannot advance our understanding of the MN pool activity at the level of single MNs. Our experimental knowledge on the remaining MN membrane properties in mammals is therefore obtained from in vitro and in situ experiments on animals [20]. The scalability of these mechanisms to humans is debated [21] due to a systematic inter-species variance in the MN electrophysiological properties in mammals [22]. Decomposition of high-density EMG (HDEMG) or intramuscular EMG (iEMG) signals [23,24] allows the in vivo decoding in human muscles of the firing activity of individual active motoneurons during voluntary contractions and provide a direct window on the internal dynamics of MN pools. Specifically, the non-invasive EMG approach to MN decoding has recently advanced our physiological understanding of the neurophysiology of human MU pools and of the interplay between the central nervous system and the muscle contractile machinery [25,26] Yet, the activity of all the MNs constituting the complete innervating MN population of a muscle cannot be identified with this technique. High-yield decomposition typically detects at most 30-40 MNs [27], while MU pools typically contain hundreds of MUs in muscles of the hindlimb [20]. The small sample of recorded MNs is besides not representative of the continuous distribution of the MN electrophysiological properties in the complete MN pool with a bias towards the identification of mainly high-threshold MUs. The samples of spike trains obtained from signal decomposition therefore provide a limited description of the role of individual MNs in the firing activity of the complete MN pool.
To allow the investigation and description of specific neurophysiological mechanisms of the complete MN population, some studies have developed mathematical frameworks and computational models of pools of individual MNs. These MN pool models have provided relevant insights for interpreting experimental data [28,29], investigating the MN pool properties and neuromechanics [30][31][32], neuromuscular mechanisms [12,33], and the interplay between muscle machinery and spinal inputs [34]. However, none of these MN pool models have been tested with experimental input data, instead either receiving artificial gaussian noise [32], sinusoidal [31] or ramp [29] inputs, inputs from interneurons [30] or feedback systems [28]. These MN pool models have therefore never been tested in real conditions of voluntary muscle contraction. The forward predictions of MN spike trains or neural drive to muscle obtained in these studies were consequently not or indirectly validated against experimental recordings.
The MN-specific recruitment and firing dynamics of these MN pool models are usually described with comprehensive or phenomenological models of MNs. The biophysical approaches [30][31][32][33], which rely on a population of compartmental Hodgkin-Huxley-type MN models provide a comprehensive description of the microscopic MN-specific membrane mechanisms of the MN pool and can capture complex nonlinear MN dynamics [35]. However, these models are computationally expensive and remain generic, involving numerous electrophysiological channel-related parameters for which adequate values are difficult to obtain in mammalian experiments [36,37] and must be indirectly calibrated or extrapolated from animal measurements in human models [38]. On the other hand, phenomenological models of MNs [10,28,39] provide a simpler description of the MN pool dynamics and rely on a few parameters that can be calibrated or inferred in mammals including humans. They are inspired from the Fuglevand's formalism [29], where the output MN firing frequency is the gaussian-randomized linear response to the synaptic drive with a MN-specific gain. However, these phenomenological models cannot account for the MN-specific nonlinear mechanisms that dominate the MN pool behaviour [20,35,40,41]. MN leaky integrate-and-fire (LIF) models, the parameters of which can be defined by MN membrane electrophysiological properties for which mathematical relationships are available [42], are an acceptable trade-off between Hodgkin-Huxley-type and Fuglevand-type MN models with intermediate computational cost and complexity, and accurate descriptions of the MN macroscopic discharge behaviour [43] without detailing the MN's underlying neurophysiology [44]. While repeatedly used for the modelling and the investigation of individual MN neural dynamics [45][46][47], MN LIF models are however not commonly used for the description of MN pools.
To the authors' knowledge there is no systematic method to record the firing activity of all the MNs in a MN pool, or to estimate from a sample of experimental MN spike trains obtained from signal decomposition the firing behaviour of MNs that are not recorded in the MN pool.
There is no mathematical model of a MN pool that (1) was tested with experimental neural inputs and investigated the neuromechanics of voluntary human muscle contraction, (2) involves a cohort of MN models that relies on MN-specific profiles of inter-related MN electrophysiological properties, (3) is described by a physiologically-realistic distribution of MN properties that is consistent with available experimental data. This limits our understanding of the neural strategies of the whole MN pool and of how the firing activity of each MN contributes to the neural drive to muscle.
In this study, a novel four-step approach is designed to predict, from the neural information of N r MN spike trains obtained from HDEMG signal decomposition, the recruitment and firing dynamics of the N−N r MNs that were not identified experimentally in the investigated pool of N MNs. The model of the MN pool was built upon a cohort of N single-compartment LIF models of MNs. The LIF parameters are derived from the available HDEMG data, are MN-specific, and account for the inter-relations existing between mammalian MN properties [42]. The distribution of N MN input resistances hence obtained defines the recruitment dynamics of the MN pool. The MN pool model is driven by a common synaptic current, which is estimated from the available experimental data as the filtered cumulative summation of the N r identified spike trains. The MN-specific LIF models phenomenologically transform this synaptic current into accurate discharge patterns for the N r MNs experimentally identified and predict the MN firing dynamics of the N−N r unidentified MNs. The blind predictions of the spike trains of the N r identified MNs and the effective neural drive to the muscle, computed from the firing activity of the complete pool of N virtual MNs, are both successfully validated against available experimental data.
Neuroscientists can benefit from this proposed approach for inferring the neural activity of MNs that cannot be observed experimentally and for investigating the neuromechanics of MN populations. Experimenters can use this approach for better understanding their experimental dataset of N r identified discharging MNs. Moreover, this approach can be used by modelers to design and control realistic neuromuscular models, useful for investigating the neural strategies in muscle voluntary contractions. In this study, we provide an example for this application by using the simulated discharge patterns of the complete MN pool as inputs to Hill-type models of muscle units to predict muscle force.

Overall approach
The spike trains sp sim j ðtÞ elicited by the entire pool of N MNs were inferred from a sample of N r experimentally identified spike trains sp exp i ðtÞ with a 4-step approach displayed in Fig 1. The N r experimentally-identified MNs were allocated to the entire MN pool according to their recorded force recruitment thresholds F th i (step 1). The common synaptic input to the MN pool was estimated from the experimental spike trains sp exp i ðtÞ of the N r MNs (step 2) and was linearly scaled for simplicity to the total postsynaptic membrane current I(t) responsible for spike generation. A cohort of N r single-compartment leaky-and-fire (LIF) models of MNs, the electrophysiological parameters of which were mathematically determined by the unique MN size parameter S i , transformed the input I(t) to simulate the experimental spike trains sp exp i ðtÞ after calibration of the S i parameter (step 3). The distribution of the N MN sizes S(j) in the entire MN pool, which was extrapolated by regression from the N r calibrated S i quantities, scaled the electrophysiological parameters of a cohort of N LIF models. The N calibrated LIF models predicted from I(t) the spike trains sp sim j ðtÞ of action potentials elicited by the entire pool of N virtual MNs.
The following assumptions were made. (1) The MU pool is idealized as a collection of N independent MUs that receive a common synaptic input and possibly MU-specific independent noise. (2) In a pool of N MUs, N MNs innervate N muscle units (mUs). (3) In our notation, the pool of N MUs is ranked from j = 1 to j = N, with increasing recruitment threshold. For the N MNs to be recruited in increasing MN and mU size and recruitment thresholds according to Henneman's size principle [20,[48][49][50][51][52][53], the distribution of morphometric, threshold and force properties in the MN pool follows where S is the MN surface area, I th the MN current threshold for recruitment, IR the MU innervation ratio defining the MU size, F th is the MU force recruitment threshold, and f max iso is the MU maximum isometric force. (4) The MN-specific electrophysiological properties are mathematically defined by the MN size S [42]. This extends the Henneman's size principle to: Where C is the MN membrane capacitance, R the MN input resistance and τ the MN membrane time constant.

Experimental data
The four sets of experimental data used in this study, named as reported in the first column of Table 1, provide the time-histories of recorded MN spike trains sp exp i ðtÞ and whole muscle force trace F(t) (left panel in Fig 1), and were obtained from the studies [26,27,54,55], as opensource supplementary material and personal communication, respectively. In these studies, the HDEMG signals were recorded with a sampling rate f s = 2048Hz from the Tibialis Anterior (TA) and Gastrocnemius Medialis (GM) human muscles during trapezoidal isometric contractions. As displayed in Fig 2, the trapezoidal force trajectories are described in this study by Step (1): according to their experimental force thresholds F th i , each MN, ranked from i = 1 to i = N r following increasing recruitment thresholds, was assigned the N th i location in the complete pool of MNs (i!N i mapping).
Step (2): the current input I(t) common to the MN pool was derived from the N r spike trains sp exp i ðtÞ.
Step (3): using I(t) as input, the size parameter S i of a cohort of N r leaky-and-fire (LIF) MN models was calibrated by minimizing the error between predicted and experimental filtered spike trains. From the calibrated S i and the MN i!N i mapping, the distribution of MN sizes S j in the entire pool of virtual MNs was obtained by regression.
Step (4) The train of instantaneous discharge frequency IDF i (t) of the i th identified MN was computed between firing times ft k i and ft kþ1 i as: The IDFs were moving-average filtered by convolution with a Hanning window of length 400ms [56], yielding the continuous filtered instantaneous discharge frequencies (FIDFs) for all N r identified MNs.

Approximation of the TA and GM MU pool size
The typical number N of MUs was estimated for the TA muscle from cadaveric studies [57], statistical methods [58], decomposed-enhanced spike-triggered-averaging (DESTA) approaches [59][60][61][62][63][64], and adapted multiple point stimulation methods [65] in 20-80-year-old human subjects. Because of method-specific limitations [66], results across methods varied substantially, with estimates for N of, respectively, 445,194,190,188 and 300 MUs for the TA muscle. DESTA methods systematically underestimate the innervation ratio due to the limited muscle volume covered by the surface electrodes. Cadaveric approaches rely on samples of small size and arbitrarily distinguish alpha from gamma MNs. Twitch torque measurements are an indirect method for estimating N. Accounting for these limitations, we estimated the potentially conservative N TA = 400 MUs in a typical adult TA muscle. Assuming 200,000 fibres in the TA muscle [67,68], N TA = 400 yields a mean of 500 fibres per TA MU, consistently with previous findings [68]. In two cadaveric studies [57,69], the estimate for the GM was N GM = 550 MUs, which is consistent with N TA = 400 as the GM muscle volume is typically larger than TA's [70]. A sensitivity analysis on the performance of the method presented in Fig 1 upon variations on the value of N TA was included in this study.

PLOS COMPUTATIONAL BIOLOGY
Estimation of the firing behaviour of a complete motoneuron pool

Step (1): MN mapping
In the first step of the approach overviewed in Fig 1, the N r experimentally identified MNs were allocated to the entire pool of N MNs according to their recorded force recruitment thresholds F th i . Three studies measured in the human TA muscle in vivo the force recruitment thresholds F th of MUs, given as a percentage of the maximum voluntary contraction (%MVC) force, for 528 [64], 256 [71], and 302 [72] MUs. Other studies investigated TA MU pools but reported small population sizes [73] and/or did not report the recruitment thresholds [74][75][76].
We digitized the scatter plot in Fig 3 in [64] using the online tool WebPlotDigitizer [77]. The normalized MU population was partitioned into 10%-ranges of the values of F th (in % MVC), as reported in [71,72]. The distributions obtained from these three studies were averaged. The normalized frequency distribution by 10%MVC-ranges of the F th quantities hence obtained was mapped to a pool of N MUs, providing a step function relating each j th MU in the MU population to its 10%-range in F th . This step function was least-squares fitted by a linear-exponential trendline (Eq 3) providing a continuous frequency distribution of TA MU recruitment thresholds in a MU pool that reproduces the available literature data.
Simpler trendlines, such as F th j , returned fits of lower r 2 values. According to the three studies and to [20], a Δ F = 120-fold range in F th was set for the TA muscle, yielding F th (N) = 90%MVC = Δ F �F th (1), with F th (1) = 0.75%MVC. Finally, the equation assigned the N th i locations of the complete pool of N MUs ranked in order of increasing F th : When considering a 100-ms electromechanical delay between MN recruitment time and onset of muscle unit force, the mapping i 2 ½ ½1; N r � � ! N i 2 ½ ½1; N� � did not substantially change. It was therefore simplified that a MN and its innervated muscle unit were recruited at the same time, and the i 2 ½ ½1; N r � � ! N i 2 ½ ½1; N� � mapping derived for muscle units was extrapolated to MNs.
Considering that typically less than 30 MUs (5% of the GM MU pool) can be currently identified by HDEMG decomposition in GM muscles [27], and that the few papers identifying GM MUs with intramuscular electrodes either did not report the MU F th [78][79][80] or identified less than 24 MUs up to 100% MVC [81,82], a GM-specific F th (j) distribution could not be obtained from the literature for the GM muscle. The F th (j) distribution obtained for the TA muscle was therefore used for the simulations performed with the GM muscle, which is acceptable as an initial approximation based on visual comparison to the scattered data provided in these studies.
Step ð2Þ: Current input IðtÞ In the second step of the approach in Fig 1, the common synaptic input to the MN pool was first estimated from the experimental data. The cumulative spike train (CST) was obtained as the temporal binary summation of the N r experimental spike trains sp exp i ðtÞ.
CST ¼ The effective neural drive (eND) to the muscle was estimated by low-pass filtering the CST in the bandwidth [0; 10]Hz relevant for force generation [47]. As the pool of MNs filters the MN-specific independent synaptic noise received by the individual MNs and linearly transmits the common synaptic input (CSI) to the MN pool [31,83], the CSI was equalled to the eND in arbitrary units: The common synaptic control (CSC) signal was obtained by low pass filtering the CSI in [0; 4]Hz.
This approach, which estimates the CSI from the N r experimental spike trains is only valid if the sample of N r MNs is 'large enough' and 'representative enough' of the complete MN pool for the linearity properties of the population of N r MNs to apply [83]. To assess if this approximation holds with the N r MNs obtained experimentally, the following two validations were performed. (1) The coherence coherN r  [47]. (2) The time-histories of the normalized force FðtÞ and common synaptic control CSC, which should superimpose if the linearity properties apply (see Fig 6 in [83]), were compared with calculation of the normalized root-mean-square error (nRMSE) and coefficient of determination r 2 . If coher N r > 0:7, r 2 >0.7 and nRMSE<30%, it was assumed that the sample of N r MNs was large and representative enough of the MN pool for the linearity properties to apply, and the eND was confidently assimilated as the CSI to the MN pool in arbitrary units. It must be noted that if coher N r < 1, the linearity properties do not fully apply for the sample of N r MNs, and the CSI computed from the N r MNs is expected to relate to the true CSI with a coherence close but less than coher N r . The CSI hence obtained reflects the net excitatory synaptic influx common to the MN population, which triggers the opening of voltage-dependent dendritic channels and the activation of intrinsic membrane currents, such as persistent inward currents (PICs), which contribute to the total dendritic membrane current I(t) responsible for spike generation. As the single-compartment LIF model considered in this study does not describe the dendritic activity, the PICrelated nonlinearity in the CSI−I transformation was excluded from the model and I(t) was simplified to be linearly related to CSI with a constant gain G across the MN pool. To avoid confusion, I(t) is referred as 'current input' in the following. The limitations of this simplification are investigated in the Limitations Section (Limitation 3). To identify G, it must be noted that the CSI, which was computed from a subset of the MN pool, does not capture the firing activity of the MNs that are recruited before the smallest identified N th 1 MN, which is recruited at time ft 1 N 1 . It non-physiologically yields CSIðt < ft 1 N 1 Þ ¼ 0. Accounting for this experimental limitation, I(t) was defined to remain null until the first identified MN starts firing at ft 1 N 1 , and to non-continuously reach I th In (Eq 6), the rheobase currents of the first and last identified MNs I th N 1 and I th N r were estimated from a typical distribution of rheobase in a MN pool I th j , obtained as for the distribution of F th , from normalized experimental data from populations of hindlimb alpha-MNs in adult rats and cats in vivo [84][85][86][87][88][89][90][91][92]. It is worth noting that the experimental I th were obtained by injecting current pulses directly into the soma of identified MNs, thus bypassing the dendritic activity and most of the aforementioned PIC-generated nonlinear mechanisms [40]. A Δ I = 9.1-fold range in MN rheobase in [3.9; 35.5]nA was taken [21,42], while a larger Δ I is also consistent with the literature ( [42], Table 4) and larger values of I th can be expected for humans [21].

Step (3): LIF model -parameter tuning and distribution of electrophysiological properties in the MN pool
In the third step of the approach in Fig 1, single-compartment LIF MN models with current synapses are calibrated to mimic the discharge behaviour of the N r experimentally identified MNs.

MN LIF model: description
A variant of the LIF MN model was chosen in this study for its relative mathematical simplicity and low computational cost and its adequacy in mimicking the firing behaviour of MNs. The LIF model described in (Eq 7) describes the discharge behaviour of a MN of rheobase I th and input resistance R as a capacitor charging with a time constant τ and instantaneously discharging at time ft when the membrane voltage V m meets the voltage threshold V th , after which V m is reset and maintained to the membrane resting potential V rest for a time duration called 'inert period' IP in this study. For simplicity without loss of generalisation, the relative voltage threshold was defined as ΔV th = V th −V rest >0, and V rest was set to 0.
The model is described by the following set of equations: The differential equation was solved with a time step dt�0.001s as: This model includes 5 electrophysiological parameters: R, C, τ, ΔV th and IP. The passive parameters R and C were mathematically related to the MN surface area S as R ¼ k R S 2:43 and C = C m �S after an extensive meta-analysis of published experimental data on hindlimb alpha-MNs in adult cats in vivo [42], that sets the constant value of C m to 1.3μF�cm 2 , supports the equality τ = RC and the validity of Ohm's law in MNs as I th ¼ 2:7�10 À 2 R ¼ DV th R , setting the constant value ΔV th = 27mV in this study. The model was thus reduced to the MN size parameter S and to the IP parameter.

MN LIF model: the IP parameter
Single-compartment LIF models with no active conductances receiving current inputs can predict non-physiologically large values of MN firing frequency (FF) because of a linear FF-I gain at large but physiological current inputs I(t) [44]. The saturation in the FF-I relation typically observed in the mammalian motor unit firing patterns is primarily mediated by the voltage-dependent activity of the PIC-generating dendritic channels [40], that was overlooked in the CSI−I transformation in Step (2), and which decreases the driving force of the synaptic current flow as the dendritic membrane depolarizes. While a physiological modelling of this saturation could be achieved with a LIF model with conductance synapses, as discussed in Limitation 3 in the Limitations section, this nonlinear behaviour could also be captured in a simple approach by tuning the phenomenological IP parameter in (Eq 7). Considering a constant supra-threshold current I�I th input to a LIF MN, the steady-state firing frequency FF predicted by the LIF model is: As I and C typically vary over a 10-fold and 2.4-fold range respectively [42], the FF predicted by the LIF model is dominantly determined by the value of IP as the input current increasingly overcomes the MN current threshold: In the previous phenomenological models of MNs [9,10,29], a maximum firing rate FF max was defined and a non-derivable transition from FF(I) to constant FF max was set for increasing values of I(t). Here, IP was integrated to the dynamics of the LIF model and was derived from experimental data to be MN-specific in the following manner. For each of the N r identified MNs, the time-course of the MN instantaneous discharge frequencies (IDFs) was first smoothed, as performed in [55], with a sixth-order polynomial trendline (IDF trend ) to neglect any unexplained random noise in the analysis. The mean M of the trendline values during the plateau of force ½t r 2 ; t r 3 � was then obtained. If 9t 2 ½t r 1 ; t r 2 À 1�; IDF trend ðtÞ > 0:9M, i.e. if the MN reached during the ramp of I(t) in ½t r 1 ; t r 2 À 1� an IDF larger than 90% of the IDF reached one second before the plateau of force, the MN was identified to 'saturate'. Its IP parameter was set to IP ¼ 1 maxðIDF trend Þ , which constrains the MN maximum firing frequency for high input currents to FF max � 1 IP ¼ maxðIDF trend Þ. A power trendline IP(j) = a�j b was finally fitted to the pairs ðN i ; IP N i Þ of saturating MNs and the IP N i values of the non-saturating MNs were predicted from this trendline. To account for the residual variation in FF observed to remain at high I(t) due to random electrophysiological mechanisms, the IP parameter was randomized at each firing time, taking the value IP+o, where o was randomly obtained from a normal gaussian distribution of IP 10 standard deviation.

MN LIF model: MN size parameter calibration
The remaining unknown parameter -the MN size S -defines the recruitment and firing dynamics of the LIF model. The size S i of the i th identified MN was calibrated by minimizing over the time range t tr 0 ; Table 1) the cost function J(S i ) computed as the root-meansquare difference between experimental FIDF exp i ðtÞ and LIF-predicted FIDF sim i ðtÞ filtered instantaneous discharge frequencies: ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffiffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi ffi dt t tr 2 þ t tr 3 2 To assess how well the calibrated LIF models can replicate the available experimental data, the normalized RMS error (nRMSE) (%) and coefficient of determination r 2 between FIDF exp i ðtÞ and FIDF sim i ðtÞ, and the error in seconds between experimental and predicted recruitment times were computed for the N r MNs. Finally, a power trendline in (Eq 11) was fitted to the pairs ðN i ; S N i Þ, and the continuous distribution of MN sizes in the entire pool of N MNs was obtained.
Δ S = 2.4 was taken in [42]. The S(j) distribution defines the continuous distribution of the MN-specific electrophysiological properties across the MN pool ( [42], Table 4).

MN LIF model: parameter identification during the derecruitment phase
The time-range ½t r 3 ; t r 5 � over which the MNs are being derecruited was not considered in the S calibration in (Eq 10) because the MN's current-voltage relation presents a hysteresis triggered by long-lasting PICs, as discussed in [40,93,94]. This hysteresis, which leads MNs to be derecruited at lower current threshold than at recruitment, can be phenomenologically interpreted in the scope of a single-compartment LIF model with no active conductance as an increase in the 'apparent' MN resistance R during derecruitment to a new value R d . To determine R d , a linear trendline I dth ¼ k dth th � I th (k dth th < 1) was fitted to the association of experimental MN recruitment I th and recruitment I dth current thresholds for the N r recruited MNs, and the distribution of MN input resistance R ¼ k R S 2:43 was increased over the derecruitment time range ½t r 3 ; t r 5 � to: As reviewed in [20,40], the current-voltage hysteresis also explains the typically lower MN discharge rate observed at derecruitment than at recruitment. For purely modelling perspectives, this phenomenon was phenomenologically captured by increasing over the time-range ½t r 3 ; t r 5 � the value of the C m property, which is the main parameter influencing the precited MN FF in (Eq 7) for current inputs close to I th . The simulated FIDF sim i traces were iteratively simulated for the N r MNs over the time-range ½t r 3 ; t r 5 � with the N r S i −IP i −calibrated LIF models obtained from Step (3), for 0.1μF�cm 2 incremental changes in the value of the membrane specific capacitance C m . The FIDF sim i results were compared to the FIDF exp i traces with nRMSE and r 2 values. The 'apparent' C m value returning the lowest output for (Eq 13) was retained and was renamed C d m .
In the following, the individual spike trains sp sim i ðtÞ were predicted with C m over the ½t r 0 ; t r 3 ½ time range and C d m over the ½t r 3 ; t r 5 � time range. It must be noted that, although this approach can accurately capture the observed MN nonlinear input-output behaviour and is coherent with the modelling constraints imposed by the LIF formulation in (Eq 7), as the C m and R parameters mainly affect the discharge and recruitment properties of the LIF model respectively in this model, the R and C m parameters are passive MN properties which remain independent from the MN synaptic and membrane activity in actual MNs, while the R-to-R d and C m -to-C d m transformations are merely phenomenological interpretations of complex neurophysiological mechanisms, and do not provide any insights into the working principles underlying MN recruitment. This limitation is further discussed in Limitation 3.

Step (4): Simulating the MN pool firing behaviour
The firing behaviour of the complete MN pool, i.e. the MN-specific spike trains sp j (t) of the N virtual MNs constituting the MN pool, was predicted with a cohort of N LIF models receiving the common synaptic current I(t) as input. The inert period IP j and MN size S j parameters scaling each LIF model were obtained from the distributions IP(j) and S(j) previously derived at Step (3).

Validation
Validation 1. We assessed whether the N r experimental spike trains sp exp i ðtÞ were accurately predicted by this 4-step approach (Fig 1). Steps (2) and (3) (Fig 1) were iteratively repeated with samples of N r −1 spike trains, where one of the experimentally recorded MNs was not considered. At each i th iteration of the validation, the i th identified spike train sp exp i ðtÞ was not used in the derivation of the synaptic current I(t) (step 2) and in the reconstruction of the MN size and IP distributions S(j) and IP(j) in the MN pool (step 3). As in Step (4) Step 4 for the entire MN pool from the N r identified trains sp exp i ðtÞ, accurately predicted the effective neural drive (eND N ) to the muscle. The eND N was computed as the ½0; 4�Hz low-pass filtered cumulative spike train of the N predicted spike trains, CST N ¼ P N j¼1 sp sim j ðtÞ. As suggested for isometric contractions [83], the normalized eND N was compared for validation against the normalized experimental force trace FðtÞ with calculation of the nRMSE and r 2 metrics. FðtÞ was also compared with nRMSE and r 2 to the normalized effective neural drive eND N r computed directly from the N r experimentally identified MN spike trains. The added value of the presented workflow (Steps 1 to 4) in predicting the neural drive to muscle was finally assessed by comparing the (nRMSE, r 2 ) pairs obtained for the eND N r and eND N traces. Considering the uncertainty on the value of N in the literature, as previously discussed, the performance of the method with reconstructed populations of N ¼ fN r ; 100; 200; 300; 400g MNs was investigated. This sensitivity analysis provides different scaling factors N for the normalized MN mapping (Step 1) and IP and S MN parameter distributions (Step 3) and constrains the number of firing MNs involved in the computation of eND N (Step 4).

Application to MN-driven muscle modelling
The N MN-specific spike trains sp sim j ðtÞ predicted in step (4) were input to a phenomenological muscle model to predict the whole muscle force trace F sim N in a forward simulation of muscle voluntary isometric contraction. As displayed in Fig 3, the muscle model was built as N in-parallel Hill-type models which were driven by the simulated spike trains sp sim j ðtÞ and replicated the excitation-contraction coupling dynamics and the contraction dynamics of the N MUs constituting the whole muscle. In brief, the binary spike train sp sim j ðtÞ triggered for each j th MU the trains of nerve and muscle fibre action potentials that drove the transients of calcium ion Ca 2+ and Ca 2+ -troponin complex concentrations in the MU sarcoplasm, yielding in a last step the time-history of the MU active state a j (t). The MU contraction dynamics were reduced to a normalized force-length relationship that scaled nonlinearly with the MU active state [95] and transformed a j (t) into a normalized MU force trace f j ðtÞ. The experiments being performed at constant ankle joint (100˚, 90˚being ankle perpendicular to tibia) and muscle-tendon length, it was simplified, lacking additional experimental insights, that tendon and fascicle length both remained constant during the whole contractile event at optimal MU length l opt;j ¼ 1. The dynamics of the passive muscle tissues and of the tendon and the fascicle force-velocity relationships were therefore neglected. Finally, the MU-specific forces f j (t) were derived with a typical muscle-specific distribution across the MU pool of the MU isometric tetanic forces f max iso ðjÞ [64,71,74]. The whole muscle force was obtained as the linear summation of the MU forces F sim N ðtÞ ¼ P j f j ðtÞ. To validate F sim N ðtÞ, the experimental muscle force F exp (t) was first approximated from the experimental force trace F(t), which was recorded at the foot with a force transducer [26]. The transducer-ankle joint and ankle joint-tibialis anterior moment arms L 1 and L 2 were estimated using OpenSim [96] and a generic lower limb model [97]. Using the model muscle maximum isometric forces, it was then inferred the ratio q of transducer force F(t) that was taken at MVC by the non-TA muscles spanning the ankle joint in MVC conditions. The experimental muscle force was estimated as and was compared with calculation of normalized maximum error ðnMEÞ; nRMSE and r 2 values against the muscle force F sim N ðtÞ predicted by the MN-driven muscle model from the N neural inputs.
The whole muscle force F sim N r ðtÞ was also predicted using the N r experimental spike trains sp exp i ðtÞ as inputs to the same muscle model of N r in-parallel Hill-type models (Fig 3). In this case, each normalized MU force trace f j ðtÞ was scaled with the same f max iso ðjÞ distribution, however assuming the N r MNs to be evenly spread in the MN pool. F sim N r ðtÞ was similarly compared to F exp (t) with calculation of nME, nRMSE and r 2 values. To assess the added value of the step (1-4) approach in the modelling of MN-driven muscle models, the (nME, nRMSE, r 2 ) values obtained for the predicted F sim N r ðtÞ and F sim N ðtÞ were compared.

Experimental data
As reported in Table 1, the experimental datasets DTA 35  The N r MNs were globally derecruited at relatively lower force thresholds ( Fig 4B) and generally discharged at a relatively lower firing rate (Fig 4C) at derecruitment than at recruitment.

Step (1): MN mapping
The N r = 32 MNs identified in the dataset DTA 35 were allocated to the entire pool of N = 400 MNs according to their recruitment thresholds F th i (%MVC). The typical TA-specific frequency distribution of the MN force recruitment thresholds F th , which was obtained from the literature and reported in the bar plot in Fig 5A, was approximated (Fig 5B) [20].
From the F th (j) distribution, the N r identified MNs were mapped to the complete MN pool (blue crosses in Fig 5B) according to their recorded force recruitment thresholds F th i (ordinates in Fig 4B). As shown in Fig 5C, the N r MNs identified experimentally were not homogeneously spread in the entire MN pool ranked in the order of increasing force recruitment thresholds, as two MNs fell in the first quarter of the MN pool, 5 in the second quarter, 18 in the third quarter and 5 in the fourth quarter. Such observation was similarly made in the three other experimental datasets, where no MN was identified in the first quarter and in the first half of the MN pool in the datasets HGM 30 and HTA 50 respectively (second column of Table 2). In all four datasets, mostly high-thresholds MNs were identified experimentally.

Step (2): Common current input I(t)
To approximate the common synaptic input CSI(t) to the MN pool, the CST and the eND to the MN pool were obtained in Fig 6 from the N r MN spike trains identified experimentally using (Eq 4) and (Eq 5). After 20 random permutations of N r 2 complementary populations of MNs, an average coherence of coherN r 2 ;mean ¼ 0:56 was obtained between N r 2 -sized CSTs of the DTA 35 dataset. From Fig 2 in [47], a coherence of coher N r ¼ 0:8 is therefore expected between the CST in Fig 6A and a typical CST obtained with another virtual group of N r = 32 MNs, and by extension with the true CST obtained with the complete MN pool. The normalized eND and force traces (black and green curves respectively in Fig 6B) compared with r 2 = 0.92 and nRMSE = 20.0% for the DTA 35 dataset. With this approach, we obtained coher N r > 0:7, r 2 >0.7 and nRMSE<30% for all four datasets, with the exception of HGM 30 for which coher N r < 0:7 (third column of Table 2). For the TA datasets, the sample of N r identified MNs was therefore concluded to be sufficiently representative of the complete MN pool for its linearity property to apply, and the eND (red curve in Fig 6B for the dataset DTA 35 ) in the bandwidth [0,10]Hz was confidently identified to be the common synaptic input (CSI) to the MN pool. As observed in Fig 6B, a non-negligible discrepancy, which partially explains why coher N r 6 ¼ 1, was however systematically obtained between the eND and force traces in the regions low forces, where mostly small low-threshold MNs are recruited. As discussed in the Discussion Section, this discrepancy reflects the undersampling of small MUs identified from decomposed HDEMG signals and a bias towards mainly identifying the large high-threshold units which are recruited

Identified population
CST close to the plateau of force. From Fig 2 in [47], it is worth noting that the computed CSI in Fig  6B (red curve) accounts for 60% of the variance of the true synaptic input, which is linearly transmitted by the MN pool, while the remaining variance is the MN-specific synaptic noise, which is assumed to be filtered by the MN pool and is neglected in the computation of the eND in this workflow. As discussed in the Methods section, the normalized CSI (red curve in Fig 6B) was simplified to be linearly related to the total dendritic membrane current I(t). To do so, the scaling factor was determined with a typical distribution of the MN membrane rheobase I th (j) in a cat MN pool, which was obtained (Fig 7A)  Using (Eq 6), the time-history of the current input I(t) was obtained ( Fig 7B): 3:9 � 10 À 9 þ 6:0 � 10 À 8 � CSIðtÞ else ðEq 16Þ Step (

3): LIF model -MN size calibration and distribution
Because of the modelling choices made for our MN LIF model, the MN inert period (IP) and the MN size S parameters entirely define the LIF-predictions of the MN firing behaviour. The IP i parameters of the N r MNs in the dataset DTA 35 (Fig 8) were obtained from the maximum firing frequency of the 20 MNs identified to 'saturate', from which the distribution of IP values in the entire pool of N MNs was obtained: IP(j)[s] = 0.04�j 0.05 . With this approach, the maximum firing rate assigned to the first recruited and unidentified MN is 1 The IP distributions obtained with this approach for the three other datasets are reported in the fourth column of Table 2 and yielded physiological approximations of the maximum firing rate for the unidentified lowest -threshold MN for all datasets, with the exception of the dataset HTA 50 , which lacks the information of too large a fraction of the MN pool for accurate extrapolations to be performed.
The size parameter S i of the N r LIF models was calibrated using the minimization function in (Eq 14) so that the LIF-predicted filtered discharge frequencies FDIF sim i ðtÞ of the N r MNs replicated the experimental FDIF exp i ðtÞ, displayed in Fig 9A and 9B (blue curves). As shown in Fig 9C, the recruitment time ft 1 of three quarters of the N r identified MNs was predicted with an error less than 250ms. The calibrated LIF models were also able to accurately mimic the firing behaviour of the 27 lowest-threshold MNs as experimental and LIF-predicted FIDF traces compared with r 2 >0.8 and nRMSE<15% (Fig 9D and 9E). The scaled LIF models reproduced the firing behaviour of the five highest-threshold MNs (Ni>300) with moderate accuracy, with The cloud of N r pairs fN i ; S N i g of data points (black crosses in Fig 10A), obtained from MN mapping (Fig 5B and 5C) and size calibration for the dataset DTA 35 , were least-squares fitted (red curve in Fig 10A, [42,[102][103][104][105][106][107], even if slightly in the lower range. The low value for the c coefficient observed in Table 2 for the TA 50 and GM 30 datasets suggests that the typical skewness in MN size distribution reported in the literature [42] was not captured for these two datasets, in which low-threshold MNs in the first quarter of the MN population were not identified from HDEMG signals. As displayed in Fig 9A and 9B, the phenomenological adjustment of the R and C m parameters described in (Eq 12) and (Eq 13) successfully captured the main effects of the hysteresis mechanisms involved with MN derecruitment, where MNs discharge at lower rate and stop firing at lower current input than at recruitment. In all four datasets, J C m ð Þ ¼ nRMSEÀ r 2 2 was parabolic for incremental values of C m and a global minimum was found for C d m ¼ 2:0; 2:0; 1:6 and 2.2μF�cm 2 for the datasets DTA 35 , HTA 35 , HTA 50 and HGM 30 , respectively. These values for the parameter C d m were retained for the final simulations (Step (4)) over the ½t r 3 ; t r 5 � time range.
Step (4): Simulating the MN pool firing behaviour As displayed in Fig 10B for the dataset DTA 35 , the S(j) distribution determined the MN-specific electrophysiological parameters (input resistance R and membrane capacitance C) of a cohort of N = 400 LIF models, which predicted from I(t) the spike trains of the entire pool of N MNs ( Fig 10C). As displayed in the plot of smoothed FIDFs in Fig 10C, the output behaviour of the reconstructed MN pool agrees with the Onion-skin scheme [108] where lower-threshold MNs start discharging at higher frequencies and reach higher discharge rates than larger recruited units.

Validation
The four-step approach summarized in  (2) and (3) therefore affects the quality of the predictions more than ignoring the information of MNs that are representative of a small fraction of the MN population. In all datasets, nRMSE>20% and r 2 <0.8 was mainly obtained for the last-recruited MNs (4 th quarter of each plot in Fig 12) that exhibit recruitment thresholds close to the value of the current input I(t) during the plateau of constant force in the time range ½t tr 2 ; t tr 3 �. In all datasets, the predictions obtained for all other MNs that have intermediate recruitment thresholds and are the most identified MNs in the datasets, were similar and the best among the pool of N r MNs.
Validation 2. The effective neural drive predicted by the 4-step approach summarized in Fig 1 was validated by comparing for isometric contractions, the normalized eND N (orange traces in Fig 13) computed from the N predicted spike trains sp sim j ðtÞ to the normalized force trace FðtÞ (green trace in Fig 13). The eND N was accurately predicted for the datasets DTA 35 and HTA 35 , with r 2 = 0.97−0.98 and nRMSE<8% ( Table 3). The results obtained from the dataset HGM 30 returned r 2 = 0.89 and nRMSE = 18.2% (Table 3), accurately predicting the eND N for the positive ramp and first half of the plateau of force and then underestimating the true neural drive (HGM 30 , Fig 13). The eND N predicted for the dataset HTA 50 returned results . One out of ten IDFs is displayed for clarity. The blue-to-red gradient for smallto-large MNs shows that the output behaviour of the reconstructed MN population is in agreement with the Onion-Skin scheme.
https://doi.org/10.1371/journal.pcbi.1010556.g010 of lower accuracy (r 2 = 0.88 and nRMSE = 15.1%) compared to the three other datasets. A null eND was predicted for one third of the simulation where a muscle force up to 20%MVC is generated, and the eND was overestimated for the rest of the (de)recruitment phase (HTA 50 , orange trace, Fig 12). As detailed in the discussion, the latter is explained by an inadequate IP (j) distribution (Table 2), due to a lack of experimental information in the dataset HTA 50 , which returns non-physiological maximum firing rates for the low-thresholds MNs. Two data points (16; 0.032) and (118; 0.037), obtained from the experimental data in the dataset DTA 35 with a scaling factor of 35 50 applied to the IP values, were appended to the list of N r data points ðN i ; IP N i Þ to describe the maximum firing behaviour of the first half of the MN pool for which no MN was identified for the dataset HTA 50 (Table 2). A new IP(j) distribution was obtained, returning an improved estimation of the eND N (HTA 50 , purple trace, Fig 12) with respectively lower and higher nRMSE and r 2 metrics (Table 3). With three times lower nRMSE and higher r 2 values for all four datasets (Table 3), the eND N predicted from the 4-step approach (orange dotted traces in Fig 13) was a more accurate representation of the real effective neural drive than the eND N r (blue dotted traces in Fig 13) computed from the N r experimental spike trains,  (2)). I(t) is input to N r LIF models of MN to derive, after a one-parameter calibration step minimizing the error between experimental and LIF-predicted FIDFs, the distribution of MN sizes across the complete MN pool (Step (3)). The distribution of MN sizes, which entirely describes the distribution of MN-specific electrophysiological parameters across the MN pool, scales a cohort of N = 400 LIF models which transforms I(t) into the simulated spike trains of the N MNs of the MN pool (Step (4)). The effective neural drive to muscle is estimated from the N simulated spike trains.
https://doi.org/10.1371/journal.pcbi.1010556.g011 especially in the phases of MN (de)recruitment where the real effective neural drive was underestimated by eND N r . The predicted eND N also produced less noise than the eND N r trace during the plateau of force. With r 2 >0.85 and nRMSE<20%, this four-step approach is valid for accurately reconstructing the eND to a muscle produced by a collection of N simulated MN spike trains, which were predicted from a sample of N r experimental spike trains. It is worth noting that the accuracy of the eND N predictions, performed in Table 3  MNs systematically compared to the real effective neural drive with r 2 = 0.98, and respectively returned nRMSE = 7.7, 6.9, 6.5, 6.1 and 5.9%, with slightly more accurate predictions of eND N with an increasing N at low force during derecruitment.  35 , for which the eND N was the most accurately predicted among the four datasets (Fig 13, Table 3), the time evolutions of the whole muscle force F(t) recorded experimentally (green curve) and the whole muscle forces predicted with the MN-driven muscle model described in Fig 3 using N r experimental (F N r ðtÞ, blue curves) and N simulated (F N (t), red curves) spike trains as inputs. The muscle force trace F N (t) obtained from the N simulated spike trains derived in steps (1-4) compared with the experimental force F(t) with r 2 = 0.95 and nRMSE = 11% for the dataset DTA 35 and r 2 = 0.97 and nRMSE = 11% for the dataset HTA 35 . These results validate the approach provided in this study for predicting the whole muscle force as a collection of MU force contributions from the   50 (1) were obtained with the standard approach, while those for HTA 50 (2)  MN-specific contributions to the effective neural drive. The muscle force trace F N (t) obtained from the N MN spike trains returned more accurate predictions than F N r ðtÞ (r 2 = 0.88 and nRMSE = 27% for the dataset DTA 35 and r 2 = 0.79 and nRMSE = 33% for the dataset DTA 35 ), in which case the reconstruction of the MN pool in steps (1-4) was disregarded and the N r experimental spike trains were directly used as inputs to the muscle model in Fig 3.

Discussion
This study reports a novel four-step approach, summarized in Fig 1 and displayed in detail in Fig 11, to reconstruct the recruitment and firing behaviour of a complete human pool of N MNs from a sample of N r experimental spike trains obtained from the decomposition of HDEMG or intramuscular recordings during voluntary contractions. This approach can help neuroscientists, experimentalists, and modelers to investigate MN pool neuromechanics, better understand experimental datasets, and control more detailed neuromuscular models to advance our understanding of the neural strategies taken by the human central nervous system to control voluntary muscle contractions. The three first steps of our approach identify from a sample of N r experimental spike trains a distribution of the MN electrophysiological properties across the MN pool. The N r MN spike trains are used to approximate the common synaptic input CSI(t) to the complete MN pool (Fig 7). For simplicity, the CSI is linearly related to the post-synaptic total dendritic membrane current I(t), which is input to a cohort of N r single-compartment LIF models with current synapses. The LIF firing behaviour is entirely described by the phenomenological MN Inert Period parameter IP i , derived from the experimental data (Fig 8), and by the MN size S i parameter to which all the MN electrophysiological properties are related, according to the relations provided in [42]. After calibration of the S i parameter, the N r LIF models accurately mimic the filtered discharge behaviour and accurately predict the recruitment dynamics of the N r experimental MNs (Fig 9). The N r MNs are allocated into the complete MN pool (Fig 5B  and 5C) according to their recorded force recruitment thresholds F th i ( Fig 4B) and a typical species-and muscle-specific F th (j) distribution (Fig 5). From the previous findings, the continuous distribution of MN sizes S(j) (Fig 10A) is derived for the complete pool of N MNs. S(j) defines the electrophysiological properties [42] of the MNs constituting the complete MN pool. The neural behaviour of the complete pool of N MNs is predicted in the 4 th step with a cohort of N S j -IP j -scaled LIF models and the application of the current drive I(t) (Fig 10).

Validation of the approach
The approach was successfully validated both for individual MNs and for the complete pool of N MNs. By blindly scaling the LIF models in steps (1-3) with permutations of N r −1 input experimental spike trains, the filtered discharge behaviour and the recruitment dynamics of the N r individual MNs were accurately predicted for the four investigated datasets (Fig 12). The effective neural drive (eND N ) to muscle elicited by the complete pool of N MNs was also accurately predicted by the 4-step approach for the HTA 35 , DTA 35 and HGM 30 datasets ( Fig  13, Table 3). The latter result suggests that the recruitment dynamics and the normalized distribution of firing rates across the firing MN pool were accurately predicted for the non-identified population of N−N r MNs. The accuracy of the eND N predictions was not sensitive to the size N of the reconstructed population, with acceptable eND N predictions (nRMSE<20%, r 2 >0.9) obtained with as few as ten distributed MNs. This suggests that the MN mapping (Step 1) and the derivation of the MN properties distribution (Step 3) in Figs 8 and 10A are the key contributions for describing the complete MN pool behaviour from the input experimental data.
Despite the modelling limitations, discussed in [44] and in the Limitations Section, inherent to single-compartment LIF models with current synapses and related to the phenomenological IP, C d and R d parameters, the LIF models scaled following the methodology described in Step (3) can capture the distribution of MN discharging characteristics across the MN population. For example, the firing activity of the reconstructed MN population agrees with the onionskin scheme (Fig 10D), as expected for slow low-force trajectories in human muscle contraction. While the MN size parameter calibration (Fig 9) only relies on the minimization of the root-mean square difference between experimental and predicted FIDFs, the normalized timehistory of the FIDFs (high r 2 in Fig 12) and the MN recruitment times (low Δft 1 values in Fig  12) of the N r MNs are both accurately mimicked (Fig 9) and blindly predicted (Fig 12) for all four datasets. Some realistic aspects of the MN neurophysiology are also maintained. For example, the range of MN surface areas [0.15; 0.36]mm 2 and input resistances [0.5; 4.1]MO obtained from the S(j) distribution falls into the classic range of MN properties for cats [21,42], which is consistent with the distribution of rheobase values used in Step (2) that was obtained from cat data. The distribution of MN input resistance in the MN pool, obtained from the calibrated sizes S(j), defines the individual MN rheobase thresholds and predicts that 80% of the human TA MU pool is recruited at 35%MVC (Fig 10C), which is consistent with previous findings in the literature [20,64].

Benefits of the approach
Benefit 1 -Better understanding the MN pool firing behaviour inferred from decomposed HDEMG outputs. The workflow provides a method to better understand datasets of N r experimental spike trains obtained from signal decomposition. HDEMG decomposition techniques commonly identify samples of few tens of MNs [27] at most, i.e. typically less than 10-25% of the MU pool. The mapping in Fig 5 and Table 2. EMG decomposition shows a bias towards undersampling small low-threshold MUs and identifying large high-threshold MUs because of their larger electric signal amplitude detected by the HDEMG electrodes. Consequently, the CST N r computed from the N r identified MNs is not an accurate representation of the true muscle neural drive (Fig 13, Table 3). The method in Fig 11 identifies credible values for some MNspecific electrophysiological and morphometric properties, including MN membrane surface area, input resistance, capacitance, rheobase and time constant, for the dataset of N r experimental MNs and for the N−N r MNs that cannot be recorded. Using a realistic distribution of these properties, the eND N produced by the reconstructed MN population is accurately estimated, insensitive to the value of N, because it includes the firing activity of the complete fraction of the smaller MNs that were underrepresented in the experimental dataset, especially during the force ramps in [t 1 ; t 2 ] and [t 3 ; t 4 ] (Fig 2), where the eND N r can reach 70% normalized maximum error. The eND N besides includes the activity of the complete population of highthreshold MNs, which reduces noise (orange vs blue trace in Fig 13) during the plateau of force in [t 2 ; t 3 ]. The eND N also filters the local non-physiological variations displayed by the eND N r such as the decrease in eND N r in [10; 15]s in DTA 35 and HTA 35 (top plots in Fig 13), the sudden rise and drop in eND N r at t = 10s and t = 23s in HTA 50 (bottom left plot in Fig 13) or the steeply decreasing eND N r from t = 28s in HGM 30 (bottom right plot in Fig 13). Understanding and accounting for these experimental limitations is important and has critical implications in HDEMG-driven neuromuscular modelling when the CST N r is input to a muscle model [109,110].
The workflow also provides a method to support future experimental investigations. The signal decomposition process can be refined and more MNs can be identified by iteratively running this workflow and using the N−N r simulated spike trains of the non-identified MNs as a model for identifying more MNs from the recorded signals. The method moreover provides simple phenomenological two-parameter-scaled LIF models described in step (3) to accurately replicate the recruitment and firing behaviour of the N r identified MNs (Fig 9) for further investigation.
Benefit 2 -Advances in MN pool modelling. This four-step approach advances the state-of-the-art in MN pool modelling. As discussed in the introduction, using a sample of N r experimental MN spike trains permits for the first time in human MN pool modelling to (1) estimate the true time-course of the common synaptic input to the MN pool, which cannot be measured experimentally, (2) approximate with a one-parameter calibration, available HDEMG data and knowledge from the literature, the MN-specific electrophysiological properties of all the MNs in the MN pool and a realistic distribution of these properties across the pool, which could, to date, only be obtained in specialized experimental and MN pool modelling studies in animals [39], and (3) validate the forward predictions of MN spike trains and effective neural drive to muscle to human experimental data. The pool of LIF models, the maximum firing frequency of which is obtained from the available experimental data, intrinsically accounts for the onion-skin theory [108] (Fig 10D), and better replicates the MN membrane dynamics of the MN pool than Fuglevand-type phenomenological models [10,29], where the MN-specific firing rates are predicted as a linear function of the current input I(t). Moreover, single-compartment current-input LIF models can credibly replicate most of the MN membrane dynamics [43] and allowed accurate predictions of the MN pool behaviour (Figs 12 and  13), while they provide a simpler modelling approach and a more convenient framework for MN electrophysiological parameter assignment in the MN pool than comprehensive Hodgkin-Huxley-type MN models [30,32,111,112]. This four-step approach demonstrated to be robust to systematic differences in the input experimental datasets. For example, it accurately predicted the individual MN firing behaviour of the N r MNs for all four datasets (Fig 12) despite the latter being obtained at different levels of muscle contractions, on different subjects and muscles and in different experimental approaches. The approach accurately predicted the eND N of both DTA 35 and HTA 35 datasets from N r = 32 and 21 identified MNs respectively ( Table 3), suggesting that the method is not sensitive to the number N r of identified MNs, providing that the N r identified MNs and their properties are homogeneously spread in the MN pool, as in datasets DTA 35 and HTA 35 where at least one MN is identified in the each 10%range of the rheobase domain. The accurate prediction of eND N for the dataset HGM 30 in the time range [4.7; 7.2] s (Fig 13 bottom-right plot) also suggests that the quality of the predictions is not sensitive to the hindlimb muscle investigated, providing that the F th distribution is representative of the investigated muscle.
Benefit 3 -Relevance for neuromuscular modelling. A MN-driven neuromuscular model of in-parallel Hill-type actuators describing the MUs (Fig 3) is described in the Methods section and is used to predict the experimental muscle isometric force (green trace in Fig 14) from the N r experimental or N reconstructed MN spike trains. The complete reconstruction of the discharging MN population detailed in Fig 11 is a key step towards advancing the state-ofthe-art in neuromuscular modelling on several aspects.
Firstly, the neural input to the neuromuscular model in Fig 3 is a vector of experimental spike trains that is a more easily interpretable and a more detailed description of the muscle neural drive than the rectified-normalized-filtered envelopes of recorded bipolar EMGs (bEMGs) [95,113,114] or the CST N r [109,110,115] typically used in single-input muscles models. Secondly, despite advanced multi-scale approaches in modelling whole muscles as a single equivalent MUs [98,[116][117][118][119][120], the multiple Hill-MU model in Fig 3 provides a more convenient framework to model in detail the continuous distribution of the MU excitation-contraction and force properties in the MU pool. Thirdly, a few other multi-MU models were proposed in the literature [9,10,12,39,98,121,122], some of which used the Fuglevand's formalism [29] to model the MN firing behaviour and recruitment dynamics of complete pools of N MNs, which are intrinsic to the experimental sp exp i ðtÞ, to predict whole muscle forces with collections of N MU Hill-type [9,123] Hill-Huxley-type [10] or twitch-type [39] muscle models. However, these studies lacked experimental neural control at the MN level and considered artificial synaptic inputs to their model, the predicted force of which was indirectly validated against results from other models and not against synchronously recorded experimental data, as performed in Figs 12-14 in this study. Finally, Fig 14 demonstrates that the reconstruction of the complete MN population described in this study (steps 1 to 4 in Fig 1) is a key step for accurate MN-driven neuromuscular predictions of muscle force. The N r -MU model, that receives the N r individual MN spike trains sp exp i ðtÞ, intrinsically underestimates the whole muscle activity when dominantly low-threshold MUs are recruited but are under-represented in the experimental sample (Figs 5 and 13). The N-MU model, which receives as inputs the N spike trains sp sim j ðtÞ generated by the four-step approach, allows a more realistic assignment of distributed MU properties (MU type and maximum isometric force) to the complete MU population, and returned more accurate force predictions than the N r -MU model (Fig 14). It is worth noting that the N-model did not require any parameter calibration except the MN size in step 3. The detailed modelling of the distribution of the excitation-contraction properties of the MUs makes the N-model more suitable for investigating the muscle neuromechanics than typical EMG-driven models, the neural parameters of which do not have a direct physiological correspondence and must be calibrated to match experimental joint torques [95,123,124].

Limitations of the approach and potential improvements
Despite the aforementioned good performance of the four-step workflow, the method presents 4 main limitations, for some of which potential improvements are proposed in the following.

Limitation 1 -Model validation.
The two validations of the approach (Figs 12 and 13) present some limitations. The local validation in Fig 12 only ensures that the method accurately estimates the MN firing behaviour for the fractions of the MN pool that were experimentally identified. This local validation alone does not inform on the accuracy of the predictions for the non-identified regions of the MN pool and must be coupled with a global validation of the MN pool behaviour by validating the predicted eND N . While the local validation was successful for the HTA 50 dataset (Fig 12), where less than 30% of the active MN pool is represented, it is inferred in Fig 13C that the individual MN firing rates were overestimated for the first half of the non-recorded MN pool. This local validation would be self-sufficient for experimental samples that contain a large and homogeneously spread population of identified MNs, obtained from decomposed fine-wire intramuscular electrode and HDEMG grid signals. The validation of the eND N is performed in this study against an experimental force recorded by a transducer (green traces in Fig 13), which accounts for the force-generating activity of both the muscle of interest and the synergistic and antagonistic muscles crossing the ankle. The experimental force trace measured in the TA datasets may be an acceptable validation metric as the TA muscle is expected to explain most of the recorded ankle torque in dorsiflexion. The TA muscle is indeed the main dorsiflexor muscle with a muscle volume and maximum isometric force larger than the flexor hallucis longus muscle, which moment arm is besides not aligned with the dorsiflexion direction and mainly acts for ankle inversion. However, the gastrocnemius lateralis, soleus, and peroneus muscles acts synergistically with the GM muscle for ankle plantarflexion and the experimental torque recorded in dataset HGM 30 may not be representative of the GM muscle force-generation activity and may not be a suitable validation metric for eND N . In Fig 13, the decreasing eND N and eND N r (orange and blue dotted traces) may accurately describe a gradually increasing sharing of the ankle torque between synergistic muscles initially taken by the GM muscle, which the constant validation trace does not capture. To answer such limitations, the predicted eND N should be validated against a direct measure of muscle force, which can be performed as in other recent studies [13,14,116] with ultrasound measurements of the muscle fascicle or of the muscle tendon concurrently obtained with (HD)EMG recordings of the muscle activity. Finally, these two validations do not provide any indication whether the parameters calibrated with a trapezoidal force profile would generalize, for the same subject, to another force trajectory or a trapezoidal force profile up to another force level and provide accurate predictions of the N r experimental spike trains and of the eND N . To perform such validation, the parameter identification in Step 3 ( Fig  1)  Limitation 2 -Sensitivity of the method to input HDEMG data. While the method predicts a list of simulated spike trains sp sim j ðtÞ and a eND N that more accurately describes the MN pool behaviour than the experimental sp exp i ðtÞ and eND N r , as discussed previously, the accuracy of these predictions (Figs 12 and 13) remains sensitive to the distribution in the entire MN pool of the N r MNs identified experimentally, reported in the third column of Table 2. Because of the definition of the current input I(t), the eND N onset is defined by the recruitment time tf 1 N 1 of the lowest-threshold MN N 1 identified experimentally, as shown in Fig 13, while the unknown firing behaviour for t < ft 1 N 1 of the non-identified MNs of rheobase lower than I th N 1 is not captured. This is not an issue in samples of homogeneously distributed MNs like dataset DTA 35 , where the 9 th smallest MN (first 2% of the pool of threshold-increasing MNs) is identified (Fig 5BC), Table 2) and the eND obtained from the 15 first MNs is not predicted during the short time range [2.1, 2.3]s, where the whole muscle builds only 1%MVC (top-left plot in Fig 13). However in heterogeneous or incomplete samples like dataset HTA 50 , the lowestthreshold (232 nd , Table 2) identified MN identified is recruited at t ¼ tf 1 N 1 ¼ 6:1s and the approach wrongly returns a zero muscle neural drive eND = 0 for the time period [1.6; 6.1]s where the muscle actually builds 20%MVC during the ramp of contraction (bottom-left plot in Fig 13). To tackle this issue, the normalized force trace, which is non-zero for t < ft 1 N 1 and representative of the CSI in this time region, could be scaled and used in lieu of the current definition of I(t). However, this approach, which may not be suitable for non-isometric conditions, is not applicable in forward predictions of unknown whole muscle force from neural inputs. Experimental samples of homogeneously distributed MNs are also required to derive realistic S(j) and IP(j) distributions with the four-step approach. As observed in Fig 12, ignoring in the derivation of I(t), IP(j) and S(j) the spike trains sp exp i ðtÞ of the MNs that are representative of a large subset of the MN pool affects the accuracy of the predicted recruitment and firing behaviour of the MNs falling in that subset (Fig 12). More specifically, the non-physiological distribution IP(j)[s] = 5.6�10 −4 �j 0.80 was fitted to the experimental data of the 14 high-threshold MNs identified in dataset HTA 50 , where the neural information of the 60% smallest MNs of the MN pool (Table 2) is lacking. Such IP(j) distribution overestimates the LIF-predicted maximum firing frequency of low-threshold MNs, which explains the overestimation of the eND N in Fig 13  (bottom-left plot, orange curve). Non-physiological predictions can be avoided by adding artificial data points consistent with other experiments or with the literature in the rheobase regions of the MN pool where no MNs were experimentally identified. For example, using an IP(j) relationship consistent with a dataset of homogeneously distributed MNs (DTA 35 ) constrained the predicted maximum firing rates to physiological values and returned more accurate predictions of the eND N (bottom-left plot, purple curve).
Limitation 3 -LIF MN modelling. The LIF MN model described from (Eq 7) to (Eq 13) was shown, with credible sets of inter-related parameters S MN , R, C and τ after [42], to accurately mimic (Fig 9) and blindly predict (Fig 12) most of the key phenomenological features of the N r firing MNs, including their FIDFs and time of first discharge, as well as nonlinear behaviours such as firing rate saturation (Fig 8) and hysteresis (Fig 4B and 4C). However, this MN model is mostly phenomenological and does not provide a realistic description of the actual mechanisms underlying the dynamics of action potential firing for several reasons. First, this single-compartment approach neglects the activity of the MN dendrites, which account for more than 95% of the total MN surface area and gather most of the post-synaptic receptors and MN PIC-generating channels responsible for the MN nonlinear input-output functions [40]. The inherent difference in membrane voltage dynamics between the dendrites and the soma, also partially mediated by hundred-fold differences in the value of passive electrophysiological parameters such as R which increases with somatofugal distance [125] is therefore neglected in this point model approach. Then, the nonlinear dendritic activity being overlooked, the Common Synaptic Input CSI(t) in (Eq 5), which is the common net excitatory synaptic influx, is non-physiologically assimilated with a constant gain to the post-synaptic total dendritic membrane current I(t) that depolarizes the MN soma and is responsible for spike generation. With this linear CSI−I scaling in (Eq 6), the approach neglects the realistic description of the voltage-driven dynamics of the PIC-generating channels that are responsible for some important nonlinear mechanisms, such as firing rate saturation [41] and hysteresis in the MN's current-voltage relation [93,94], which dictate the nonlinear CSI−I transformation [40]. In this study, the firing rate saturation, which in real MNs is due to a decrease in driving force for synaptic current flow as the dendrites become more depolarized [41], is captured by the phenomenological IP parameter. The observed hysteresis between the recruitment and decruitment thresholds (Fig 4B) and in the current-firing rate relation between recruitment and derecruitment phases (Fig 4C), explained by an hysteresis in the MN's current-voltage relation that arises from the activity of long-lasting PICs [20], are successfully but non-physiologically addressed by a tuning during derecruitment of the R and C m parameters respectively, although the values of the passive property R and of C m , biological constant for the MN population [126], should be independent from the MN activity [44]. For these reasons, even if the distributions of the S MN , R, C and τ parameters at recruitment take values consistent with the literature, as discussed, the MN LIF model described in (Eq 7) to (Eq 12) is not suitable for investigating the neurophysiology and the neuromechanics of individual MNs. The phenomenological tuning of the IP, R and C m parameters, although symptomatic of underlying nonlinear mechanisms, does not provide any insight into the true MN neurophysiology. If a more realistic MN model was considered in lieu of the LIF model, such as multiple-compartments Hodgkin-Huxley-based approaches [127], the four-step workflow would be an adequate tool for testing neurophysiological hypotheses and investigating some mechanisms and properties of the complete MN pool that cannot be directly observed experimentally in conditions of human voluntary contractions. A possible trade-off would be considering a single-compartment LIF model with active conductances [44], in which case (Eq 7) is re-written as (Eq 18), where g R ¼ 1 R is the constant MN passive conductance, E + >ΔV th is the constant reversal potential for excitatory channels, g(t, CSI) the time-varying total active excitatory synaptic conductance, and I ps (t, CSI, V m ) = g(t, CSI)�(E + −V m ) the post-synaptic total dendritic membrane current induced by the synaptic drive CSI, the driving force of which is decreased by the term (E + −V m ) when the dendritic membrane depolarizes, as discussed.
In the literature, g(t, CSI) = g max �T(t, CSI) where g max is the maximum active conductance of the synapse and T(t, CSI) can be interpreted as the fraction of bound synaptic receptors [30] or of opened ionic channels in the range [0; 1]or as a train of synaptic pulses [44]. Because setting T(t, CSI) = CSI(t), with CSI(t) as defined in Fig 6B, does not meet the requirements in firing rate saturation for large CSI input, the saturation in dendritic ionic channel activation for large CSI [41] must be accounted for and T(t, CSI) should be a nonlinear saturating function of CSI. In this approach, the MN-specific hysteresis and firing rate behaviour could be obtained with a tuning during recruitment and derecruitment and a distribution across the MN pool of the g max parameter and of the shape parameter of the S function based on the experimental data. Although this approach more realistically describes the neurophysiological mechanisms underlying the nonlinear MN behaviour during discharge events, it is more challenging to scale using the experimental information. Considering that this study focuses on the phenomenological behaviour of individual MNs and on the overall activity of reconstructed MN populations, the LIF model defined in (Eq 7) was judged an adequate trade-off between accuracy and modelling complexity, considering the overall accurate predictions of the MN firing behaviours (Fig 12) and its low computational cost, and no other modelling approaches such as (Eq 18) were pursued.
It must be noted that, whatever the MN modelling approach, the experimental drive to the MN currently considered in the four-step method is the CSI, which disregards the MN-specific synaptic noise SN(t), which is responsible for most of the inter-spike variability (ISV). Considering that the MN pool and the MU neuromechanical mechanisms are expected to filter the MN-specific SN(t), this simplification is adequate for accurate predictions of the eND N and of the whole muscle force [83]. Sometimes modelled as a random Gaussian-like event [29], the ISV can be obtained as the response to a SN(t) white noise added to the current definition of the synaptic input as CSI(t)+SN(t). The relative proportion of CSI and SN in % of the variance of the total synaptic current can be approximated from Fig 2A in [47]. Accounting for SN(t) might improve the accuracy of the predicted firing behaviour of the largest MNs, for which the largest Δft 1 and nRMSE and lowest r 2 values were obtained in Fig 12 for all four datasets, and the firing behaviour and recruitment dynamics of which are dominantly dictated by random fluctuations of SN(t).
Limitation 4 -Limited experimental data in the literature. The four-step approach is constrained by the limited knowledge in the literature of the characteristics of the human MU pool. Therefore, the LIF parameters and the MN rheobase distribution in the MN pool ( Fig  7A) are tuned with typical cat data from various hindlimb muscles [42] while the experimental MN spike trains were obtained from the human TA muscle. While the normalized mathematical relationships relating the MN electrophysiological parameters of the LIF model (R, C, τ, ΔV th ) can be extrapolated between mammals [42] and thus to humans, no experimental data is yet available to scale these relationships to typical human data.
Also, the typical F th (j) and I th (j) threshold distributions, derived for mapping the N r identified MNs to the complete MN pool (Fig 5) and for scaling the CSI to physiological values of I (t) (Fig 7), were obtained from studies which relied on experimental samples of MN populations of small size. These source studies therefore did not ensure that the largest and lowest values were identified and reported, and that the identification process was not biased towards a specific subset of the MN pool, such as larger MNs. In such cases, the true threshold distributions would be more skewed and spread over a larger range of values, as discussed in [42], than the distributions reported in Figs 5 and 7. The F th (j) distribution is besides muscle-specific [20] with large hindlimb muscles being for example recruited over a larger range of MVC than hand muscles. However, enough data is reported in the literature to build the F th (j) distributions for the TA and first dorsal interossei human muscles only. For these limitations, the two first steps of this approach could be made subject-species-and muscle-specific by calibrating the F th (j) and I th (j) described as the 3-parameter power functions defined in this study.
Conclusion. This study presents a four-step workflow (Figs 1 and 11) which predicts the spiking activity of the discharging MNs that were not identified by decomposed HDEMG signals. The method is driven by the common synaptic input, which is derived from the experimental data, and reconstructs, after a calibration step, the distribution across the MN population of some MN properties involved into the MN-specific recruitment and spiking behaviour of the discharging MNs (Figs 8 and 10A). The method blindly predicts the discharge behaviour of the N r experimentally identified MNs (Fig 12) and accurately predicts the muscle neural drive (Fig 13) after the complete discharging MN population is reconstructed (Fig 10C and 10D). With direct application in neuromuscular modelling (Fig 14), this method addresses the bias of HDEMG identification towards high-threshold large units and is relevant for neuroscientists, modelers and experimenters to investigate the MN pool dynamics during force generation.