Person-Specific Biophysical Modeling of Alpha-Motoneuron Pools Driven by in vivo Decoded Neural Synaptic Input

Interfacing with alpha-motoneurons (MNs) is key to understand and control motor impairment and neurorehabilitation technologies. Depending on the neurophysiological condition of each individual, MN pools exhibit distinct neuro-anatomical properties and firing behaviors. Hence, the ability to assess subject-specific characteristics of MN pools is essential for unravelling the neural mechanisms and adaptations underlying motor control, both in healthy and impaired individuals. However, measuring in vivo the properties of complete human MN pools remains an open challenge. Therefore, this work proposes a novel approach based on decoding neural discharges from human MNs in vivo for driving the metaheuristic optimization of biophysically realistic MN models. First, we show that this framework provides subject-specific estimates of MN pool properties from the tibialis anterior muscle on five healthy individuals. Second, we propose a methodology to create complete pools of in silico MNs for each subject. Lastly, we show that neural-data driven complete in silico MN pools reproduce in vivo MN firing characteristics and muscle activation profiles during force-tracking tasks involving isometric ankle dorsi-flexion, at different levels of amplitude. This approach can open new avenues for understanding human neuro-mechanics and, particularly, MN pool dynamics, in a person-specific way. Thereby enabling the development of personalized neurorehabilitation and motor restoring technologies.


I. INTRODUCTION
A LPHA-MOTONEURONS (MNs) are regarded as the final common pathway of central nervous system's motor and sensory pathways. Therefore, getting insight into MN behavior in vivo would be key for unravelling the neural mechanisms underlying motor control, both in healthy and neurologically impaired individuals. Depending on the specific neuro-physiological characteristics of an individual, including age [1], level of training [2], [3], severity of motor disorder [4], [5] and neural injury [6], pools of MNs exhibit distinctive neuro-anatomical properties and dynamics. For this reason, the ability to assess subject-specific MN pool characteristics is essential to understand motor impairment and enable the development of tailored neurorehabilitation technologies.
Signal-based approaches relying on high-density electromyography (HD-EMG) decomposition [7], [8], [9] enable identifying the firing output of MNs in vivo [10], [11]. These allowed the estimation of motor unit properties as innervation zone and conduction velocity of muscle fibers [8]. However, there is currently no existing method for characterizing the electrophysiological properties of complete MN pools in vivo.
In this context, biophysical models based on Hodgkin and Huxley's description [12] enable establishing a link between a neuron's firing response to an arbitrary input and its electrophysiological properties (i.e., membrane resistance, capacitance, voltage threshold, conductance of ionic channels, etc.). Previous implementations of biophysical MN models in motor control [13], [14] investigated the neural mechanisms of force steadiness [15], postural balance [16] and voluntary movement [17], suggesting that a pool of MNs receive a common synaptic input (CSI) that is linearly transformed into the neural drive to muscle [17]. However, these computational approaches relied on neuronal data derived from animal preparations [18], [19]. As such, current models are not able to capture electrophysiological differences among human individuals, nor to reproduce experimental firing patterns of human MNs in vivo during the execution of a motor task, thus limiting the clinical implications of these findings.
Metaheuristic optimization [20], [21] has been proposed for identifying neuron model parameters that reproduce experimental spike trains, including ionic properties of cortical interneurons [22] and cerebellar granular neurons [23]. Despite successfully characterizing neuron dynamics, these methods have only been implemented in vitro, where firing responses to controlled input currents injected into the soma can be recorded.
This work proposes merging HD-EMG decomposition, biophysical modelling and metaheuristic optimization to create This work is licensed under a Creative Commons Attribution 4.0 License. For more information, see https://creativecommons.org/licenses/by/4.0/ Fig. 1. Overview of neural-data driven framework for generating in silico motoneuron (MN) pools: in vivo spike trains are decoded from highdensity electromyography recordings via convolution kernel compensation blind source separation decomposition [9]. A subject-specific excitability constant (∆IF) is determined using multiple objective optimization, and the common synaptic input (CSI) received by the MN pool is computed as the product of ∆IF and the neural drive derived from in vivo MNs. Parameter optimization of in silico MNs driven by the CSI is then performed to minimize recruitment time error and frequency-corrected temporal spike match. Based on probability density functions of the identified parameters, additional in silico MNs are generated to complete an in silico pool of 200 MNs. a novel neural-data driven framework. This approach consists of decoding the activity of human MNs in vivo via HD-EMG decomposition, and using it to drive the metaheuristic optimization of biophysically realistic MN models. The outcome are physiologically realistic estimates of complete MN pool properties (i.e. distributions of soma sizes, ionic channel dynamics, pool excitability) from specific subject's muscles, which subsequently enable the creation of in silico MN pools that retain the same statistical behavior of their in vivo counterparts. Our results evidence the ability of this method for predicting in vivo MN firing characteristics and muscle activation profiles throughout isometric force-tracking tasks at varying levels of amplitude from the tibialis anterior (TA) muscle on a group of five healthy individuals.
This approach opens new avenues for understanding the neuronal mechanisms underlying human movement in a person-specific way, thereby bridging the gap limiting the development of personalized neurorehabilitation therapies and motor restoring technologies.

II. METHODS
This framework consists of three main components ( Fig. 1): 1) approximation of in vivo CSI from HD-EMG decomposed MN spike trains, 2) parameter optimization of in silico MN models driven by in vivo CSI to match recruitment time and firing frequency of in vivo MNs, 3) generation of extended in silico MN pools following subject-specific experimental distributions of biophysical properties. These steps will be described along study procedures from section II-C. onwards.

A. Participants Overview
Five healthy subjects were recruited (age: 27.4 ± 2.07 years, weight: 70 ± 12.34 kg, height: 173.6 ± 10.06 cm), who provided a written informed consent to participate in this experiment. All procedures were executed in compliance with the medical ethical committee at Radboudumc, under the research protocol NL73230.091.20.

B. Experimental Setup
After localization of TA's muscle belly via palpation, the skin was shaved and lightly abraded using abrasive paste (Meditec,Pharma). An 8 × 8 electrode grid with conductive gel was placed as described previously [24]. Subjects seated on a Biodex chair (M4 Biodex Medical Systems Inc., Shirley, NY, USA), with the right leg fixed with a knee joint angle of 140 • flexion and the right foot tightly strapped to the dynamometer keeping the ankle joint angle in a neutral position. Isometric ankle dorsi-flexion torque was acquired at a sampling frequency of 512 Hz using National Instruments Data Acquisition card (NI DAQ), while HD-EMG data was simultaneously recorded via an in-house developed data acquisition interface using TMSi Refa multichannel amplifier (TMS International B. V., Oldenzaal, The Netherlands) at a sampling frequency of 2048 Hz.
To determine the maximum voluntary contraction (MVC), subjects were asked to apply maximal dorsi-flexion force for a period of 5 s. This was repeated three times, with a resting time of 1 min between each trial.

C. Study Protocol
Subjects performed a ramp-and-hold task at five different target % MVC amplitudes: 10, 20, 30, 40, and 50% MVC. The holding period of the ramp was set to 5s, and the speed of the ascending-descending portions was set to 20% MVC/s. Throughout the experiment, all five target %MVCs were sorted randomly and five consecutive repetitions, for each amplitude, were presented to the subjects through a visual interface providing real-time feedback of the force-tracking performance. To avoid the effects of fatigue, subjects had a resting time of 5, 10, 15, 20 or 25s between each of the five repetitions, depending on if the current target %MVC amplitude was 10, 20, 30, 40 or 50%, respectively. Additionally, there was a 120s resting period between the end of the five repetitions and the next target %MVC amplitude.
Prior to recording, subjects underwent a learning trial to familiarize with the force-tracking task: two repetitions for each target %MVC presented randomly with the same resting rules.

D. High-Density Electromyography Decomposition
HD-EMGs were band-pass filtered with a 20 to 500 Hz zero-phase Butterworth filter, and decomposed into individual spike trains using convolution kernel compensation blind source separation [9]. Each decomposed experimental spike train was as a binary vector where '1' represented a spike event and '0' meant no firing activity. Decomposed spike trains underwent a quality-control algorithm (QC-MN), as previously proposed [24], to ensure pulse-to-noise ratio (PNR) > 20 dB, coefficient of variation < 0.3, and discharge rates < 30 Hz [24]. These thresholds were chosen to remove non-physiological spike trains while keeping a significative amount per trial. MNs not fulfilling this criterion were excluded.

E. Biophysical Motoneuron Model
MNs were described as single-compartment conductancebased models including leakage (g l ) sodium (g N a ), slow (g K s ) and fast potassium (g K f ) channels, which have been identified to reflect experimental neuronal behavior [13]. Ionic current was determined according to (1), where E N a and E K correspond to the reverse potentials of sodium and potassium, respectively. Voltage-dependent rate constants (m, h, n, q) were computed using a pulse-based model [25], which describes action potentials as rectangular pulses and defines each rate constant based on parameters beta (β i ) and alpha (α i ), where subscript i denotes each rate constant [25].
MN's leakage conductance (2) and capacitance (C) (3) were adjusted by the soma diameter (D s ), where R m and C m are membrane specific resistance and capacitance, respectively. Membrane voltage (V m ) was then computed according to (4), where E L represents leakage Nernst voltage and I in j the current injected in the soma.

F. Common Synaptic Input Estimation
The in vivo CSI current (in A) was computed as the product between the net discharge rate of the MN pool (in Hertz) and a subject-specific gain ( I F) . This way, in vivo CSIs enabled recruiting in silico MNs in a range of D s matching experimental measurements (i.e., 15 < D s < 220 µm) [26]. The net discharge rate of the MN pool (i.e., neural drive) was derived as the low-pass filtered sum of trialspecific MN spike trains (i.e., cumulative spike train filtered by 200 ms Hann window) [17], normalized by the average discharge rate of summated MNs. I F was determined via multi-objective optimization of recruitment error (5) in two conditions: a) earliest recruited in vivo MN and in silico MN with lowest soma diameter, and b) latest recruited in vivo MN and in silico MN with highest soma diameter. All remaining parameters were kept as previously proposed [14], since prior work [27] demonstrated dominance of D s over recruitment time. Optimization was executed using MATLAB's genetic algorithm (The MathWorks, Inc., Natick, MA, USA) for a Gaussian distributed initial population size of 200 I Fs in the range of 0.1 to 1, an elite percentage of 10%, 70% cross-over, and a function tolerance of 0.01. Gray relational analysis [28] was subsequently applied on the resulting pareto set to automatically determine the I F that best minimized both error functions. This process was performed on every participant, and each subject-specific I F was kept constant throughout the entire study.

G. Motoneuron Parameters Identification
MN parameter identification was performed on one trial per target %MVC (Section II-B), which yielded the highest number of decomposed MNs. The remaining trials were kept for validation of force-profile tracking in silico MNs (Section II-I).
in silico MN models were created and paired to each in vivo decomposed MN. Subsequently, the in silico MNs were driven by the CSI derived from each corresponding target % MVC trial (see section II-F). We implemented a double singleobjective optimization approach for robust sampling of the optimal pareto space without increasing computational load [29]. First, relative recruitment error (5) between in silico and in vivo MNs was minimized by optimizing D s within the anatomical range of 15 to 220µm [26]. Subsequently, spikematch (i.e. gamma-factor [30]) error (6), was minimized by the optimization of β Q and α Q , as early work [27] indicated slowpotassium dynamics to largely determine MNs' firing features. 2 2 where f e and f m represent the mean firing rate, N e and N m the number of spike events, and spikes e and spikes m the spike trains of experimental (i.e. in vivo) and model (i.e. in silico) MNs, respectively. N c is the number of coincident spikes within the time window δ = 2ms [30]. Optimization of both objective functions implemented MATLAB's genetic algorithm set to a population size of 200, elite percentage of 20%, cross-over of 70%, and function tolerance of 1e −5 . Furthermore, the initial populations for every parameter were generated randomly with uniform distributions. To constrain the framework for providing physiologically realistic solutions, the lower and upper boundaries of β Q and α Q , as well as non-optimized parameters of the model were kept the same as proposed previously [14].
After optimization, identified sets of MN parameters from all %MVC amplitudes were merged together, and the probability density function (PDF) of each parameter was computed using boundary-corrected kernel density estimator [31].

H. Generation of Complete in Silico Motoneuron Pools
As indicated by Duchateau and Enoka [32], MN pool size estimations of TA taken from human cadavers (i.e., 445 MNs) may be overestimated, as electrophysiological approaches [33] indicate a mean size of 200 ± 61 MNs. Since the usual number of identifiable MNs via HD-EMG decomposition is limited to few tens of MNs at forces up to 50% [7], [8] and 70% MVC [3], additional in silico MNs based on subjectspecific parameter's PDFs were generated to create customized pools of 200 MNs for every participant. First, D s values were randomly generated from the corresponding PDF. Following previously proposed size-based ranges of β Q and α Q for different MN-types [14], k-means clustering [34] was subsequently performed to identify three clusters of parameters, thereby defining cluster-specific PDFs. For each generated D s , the corresponding cluster was identified computing the minimal distance to their centroids, and the values of β Q and α Q were generated according to the cluster-specific PDF.
PDF-derived in silico MN underwent a pool quality-control algorithm to exclude sets of parameters resulting in unrealistic discharge rates (i.e., f m > max( f e ) + 3Hz). This consisted of iterative simulations of MNs driven by the CSI derived from the highest %MVC amplitude (i.e., where higher discharge rates can be found). For each iteration, PDF-derived in silico MNs not fulfilling discharge rate criteria were considered nonrealistic and discarded, the rest were approved and kept in the pool. New sets of MN parameters were created following the same rules, and the same procedure was repeated until all in silico MNs in the pool fulfilled the discharge rate constrain.

I. Force-Tracking in Silico Motoneuron Pools
The relation between force produced by TA's motor units and the percentage of recruited MNs in the pool has been described by an exponential fit [32]. As the force required to reach 100% MN recruitment in TA is reportedly decreased at high rates of force development (i.e., from ∼90% to ∼50% MVC) [35], [36] we modified the fit function [32] accordingly (7) and implemented it to control the percentage of MN recruitment as a function of % MVC throughout our simulations.

J. Validation Procedures
Test 1 assessed whether optimization of proposed objective functions enabled finding unique solutions for creating in silico MN copies that fire similarly to in vivo MNs. For this, we quantified the costs values of objective functions (5) and (6) before optimization (i.e. using a set of initial parameters D s 0 , β Q 0 , α Q 0 = [30, 0.03, 1.5] corresponding to an average small size MN, reportedly dominant in TA [37] and traditionally associated to slow-type motor units [14]), and after optimization. Additionally, we tested the repeatability of the solutions on a representative set of 20 in vivo MNs derived from a 50% MVC trial (i.e., including MNs from all recruitment thresholds). For each MN, we performed ten consecutive optimizations and measured the standard deviation (STD) of the resulting parameters.
Test 2 evaluated inter-subject variability in the identified MNs PDFs. Assessed metrics included mean and STD among MN parameters and cluster thresholds.
To validate whether in silico MNs driven by a trialspecific CSI produced firing outputs of similar characteristics as their in vivo counterparts, test 3 quantified absolute error in recruitment time and mean discharge rate (during plateau) between in vivo and in silico MNs. This was performed on two conditions: before and after parameter optimization.
Test 4 assessed the ability for creating extended in silico MN pools from subject-specific parameter's PDFs. This analysis included quantification at intra-subject level of median and STD, as well as visual PDF inspection, of the parameters describing optimized in silico MNs (i.e., the total number of decomposed MNs) and extended in silico pools including PDF-derived in silico MNs (i.e., 200 MNs).
Test 5 evaluated whether the use of extended in silico MN pools driven by trial-specific CSIs improved the estimates of force profiles (i.e., neural drive) throughout all conditions. Two metrics were used for this: 1) Coefficient of determination (R 2 ) between measured torque and neural drive. 2) Coefficient of variation (CoV) of neural drive. Statistical differences between R 2 and CoV of in vivo MNs, their corresponding in silico copies, and the extended in silico MN pool (including PDF-derived in silico MNs) were tested using non-parametric Friedman test.

III. RESULTS
After HD-EMG decomposition and QC-MN, a total of 413 MNs (PNR = 31.11 ± 7.2) from all subjects and conditions were included in the study. MN count per subject is shown in Fig. 3.
For all participants, the determined I F values were similar (Table I), with an inter-subject mean I F = 2.049 ± 0.195 H z A . Test 1: After optimization, normalized cost values of both objective functions (5) and (6) decreased by one order of Fig. 3. Subject-specific in vivo motoneuron (MN) properties identified by the optimization framework. Each MN is described by three parameters: beta (β Q ), alpha (α Q ) and soma diameter (D s ). A) K-means clustering consistently identified three clusters of parameters across all subject. These clusters were characterized for having a diameter threshold of D s = 79.44 ± 5.28 µm between C 1 and C 2 , and D s = 133.88 ± 2.54 µm between C 2 and C 3 . The three clusters are color-coded in blue, red and green, corresponding to C 1 , C 2 and C 3 . Additionally, the percentage of each cluster within the MN pool is depicted on the upper right corner of each subject's specific parameter characterization. Lastly, we show a pie chart depicting the averaged cluster proportions and D s threshold from all subjects. B) Each column shows the subject-specific distribution found for each MN parameter.

SUBJECT-SPECIFIC EXCITABILITY CONSTANTS (∆IF) Hz
A magnitude for all subjects (Fig. 2). For the recruitment time error function (5), cost values approached zero in all cases. Repeatability test (Fig. 2) showed no difference in estimated D s for any MN after consecutive optimizations (average STD = 0.018). For β Q and α Q , the average STDs were 1.511 and 1.212, respectively.
Test 2: The inter-subject mean value of D S was 113.53 ± 5.2 µm. Throughout all subjects, three clusters (C i ) of MN parameters were identified (Fig 3.A), with both β Q and α Q showing increasing amplitude and spread as D s increases. All five subjects showed D s thresholds separating three clusters of parameters. The inter-subject threshold between C 1 and C 2 was found at D s = 79.44 ± 5.28 µm, whereas the threshold between C 2 and C 3 was D s = 133.88 ± 2.54 µm. Based on to these cluster thresholds, the proportions of MN sizes found in the pool included: C 1 = 28.6 ± 5.59%, C 2 = 39.6 ± 4.72%, and C 3 = 31 ± 4.85% (Fig. 3.B). Similarities in inter-subject PDFs were found for the three MN parameters (Fig. 3.C), with the largest differences found in D s , where seemingly multimodal distributions can be observed.
Test 3: Fig. 4 shows a representative instance of the spike trains produced by in silico MNs (from subject 2) before and after optimization. Recruitment time of in silico MNs before optimization remained unaltered for all cases (Fig. 4.A), and the net discharge rate exhibited a limited range that varied from 6.78 to 7.42 Hz as %MVC increased (Fig. 4.C). In contrast, in silico MNs after optimization matched recruitment time of in vivo MNs with an error of < 0.01 s (Fig. 4.B), and enabled replicating the full experimental net discharge rate range of 9.45 to 14.69 Hz (Fig. 4.C). The firing characteristics of identified in vivo MNs and their in silico counterparts (before and after optimization) from all trials are depicted per subject in Fig. 5. As shown by the boxplots in Fig. 6, the absolute errors in recruitment time and mean firing frequency (during plateau) were significantly reduced across all conditions after parameter optimization: the absolute recruitment time error median decreased from 365.7 ms to 0.48 ms, and the maximum outlier value from 2.69 s to 0.017 s. For absolute discharge rate error, optimization decreased the median value Example of in silico MNs spike trains and estimated net discharge rate of the MN pool, before (i.e., generic) and after (i.e., optimized) parameter optimization (subject 2). A) in vivo (blue) and generic in silico (red) MN spike trains. B) Optimized in silico MN spike trains (red) (matching in vivo MNs (blue). C) Estimated net discharge rate of MNs (i.e., low pass filtered cumulative spike train, normalized by average discharge rate of all MNs) for each % MVC amplitude using spike trains from in vivo (blue), generic in silico (orange) and optimized in silico (red) MNs.Error bars depict standard deviation. from 5.8858 Hz to of 0.208 Hz, and the maximum outlier from 13.86 HZ to 5.75 Hz.
Test 4: Our proposed algorithm for extending the number of in silico MNs to complete a pool size of 200 MNs did not alter the original parameter's PDFs (Fig. 7). This held true regardless of the number of in vivo MNs identified for each subject (e.g., for subjects one and four, only 18% and 30% of the MNs in the pool were identified, respectively, whereas 64% was identified for subject three). In all cases, the PDFs remained largely unaltered. There were no considerable differences in the median or STD of the subject-specific PDFs of the extended in silico MN pools (Table II).
Test 5: Fig. 8 shows an example of a subject-specific in silico MN pool (from subject four) driven by the CSI corresponding to each target %MVC amplitude. As depicted by the raster plots (Fig. 8.A) showing the firing activity of the in silico MNs, the percentage of pool recruitment increased logarithmically as torque increased, going from ∼60% recruitment at 10 % MVC, to ∼100% recruitment at 50%MVC. For all %MVC amplitudes (Fig. 8.B), the neural drive produced by the in silico MN pool closely resembled the corresponding force profile, achieving an R 2 = 0.86 ± 0.04 with p-values < 0.005, and CoV = 68.12 ± 7%. In contrast, the neural drive estimated from in vivo MNs (Fig. 8.C) resulted in R 2 = 0.723 ± 0.08 and CoV = 75.47 ± 11.78%. Fig. 9 shows R 2 and CoV of in vivo MNs, optimized in silico MNs, and extended in silico pool across all subjects and conditions. The non-parametric Friedman test showed no significant difference in either R 2 or CoV between in vivo and in silico MNs, indicating that in silico MN pools generated through this framework closely mimic the firing characteristics of their in vivo counterparts. Between in vivo MNs and extended in silico pool, however, non-parametric Friedman test showed statistically significant improvement in R 2 , which increased from 0.7997 ± 0.08 to 0.8716 ± 0.06. This represents a relative improvement of 9.66 ± 7.41% using the in silico MN pool, in comparison to only in vivo MNs. Although not statistically significant, the CoV corresponding to the extended in silico pool decreased relatively to the in vivo MNs by 5.16 ± 4.49 %.

IV. DISCUSSION
This work combined HD-EMG decomposition, biophysical neuronal modelling and metaheuristic optimization to characterize the subject-specific physiological properties of complete pools of MNs innervating TA in healthy individuals. For this, we established in silico MN models driven by in vivoderived CSI, and optimized their parameters to reproduce the firing behavior of in vivo MNs. Subsequently, we used PDFs from the identified parameters as blueprints to create complete subject-specific in silico MN pools capable of estimating MN pool activity and muscle activation profiles for TA throughout multiple levels of isometric activation.
The excitability of a single MN is described by its currentfrequency slope [38]. Analogously, here we proposed I F as the slope between CSI and the net firing activity of the MN pool. Hence, I F can be interpreted as a subject-specific parameter reflecting spinal excitability. We found relatively similar values of I F for all subjects (table I). This may be due to the inclusion of only healthy subjects from a population of similar ages and sizes. Notably, the identified values of I F matched with experimental current-frequency slopes reported for single MNs (1.7556 ± 0.4953 Hz A ) in the primary firing range (<40 Hz) [38]. Additional studies are required to assess whether I F may reflect the averaged current-frequency slopes of single MNs within the pool.
As this framework relies on the availability of in vivo MNs representing the entire pool, the experimental protocol is key to ensure the activation of different types of motor units. Here, we used trapezoidal ramps at a constant rate of force development where target forces were reached in 0.5 < t < 2.5 s. Evidence of decreased recruitment thresholds in TA's MNs at this rates [35] suggests that MNs from all thresholds of recruitment could be activated for the levels of forces included in this study [39]. However, HD-EMG decomposition does not guarantee the identification of every MN, particularly from deeper motor units or in the upper/lower end of recruitment order. Due to the importance of the latter for estimating I F, not identified MNs of the earliest/latest recruitment order in the muscle may lead to over/under estimated I Fs, respectively. Additional studies are necessary to assess the impact that this may have over the resulting parameter's PDFs.
In agreement with Henneman's size principle [40], we found that optimization of D S substantially minimized recruitment error (5) (Fig. 2), evidencing the key role of MN soma size over recruitment order. Spike-match (6), minimized through optimization of β Q and α Q , yielded relatively larger cost values in comparison. This may be due to the short time window used for determining coincident spikes (δ = 2ms), as well as the exclusion of other ionic channels. Spikematch could be further enhanced by optimizing additional ionic channel properties, thought it would increase computational load and parameter space. In terms of repeatability, D S showed neglectable variability throughout repeated optimizations (Fig. 2). This is not the same for β Q and α Q , which can be explained by the intrinsic randomness behind temporal MN firing patterns and metaheuristic optimization. However, although both parameters showed variability across repeated optimizations, each oscillated around a mean value seemingly related to D S . As shown in Fig. 3, we found that value and spread of β Q and α Q increased with D S , suggesting proportionality be between voltage-dependent rates and MN soma size.  Calculated R 2 and CoV for in vivo MNs (blue), optimized in silico MNs (red), and extended in silico pool including PDF-derived MNs (green) across all subjects and conditions. Statistical difference between in vivo MNs and in silico pool, as determined by non-parametric Friedman test, indicated by a star.
In general, we observed smaller β Q values associated to higher discharge rates and, thus, earlier recruited (i.e., smaller size) MNs. This is consistent with the size-based parameters proposed by Kohn and Abdala for different types of MNs [14]. This framework assumed the same absolute ranges of MN parameters for every subject (i.e., inter-subject differences are given by the PDFs). Here, we found disctinctive PDFs for every subject (test 2, Fig. 3), with largest inter-subject variability observed in D S . For subject 1, we observed a distinctive D S distribution, which could be explained by the limited number of decoded MNs in comparison with other subjects. However, since our framework did not include MN-tracking to control for repeated units (i.e., same MN identified in multiple contractions), we cannot discard that PDFs derive from larger MN counts may be biassed by repeated units, particularly in the lower-end of recruitment order. Future work should implement MN-tracking and assess the impact that number of decomposed MNs has on the estimation of PDFs represantive of the entire MN pool. Furthermore, future work should test wheter fixed ranges of MN parameters enable characterizing pool properties in elderly and imparied subjects, as reported adaptations in motor unit sizes and contraction speeds [41] have been suggested to be a consequence of age-related losses of larger size MNs [1], while altered MN soma sizes have been related to neuronal lesion and motor disorders [42]. Traditionally, motor unit types (i.e., slow, fatigue resistant and fast fatigable) have been associated to MN soma size. Studies of human autopsies [37] reported that approximately 70 % of the motor units innervating TA are associated to small size MNs [18]. Our results showed that, for every subject, the dominant clusters C 1 and C 2 comprised together ≈ 70 % of the MN pool. Future work will systematically assess whether this cluster-based classification of MN parameters could be used for in vivo identification of distinct motor unit types.
HD-EMG decomposition revealed that the firing frequency of in vivo MNs increased proportionally to % MVC (Fig. 4.C). Although we did not control for coactivation of antagonist muscles (i.e., MN activity may be overestimated, particularly at higher % MVCs), this is consistent with previous findings [39] suggesting that rate coding has a dominant role in the modulation of force during high speed isometric contractions. After parameter optimization, in silico MNs driven by in vivo CSI reproduced the same degree of frequency modulation (Fig. 4.C) and matched both recruitment threshold and discharge rate (Fig. 5) of in vivo MNs for all subjects and conditions (Fig. 6). Although these results do not validate the underlying in silico MN properties, especially considering the limited set of isometric contractions included in the experiment, this is the first work able to create neural-data driven models that capture subject-specific dynamics of in vivo MNs. Moreover, parameter ranges for optimization were derived from the closest approximations available to in vivo human properties [26]. Given the technical challenge of measuring human MN properties in vivo, matching recruitment threshold and discharge rates is currently the only viable alternative for in vivo validation. Future work will aim at implementing a neuromuscular model for validation at the level of generated torque.
The resulting parameter's PDFs enabled extending the number of in silico MNs to complete a pool of 200 MNs (Fig. 7) while retaining the same PDFs as initially identified. Since identified PDFs of β Q were skewed to the lower end of the distribution for every subject (Fig. 7), given the inverse relationship we found between β Q and discharge rate, we implemented a quality-control algorithm to prevent PDF-derived in silico MNs from firing at discharge rates significatively larger than in vivo MNs. Here, we defined a tolerance of 3Hz. However, future work should systematically assess whether this threshold and rejection criteria are enough to compensate for not identified in vivo MNs. In comparison to the neural drive estimated from in vivo MNs (Fig. 8. D), this approach of driving in silico MN pools with trial-specific in vivo CSIs resulted in better estimations of force profiles as indicated by the higher R 2 achieved by the in silico pools ( Fig. 8.C). This may be due to the availability of neural data from entire MN pools, which is not accessible via HD-EMG alone, as larger motor units (i.e. producing action potentials of higher amplitude) overshadow smaller units [43] and prevent identifying earlier recruited MNs. This is particularly noticeable at 30%, 40% and 50% MVC, where the in vivo neural drive shows large fluctuations and mismatches during the ascending/descending portions of the ramps (Fig.8.C). In contrast, additional spike trains from the in silico MN pool resulted in smoother neural drive estimates that follow more closely the torque profile ( Fig. 8.B).Moreover, we found that R 2 of the in silico pool increased with % MVC. This could be explained by the larger number of recruited MNs at higher % MVCs (7). With more MNs contributing to the neural drive, the transmission of the CSI may be improved and the CoV reduced, thereby enhancing force steadiness [44]. Regardless, statistical significant improvement in R 2 was found for all subjects and conditions using in silico MN pools (Fig. 9). In the future, this framework could be used to complement HD-EMG decomposition for deriving mechanically consistent MN spike trains (e.g., not exhibiting pauses, noise-like spikes or merging inconsistencies). However, to validate the percentage of recruited MNs from the in silico pool as a function of force magnitude, future work should look into the implementation of MN-specific twitch models.
As MN excitability is thought to be modulated during dynamic contractions and may vary with rate of force development [35], future work should also evaluate whether this approach for driving in silico pools of MNs generalize to different motor conditions. Given the challenge of predicting broader repertoires of motor tasks, future work may also look into the addition of dendritic neuromodulation to our MN models, as persistent-inward currents have been reported to play a substantial role modulating MN behavior [45].
Lastly, we emphasize that all MN pool characteristics here discussed were derived from TA, the main contributing muscle to ankle dorsi-flexion. Applying this method to other muscles may result in distinctive motor unit sizes and distributions [32]. However, although additional studies are required, we believe this same framework can be used for any muscle where the underlying neural data is available, as ranges for parameter optimization included measurements from several limb muscles [13], [26], and large inter-muscle similarities in contractile and discharge characteristics have been previously reported [32]. Future work will focus on applying this framework to different muscles and addressing the possibility of creating comprehensive subject-specific neuro-musculoskeletal models for simulating the neuro-mechanical interaction of multiple muscles.

V. CONCLUSION
This study demonstrated the ability of our neural-data driven framework for creating subject-specific in silico MN pools that enable predicting in vivo MN firing characteristics and muscle activation profiles throughout isometric force-tracking tasks at varying levels of amplitude. Thus, bridging the gap between current neuromechanical models and clinical practice.