A novel motoneuron-driven computational muscle model with motor unit resolution for subject-specific simulations of human voluntary muscle contraction Abbreviated title: A novel motoneuron-driven muscle model

Department of Civil and Environmental Engineering, Imperial College London, SW7 2AZ, UK Department of Bioengineering, Imperial College London, SW7 2AZ, UK Graduate School of Biomedical Engineering, University of New South Wales, Sydney 2052, Australia Emails: arnault.caillet17@imperial.ac.uk; d.farina@imperial.ac.uk; l.modenese@unsw.edu.au * These authors have equal contributions and share the senior authorship


23
During voluntary skeletal muscle contractions, the pool of alpha-motoneurons (MNs) and the muscle 24 fibres they innervate, constituting the muscle's population of motor units (MUs), transform in a series 25 of chemical-mechanical events the neural message from the spinal cord into molecular forces (Huxley 26 & Niedergerke, 1954). These forces sum across the multiscale architecture of the muscle to build the 27 individual MU forces and the whole muscle force. Extensive experimental investigations in animals in 28 vivo have advanced our understanding of the neuromechanical dynamics responsible for voluntary 29 muscle contraction. For example, the force generated by skeletal muscles is modulated by the 30 discharge frequencies of the innervating MNs (Heckman & Enoka, 2012), and by the number of 31 discharging MUs, which are sequentially recruited in the increasing order of their size and the force 32 they can generate (Henneman, 1981 To support experimental investigations into understanding and simulating the interplay between 45 human motor control and force generation, mathematical models of skeletal muscles were developed 46 to mimic the dynamics of muscle contraction. Hill-type neuromuscular models are popular solutions 47 TA muscle force . was predicted (G) from the experimental (t) trains in a blind approach (blue dotted lines), after location of the identified MUs into the MU pool (green dashed lines),

140
and from the completely reconstructed population of spike trains (solid red lines). The neuromuscular model was scaled with subject-specific muscle-tendon properties derived from a subject-141 specific MSK model (I), obtained from segmented MRI scans (H). The predicted TA force ,1 , ,2 , or , depending on the type of neural control, was finally validated (G) against the 142 experimental TA force (J). was estimated from the recorded ankle torque and the estimated co-activity of the TA's agonist and antagonist muscles crossing the ankle joint.
2. Experimental data and subject-specific neuromusculoskeletal quantities 144 2.1. Contraction tasks, force and EMG recordings, and subject-specific spiking activity 145 The experimental approach represented in Figure 1A was described elsewhere (Caillet et al., 2023) and 146 is briefly summarized here. One healthy male participant (age: 26 yr; height: 168 cm; body weight: 52 147 kg) volunteered to participate in the experimental session of the study. The participant sat on a 148 massage table with the hips flexed at 30°, 0° being the hip neutral position, and his knees fully extended 149 and strapped to the massage table to avoid any knee and hip motion and co-contraction of the thigh 150 muscles during the ankle dorsiflexion. The ankle torque signals ( ) were recorded with a load cell 151 (CCT Transducer s.a.s, Turin, Italy) connected in-series to the pedal of a commercial dynamometer (OT 152 Bioelettronica, Turin, Italy), to which the participant's foot was fixed at 30° in the plantarflexion 153 direction, 0° being the foot perpendicular to the shank. The effects of gravity, including the rig weight, 154 were cancelled by removing the initial measured offset during the experiments. After registering the 155 participant's maximum voluntary contraction (MVC) in dorsiflexion, the participant performed two 156 trapezoidal isometric contractions at 30% and 50% MVC with 120 s of rest in between, consisting of 157 linear ramps up and down performed at 5%/s and a plateau maintained for 20 s and 15 s at 30% and 158 50% MVC, respectively. 159 To derive subject-specific motoneuronal controls to the neuromuscular model, HDEMG signals were 160 recorded over the TA muscle during a first experimental session with a total of 256 EMG electrodes 161 ( Figure 1B), distributed between four rectangular grids of 64 electrodes, with 4 mm interelectrode 162 distance covering 36 cm 2 of the muscle surface (10 cm x 3.6 cm; OT Bioelettronica). HDEMG and force 163 signals were recorded using the same acquisition system (EMG-Quattrocento; OT Bioelettronica) with 164 a sampling frequency of 2,048 Hz. As previously described (Caillet et al., 2023), the 256-electrode grid 165 was downsampled by successively discarding rows and columns of electrodes and artificially 166 generating two new grids of lower electrode density ( Figure 1B). The two grids covered the same 167 muscle area with interelectrode distance of 8 mm and 12 mm, involving 64 and 36 electrodes, 168 respectively. After visual inspection and band-pass filtering, the HDEMG signals recorded with the 169 three grids were decomposed into motor unit spike trains using convolutive blind-source separation, 170 as previously described . After the automatic identification of the motor units, 171 duplicates were removed, and all the motor unit spike trains were visually checked for false positives 172 and false negatives (Del Vecchio et al., 2020). Only the motor units which exhibited a pulse-to-noise 173 ratio > 28 dB were retained for further analysis. The final spike trains ( ) identified 174 experimentally ( Figure 1C) were stored as binary vectors of zeros and ones identifying the sample times 175 when a discharge occurred. The ankle torques ℎ at which the MUs were recruited, also called 176 torque recruitment thresholds in the following, were calculated as the average of the recorded values 177 over a 10 ms range centred around the MUs' first identified discharge time. In the following, the 178 MUs identified experimentally were ranked in the ascending order of measured recruitment torques 179 ℎ with the index ∈ ⟦1; ⟧. According to the Henneman's size principle (Henneman, 1981;Caillet 180 et al., 2022b), these MUs were therefore also ranked in the same order of current recruitment 181 threshold ℎ , maximum isometric forces 0 , and innervation ratios according to EQ. 1. 182 ∀ , ∈ ⟦1; ⟧, < ⟺ ℎ ( ) < ℎ ( ) ⟺ ℎ ( ) < ℎ ( ) To scale the neuromuscular model to subject-specific MSK quantities, a 3D T1-weighted VIBE 194 (volumetric interpolated breath-hold examination) sequence was used to acquire high resolution 195 images (0.45 x 0.45 mm pixel size, 1 mm slice thickness and increment) from the knee joint to mid-foot 196 of the participant's dominant leg. The sequence consisted in three blocks with 50 mm overlap, during 197 which the participant was asked not to move to facilitate a successful merging of the adjacent blocks. 198 The volumetric shapes of the tibia, fibula, and foot bones, and of the TA muscle belly and tendon were 199 carefully identified by manual segmentation ( Figure 1H) using ITK-SNAP (Yushkevich et al., 2006) and 200 an anatomical atlas as support (Möller & Reif, 2007). Slices of the agonist EDL and EHL and antagonist 201 SOL, GM, and GL muscle bellies were segmented at regular intervals along the tibial length so that their 202 centroidal lines of action could be outlined. All the main tendons crossing the ankle joint could be 203 clearly identified and were fully segmented. An expert radiologist confirmed the accuracy of the 204 segmentation. The segmented shapes were manually adjusted on the transversal plane to avoid 205 inconsistency of the muscle geometry in correspondence of the MRI sequence blocks using Netfabb 206 (https://www.autodesk.com/products/netfabb/) and Meshlab (Cignoni et al., 2008). The TA muscle  207  volume  computed from the segmentation was found consistent with the relationships between  208   anthropometry and muscle volumes proposed by Handsfield et al. (2014) from a cohort of segmented  209   MRI scans. Therefore, the volume  of the EDL, EHL, SOL, GM, GL muscle bellies was estimated with  210   the volume-height-mass  To reproduce the experimental conditions, the plantarflexion angle of the final subject-specific MSK 222 model was set to 30°. At this angle, the subject-specific muscle-tendon length of the six muscles 223 and their moment arm with respect to the ankle joint were computed using OpenSim. The subject-224 muscle-specific optimal length 0 and the tendon slack length of the six muscles were estimated 225 from the muscle-specific values 0, and , proposed in a generic published model (Rajagopal et 226 al., 2016), which were scaled with EQ. 2 by maintaining their ratio with respect to the entire 227 musculotendon lengths and in the neutral ankle position. 228 Considering a specific tetanic tension of = 60 / 2 (Rajagopal et al., 2016), the subject-specific 229 maximum isometric force 0 of the six muscles was obtained from the known muscle volumes 230 with EQ. 3. 231 Combining these subject-muscle-specific properties with the bEMG signals recorded from the co-232 contracting flexor muscles during the second experimental session ( Figure 1A), the experimental force 233 ( ) developed by the TA during the first experimental session was inferred ( Figure 1J) from the 234 measured experimental ankle torque ( ) with EQ. 4. In EQ. 4, is the moment arm the TA tendon 235 makes with the ankle joint, and Δ ( ) relates the level of measured ankle torque ( ) to the algebraic 236 amount of ankle torque Δ taken by the group of agonist EHL and EDL and antagonist SOL, GM, and 237 GL muscles. The continuous Δ ( ) relationship was obtained by trendline fitting the ( ; Δ )  The MUs identified in the experimental trials were ranked with the index ( ∈ ⟦1; ⟧, with < 251 ) in the ascending order of measured ℎ . Using the ℎ ( ) distribution in EQ. 5, these experimental 252 MUs were mapped into the complete ℎ -ranked pool of 400 MUs, where they were assigned a new 253 ) , j ∈ ⟦1; ⟧ EQ. 7 The normalized 0 ̅̅̅̅̅̅ ( ) distribution in EQ. 8 was then estimated by scaling the coefficients in EQ. 7 to 264 obtain a 11-fold range for the 0 ̅̅̅̅̅̅ values in the TA MU pool, to account for the different twitch-tetanus 265 ratios obtained between slow and fast fibres in the rodent literature (Celichowski & Grottel, 1993 Figure 2D. 274 Although it is debated (Caillet et al., 2022b) whether the MU type, currently determined by the MU 275 twitch force , is also correlated or not to the MU size, i.e., to the MU Force recruitment threshold 276 in humans, it was here assumed that the small low-threshold MUs were of the 'slow' type, and the 277 large high-threshold MUs of the 'fast' type. Therefore, from the cumulative summation of the MU , 278 the 359 lowest-threshold MUs that have the lowest innervation ratios were considered to account for 279 72% of and to be of the 'slow' type in the following. 280  Figure 1F) and its passive force dynamics were neglected, {2} all the MU actuators were assigned the 304 same optimal length, set as the muscle's 0 , and the same normalized constant length, set as ̅ = 1.16 305 according the subject-specific MSK measurements, {3} the Force-Velocity properties of the MUs were 306 therefore neglected, and {4} the passive force of the CE was found negligible, so that the passive elastic 307 element (PEE, in grey in Figure 1F) was neglected. With these assumptions, and after normalizing the 308 state variables to the 0 and 0 parameters and reporting them with a bar, as described in Appendix 309 A.2, the whole muscle dynamics in isometric contractions were described with EQ. 10. In EQ. 10, at 310 time and for a population of modelled MUs, is the whole muscle force in Newtons, ̅̅̅̅̅̅ is the 311 normalized force developed by the ℎ MU ( ∈ ⟦1; ⟧), and ,0 ̅̅̅̅̅̅ is the MU-specific normalized 312 maximum isometric force that the k th MU can develop. In EQ. 10, the MU active states ( ̅ ) and Force-313 Length (FL) scaling factors ( , ̅ ) were derived from the input spike trains ( ) and the common 314 MU length ̅ with six cascading mathematical models of the MU's excitation, activation, and 315 contraction dynamics summarized in Figure 3. Finally, because the MUs build force upon the discrete 316 activity of their = • ,0 ̅̅̅̅̅̅ fibres that receive the same MN action potential with a random time 317 delay , the calculation of the normalized MU forces ̅̅̅̅̅̅ in EQ. 10 was updated to ̅̅̅̅̅̅ to account 318 for this non-simultaneous twitch activity, which acts as a low-pass filter.

EQ. 10
For improved readability, the subscript was discarded in the following of the Methods section to 320 refer to the dynamics of a representative MU of the pool. 321  to-MU multiscale assumption is here acceptable as muscles modelled as a scaled sarcomere for their 343 normalized isometric FL properties can provide accurate Hill-type model predictions for mammalian 344 whole muscles (Winters et al., 2011). Other available solutions were considered less appropriate, such 345 as torque-length and torque-angle measurements in humans which are affected by muscle-tendon 346 dynamics and muscle co-contraction, or FL measurements from fibre bundles in cat, rabbit and rat 347 specimens, that have different optimal sarcomere lengths than humans (Burkholder & Lieber, 2001). 348 Rather than a piecewise linear relationship consistent with the sliding filament theory at the sarcomere 349 scale (Gordon et al., 1966), a gaussian mathematical description in EQ. 11 was fitted to the FL 350 experimental data (red curve and black dots in Figure 4A 2022b) were fitted with EQ. 12, yielding = 90 and = 1.4 (Table 1). is here shorter ( Figure  364 4B) than the amphibian = 2 (Ebashi & Endo, 1968) used in a previous neuromuscular model 365 (Hatze, 1977). 366 Dynamics of fibre AP elicitation 368 The train of the fibre APs ( ( ) in Figure 3) was modelled in EQ. 13 as the 2 nd -order response (Hatze, 369 1977) to the train ( ) of MN APs. In EQ. 13, when the fibre receives an isolated MN discharge, 1 370 determines the peak-to-peak amplitude of the fibre AP, while 2 and 3 respectively determine its 371 time-to-peak duration and exponential decay time constant . The fibre APs simulated with the 372 tuned coefficients in Table 1  previously proposed for amphibians at low temperature (Hatze, 1977). 377 The MU excitation state ( ) was modelled to control the level ( ) of calcium concentration [ 2+ ] 381 in the MU sarcoplasm ( Figure 3) with the 2 nd -order ODE in EQ. 14, although alternative reviewed 382 approaches were previously proposed (Caillet et al., 2022c). In EQ. 14, when the fibre receives an 383 isolated fibre AP, 1 determines the peak-to-peak amplitude of the [ 2+ ] twitch, and 2 and 3 384 its time-to-peak duration and decay time constant . The equation proposed in Hatze (1977) was 385 modified to introduce the length-dependent scaling factors 1 ( ̅ ) and 2 ( ̅ ), defined in EQ. 15 and EQ. 386 16, that respectively make and nonlinearly dependent to the normalized MU length ̅ . The 387 procedure taken to tune the coefficients and derive the ( ̅ ) scaling factors to match available 388 experimental data for amphibian and mammal slow and fast fibres in vitro is reported in Appendix A.3. 389 As reported in Figure 4D for fast fibres, the calcium transients simulated with the tuned coefficients 390 (Table 1) & Hollingworth, 2003) were scaled to 35°C using the same scaling factors. As shown in 400 Figure 4C, the simulated calcium transients at 35°C for slow fibres under 67Hz stimulation compared 401 well with the scaled experimental data with similar , and of 19 and 44ms and an average 402 underestimation of of 20%. The calcium twitches simulated at optimal fibre length with EQ. 14 and 403 the parameter values in Table 1 have ten times shorter half-widths and 2.3 times lower amplitudes 404 ( Figure 4E, F) than the twitches obtained with Hatze's model (Hatze, 1977), which was calibrated on 405 frog data at 9°C. 406 407 Dynamics of Calcium-Troponin concentration 408 The sarcomeric calcium-troponin (CaTn) binding-unbinding process was described with EQ.

Dynamics of MU activation 424
The MU active state can be defined as the normalized instantaneous ratio of formed cross-bridges in 425 the force-generating state in a population of myosin heads from a representative sarcomere (Caillet et  426 al., 2022c). Although semi-phenomenological models of multi-state cross-bridge attachments (Zot & 427 Hasbun, 2016) are available, a phenomenological model (Ding et al., 2000) was used to infer the MU 428 active state from the CaTn concentration ( Figure 3) to reduce the computational cost, the number 429 of parameters to tune, and the overall model complexity. In EQ. 18, 1 , 2 and 3 respectively control 430 the amplitude of the activation twitch, its time-to-peak and its time of half-relaxation 0.5 , and 431 were tuned to match, lacking more suitable experimental measurements in the literature, the Because of the MU type-specificity of the [ 2+ ] and CaTn dynamics modelled previously, a unique 436 set of values was obtained for the tuned 1 , 2 and 3 coefficients (Table 1)  of 30% at 15Hz with the extensor digitorum brevis (EDB) muscle data (Macefield et al., 1996) and less 458 than 10% at all physiological frequencies above 5Hz with the data obtained from thenar muscles 459 ( EQ. 18, that describe the MU-specific excitation and activation dynamics, were first numerically solved 466 at each time step, yielding a vector of MU active state time-histories ( ). Figure 4I, J displays the 467 time-histories of the five state variables that describe the NE's activity at 1 and 30 Hz, respectively. An 468 estimation of the MU normalized forces and of the whole muscle force was then obtained with EQ. 10. 469    (Hatze, 1977).

484
(C-F) Time-course of the free calcium concentration in for slow (C, E) and fast (D, F) MUs. The calcium 485 transients simulated after tuning the coefficients in EQ. 14 to match experimental rodent data (Table 1)   were assigned with EQ. 19 representative maximum isometric forces ,0 ̅̅̅̅̅̅ , necessary parameter in EQ. 513 10, to account for the force-generating capacities of their neighbouring MUs, in the sense of their 514 recruitment thresholds, that were not identified experimentally. The muscle force ,1 ( ) (blue 515 dotted line in Figure 1G) was first predicted when the identified MUs were assumed to be 516 homogeneously spread across the ℎ -ranked MU pool, in which case the 0 ( ) distribution in EQ. 8 517 was evenly split in domains along the x-axis ( Figure 2B). A second TA force ,2 ( ) (green dashed 518 line in Figure 1G) was predicted after assigning to the MUs representative maximum forces with 519 EQ. 19 after a physiological mapping of these MUs to the complete ℎ -ranked MU pool using EQ. 6. 520 For example, when = 30 MUs are identified at 30% MVC, where =300 MUs are supposed to be 521 recruited (EQ. 5), the 9 th , 10 th , and 11 th identified MUs in the ℎ -ranked experimental sample are 522 mapped into the complete ℎ -ranked pool to the 90 th , 100 th , and 110 th MUs with the first method 523 The TA force ( ) (red solid line in Figure 1G) was predicted from a comprehensive cohort of = 528 400 simulated spike trains, that described the discharge activity of a completely reconstructed MU 529 population (red solid line in Figure 1E Figure 1E, 552 respectively) derived from the three grids of increasing electrode density ( Figure 1B). The subject-specific MSK model ( Figure 1I) yielded the subject-muscle-specific properties in Table 2 for 575 the TA muscle and for its agonist EDL and EHL, and antagonist muscles SOL, GM, and GL in ankle 576 dorsiflexion. By reapplying with EQ. 4 the exponential Δ ( ( )) relationship derived from the bEMG 577 and torque recordings obtained during the second experimental session to the torque traces ( ) 578 recorded during the first experimental session, the experimental TA force ( ) (solid red traces in 579

590
The black dashed traces display the TA force calculated without considering the co-contraction of surrounding 591 muscles (as if the TA alone produced the entire recorded ankle torque ).  Figure 1C) 604 were respectively identified (Caillet et al., 2023) with the EMG grids of 4, 8, and 12 mm interelectrode 605 distance ( Figure 1B). The dataset of three identified MUs was disregarded in the following results 606 because it provided a highly inaccurate description of the discharge activity of the real MU pool. The 607 histogram distribution of those identified MUs across the MU pool according to their recruitment 608 threshold ℎ in %MVC is displayed in Figure 6. Importantly, the denser the EMG grid and the lower 609 the generated force, the more homogeneous the distribution of the identified MUs across the MU pool 610 (Caillet et al., 2023). For example, 12% to 20% of the MUs decoded with the densest grid at 30% MVC 611 were systematically identified in each of the 5% ranges sampling the recruitment range (histogram in 612 Figure 6C). The normalized effective neural drive to muscle ̅ ( ) ( Figure 1E) was computed with these 613 cohorts of spike trains and was validated against the experimental normalized TA force ̅̅̅̅̅ ( ) (black 614 solid traces in Figure 6) with the three validation metrics reported in Table 3. 615

592
The normalized effective neural drive ̅ ( ) estimated from the experimental samples of MUs (blue 616 dotted lines in Figure 6) was always better approximated when the electrode density increased, 617 especially in the regions of low generated forces. For example, at 30% MVC, it was obtained Δ 1 = 618 0.07 , = 8.0 %, and 2 = 0.97 from the 81 discharging MUs identified with 4 mm 619 interelectrode distance, while the 14 MUs identified with 12 mm interelectrode distance returned 620 Δ 1 = 0.14 , = 17.5 %, and 2 = 0.89 (Table 3). The effective neural drive ̅ ( ) estimated 621 from the experimental MU samples was always better approximated at 30% MVC than at 50% MVC 622 (blue dotted traces in upper row versus lower row in Figure 6), with two to three times lower nRMSE 623 values, and strongly lower Δ 1 and higher 2 values, respectively (Table 3). 624 Complete populations of = 400 MNs ( Figure 1D) were also reconstructed from the discharge activity 625 of the identified MNs by estimating the discharge activity of the MUs not identified experimentally 626 (Caillet et al., 2022a). In this case, the ̅ ( ) estimated at 30% MVC (solid red lines, upper row, in Figure  627 6) was accurately predicted and was unrelated to the quality of the input experimental samples, with 628 < 7%, and 2 = 0.98 for the three EMG grid configurations (Table 3). At 50% MVC, the ̅ ( ) 629 derived from the MNs (solid red lines, bottom row, in Figure 6) was better estimated when the MU 630 pool was reconstructed from denser grids of electrodes, with lower nRMSE and higher 2 values (Table  631 3). As for the experimental MU samples, ̅ ( ) estimated from the = 400 MUs was better 632 approximated at 30% MVC than at 50% MVC (solid red traces in upper row versus lower row in Figure  633 6). At 50% MVC, the reconstructed populations of discharging MUs returned inaccurate estimations of 634 the onset and ascending ramp of force, with Δ 1 > 1 , > 12%, and 2 ≤ 0.95. Overall, the 635 normalized effective neural drive ̅ ( ) was systematically better approximated by the reconstructed 636 pool of = 400 MUs than by the experimental samples of identified MUs (solid red versus blue 637 dotted traces in Figure 6) with lower nRMSE and higher 2 values (Table 3). 638 Therefore, ̅ ( ) was best approximated by samples of MU spike trains of identified MUs that spanned 639 across the whole recruitment range and followed a homogeneous distribution across the MU pool. 640 Obtaining such representative description of the discharge behaviour of the real MU pool was possible 641 with EMG signals recorded at 30% MVC with high electrode density (4 mm interelectrode distance) or, 642 to some extent, with the computational reconstruction of the complete MU pool, both of which aimed 643 to correct for the systematic bias of HDEMG decomposition towards identifying the higher-threshold 644 MUs (histograms in Figure 6   Depending on the modelling approach, the modelled MUs were assigned maximum MU forces ,0 683 (black dots in Figure 7B-D) by either blindly assuming the MUs to be evenly distributed across the 684 MU pool (Figure 7B), mapping the MUs to the real MU pool with EQ. 6 ( Figure 7C), or directly 685 applying the 0 ( ) distribution in EQ. 8 to the complete MU pool ( Figure 7D). Depending on the 686 approach, the individual MUs produced maximum forces over the contraction up to 30% MVC in 687 the range 1.8-6.6 N, 0.4-16.5 N, and 0.5-1.5N (blue crosses in Figure 7B, green triangles in Figure 7C, 688 and red dots in Figure 7D

700
To do so, in B, the identified MUs were assumed to be homogeneously distributed across the ℎ -ranked MU 701 pool, while they were accurately mapped to the real pool with EQ. 6 in C. In B, the non-continuous distribution of 702 the ,0 values is explained by the fact that is not an integer value in EQ. 19.

704
The simulated MU forces were processed and linearly summed with EQ. 10 to yield the three predicted 705 whole muscle forces in Figure 8A Table 4. 709 When high numbers of MUs were homogeneously identified over the full range or recruitment ( = 710 81 MUs, 30% MVC, 4 mm interelectrode distance), i.e., when ̅ ( ) was accurately estimated from the 711 experimental samples of spike trains ( Figure 6C), was predicted by ,1 , ,2 , and with the 712 same level of accuracy ( Figure 8C) for the eight validation metrics in Table 4 with 2 ≥ 0.98, all 713 values below 10%, and ME in 36-66 %, irrespective from the three choices of 0, assignment (black 714 dots in Figure 7B-D). 715 In all other cases, ( ) was the most accurately predicted by both ( ) when the complete MU 716 pool was modelled (solid red traces in Figure 8A (Table 4). 742 As expected from the neural drive estimation in Figure 6, for all the modelling approaches and EMG 743 grid configurations, was more accurately predicted at 30% MVC than at 50% MVC (upper row 744 versus lower row in Figure 6), with two to three times lower nRMSE values, and strongly lower Δ 1 and 745 Δ 1 F and higher 2 values, respectively (Table 4)

751
The experimental TA force ( ) was estimated from the recorded ankle torque with EQ. 4. This study reports a novel approach, summarized in Figure Figure 1E). Finally, the whole muscle force predicted by the MN-driven neuromuscular model 783 was validated against an experimental muscle force ( Figure 1G, Figure 8), that was estimated from 784 measured joint torque and bEMG signals recorded from agonist and antagonist muscles ( Figure 1J, 785 type, some of which were never included in neuromuscular models. This description of the muscle as 830 a cohort of MUs also provided a convenient structure for deriving physiological assumptions and 831 simplifications to the rheological structure of the model (Appendix A.2) and opens avenues to 832 accounting for the impact of the muscle's internal architecture on its force-generating capacities, 833 especially if the neural control is mapped to an appropriately discretised anatomical model (Modenese 834 & Kohout, 2020). This approach has only rarely been modelled in Hill-type muscle modelling (Hatze, 835 1980;Woittiez et al., 1984;Kaufman et al., 1991). 836

753
This MU-scale approach was ideal for developing advanced models of the MU's excitation, activation, 837 and contraction dynamics (Figure 4). Experimental data from the literature was used in this study to 838 derive adequate mathematical equations (EQ. 11 to EQ. 18), tune physiological parameters (Table 1), 839 and validate the modelled dynamics (Figure 4). Importantly, the calibration of the parameters in EQ. 840 11 to EQ. 18 and the validation of those mathematical descriptions were performed using different 841 sets of experimental data from the literature that were measured at different discharge rates (see 842 Appendices A.3 and A.4 for details). Furthermore, we limited the amount of multiscale and inter-843 species scaling inaccuracy, which is an important limitation of single-actuator Hill-type models (Caillet 844 et al., 2022c), by considering source experimental data, in turn and decreasing order of preference, 845 from studies on individual human MUs or bundles of fibres, human fibres and sarcomeres, cat fibres, 846 rodent fibres, rodent muscles, and amphibian fibres or muscles. In doing so, the simulated MN APs, 847 fibre APs, and free calcium twitches were more physiological and consistent with mammalian dynamics 848 than previous approaches based on experimental amphibian data at low temperature ( Figure 4B The accuracy of the whole muscle force predictions in Figure 8 was mainly related to the accuracy of 863 the experimental neural control in approximating the neural drive to muscle (Figure 6). It is currently 864 impossible to identify the discharge activity of the complete MU pool of a muscle with surface EMG 865 because of the filtering effect of the volume conductor. It is however possible to increase both the 866 number of identified MUs and the ratio of identified low-threshold MUs by increasing the size and the 867 electrode density of HDEMG surface grids to obtain experimental samples of identified MUs that are 868 representative of the discharge activity of the complete MU pool (Caillet et al., 2023). 869 The experimental sample of = 81 MUs, that was obtained at 30% MVC with a dense and large grid 870 of 256 electrodes with 4 mm interelectrode distance, was representative of the discharge activity of 871 the complete MU pool because the identified MUs both spanned across the entire recruitment range, 872 with the lowest-threshold identified MU recruited at 0.6% MVC, and were homogeneously distributed 873 across the recruited MU pool, with 12% to 20% of the identified MUs in each of the 5% ranges sampling 874 the recruitment range, as displayed with the histogram in Figure 6C. Consequently, the neural drive to 875 muscle was accurately estimated with the discharge activity of the = 81 MUs ( Figure 6C, Table 3), 876 which produced the most accurate force estimations in this study when input to the neuromuscular 877 model ( Figure 8C, Table 4). Consistent with previous conclusions (Caillet et al., 2023), the experimental 878 samples obtained at 30% MVC with lower electrode density (8 and 12 mm interelectrode distance) 879 provided a less accurate description of the discharge activity of the complete MU pool. Although they 880 also spanned across the entire recruitment range, fewer MUs were identified and their distribution 881 across the MU pool shifted towards relatively larger sub-population of high-threshold MUs (histograms 882 in Figure 6A, B). Consequently, the neural drive to muscle ( Figure 6A, B, Table 3) and the predicted 883 muscle force ( Figure 8A and B, higher 1 and 2 values in Table 4) Figure 6A, B) and the muscle force ( Figure 8A, B) 891 systematically increased (Table 3, Table 4). As expected, this reconstruction step did not improve the 892 accuracy of the predictions for the experimental sample of = 81 MUs, which was already 893 representative of the discharge activity of the MU pool. 894 Because of physiological reasons discussed elsewhere (Caillet et al., 2023), the discharge activity of the 895 low-threshold MUs was not identified during high-force contractions up to 50% MVC, even with high 896 electrode density. Beyond shifting the MU identification towards the highest-threshold MUs and 897 yielding imbalanced MU distributions (histograms in Figure 6D, E), the experimental samples of 898 identified MUs at 50% MVC did not span across the entire recruitment range, the identified MUs being 899 first recruited above 5% and 20% MVC. Consequently, those experimental samples inaccurately 900 predicted the onsets of the neural drive ( Figure 6D, E) and of the whole muscle ( Figure 8D, E) with 1.2 901 to 5.2 seconds delay and initial underestimations of the muscle force up to 268 N (Δ 1 and Δ 1 F in Table  902 4). Because the reconstruction method relies on the neural drive estimated from the experimental 903 sample ( Figure 6D, E), it could not correct for this source of inaccuracy (solid red traces in Figure 6 and 904 experimentally. In such case, the distribution of simulated MU forces ( Figure 7B, C) was not 927 physiological and interpretable. We showed that locating the MUs into the real ℎ -ranked MU 928 pool with EQ. 6 before assigning representative 0, values ( Figure 7C) corrected for the non-929 homogeneous distribution of the experimental samples across the MU pool, which was previously 930 discussed, and allowed predictions accuracies of ( ) close to those obtained with the complete 931 population of MUs (green dashed versus solid red traces in Figure 8, Table 4). For example, the few 932 low-threshold MUs identified with grids of low electrode density were assigned high representative 933 0, values to represent the force-generating properties of the large population of small MUs that 934 were not identified experimentally. When the representative 0, values were assigned under the 935 assumption of a homogeneous distribution of the MUs across the MU pool ( Figure 7B), the 936 predictions of the whole muscle force were systematically less accurate, again under-evaluating the 937 force-generating activity of the low-threshold MU population. 938

Model scaling 939
The neuromuscular model developed in this study was scaled with the muscle-specific distribution of 940 MU maximum isometric forces 0 ( ) derived from the literature (Figure 2) and the muscle 941 architectural parameters (Table 2)  properties (EQ. 8 and EQ. 11 to EQ. 18) were tuned to match external experimental data independently 947 from the EMG and torque signals that were used in this study to control and validate the model. The 948 accurate prediction of the whole muscle force amplitude at 30% MVC in Figure 8E-G relied instead on 949 an estimation of muscle co-contraction during ankle dorsiflexion ( Figure 5), a physiological and 950 comprehensive description of the neural drive, which the normalization of bEMG envelopes with 951 signals recorded at MVC cannot achieve, and an adequate distribution of the MU's maximum isometric 952 forces 0, across the modelled MU population, as discussed. In cases for which subject-specific 953 measurements are not possible or muscle-specific data are not as available in the literature as for the 954 TA, which is the case for most human muscles (Duchateau & Enoka, 2022), missing muscle architectural 955 parameters or the coefficients in EQ. 8 could be included in an optimization routine aiming to minimize 956 the difference between experimental and predicted joint torques, for example. 957 Although the approach presented in Figure 1 provides a state-of-the-art approach for investigating and 989 modelling the dynamics of the MU pool in forward simulations of human voluntary muscle contraction, 990 it suffers from the following limitations. First, the muscle model was simplified to account for a rigid 991 tendon. This approach was acceptable for the TA and decreased the computational load. Yet, the force 992 equilibrium between the elastic tendon and the active muscle tissue should be considered with EQ. 28 993 in Appendix A.2 for muscles with more compliant tendons or higher 0 ratios, in which case ̅ would 994 vary and the PEE would not be neglected anymore. Second, the model is computationally more 995 expensive than the classic single-input Hill-type actuators, which prevents its use in real-time 996 investigations. As previously discussed (Caillet et al., 2022a), the computational load can be decreased 997 by considering smaller reconstructed populations of 50 MUs without loss of accuracy in the estimated 998 During the second experimental session, voluntary ankle torques (black dashed traces in Figure 9) of produced by the TA muscle despite having more than half the force-generating capacities, which is 1102 again consistent with a participant who aimed to isolate the TA during those contractions. 1103 1104 Figure   In the following, all TA quantities are normalized with respect to the MU or whole muscle optimal 1129 length or maximum isometric force and are reported with a bar following EQ. 25. In EQ. 25, ,0 , ,0 , 1130 , and ̅̅̅ are the maximum isometric force, the optimal length, the length, and the contraction 1131 velocity of the FG of MU , respectively. 1132 EQ. 27 cannot be solved without knowing, for each MU, the value of its individual fibre architecture 1139 (optimal fibre length, fibre length), which implies knowing their spatial distribution, at all frames and 1140 active state. The muscle quantities are currently unfeasible to measure and therefore we implemented 1141 a number of simplifications displayed in Figure 10 and detailed in the following. First, it was assumed 1142 that all the MUs in the MU pool shared the same optimal MU length 0 , which was taken to be the 1143 subject-specific muscle optimal length 0 derived previously in EQ. 2. Second, expecting the MU 1144 optimal lengths 0, to be related to the MU lengths , all the MUs in the MU pool were assigned a 1145 common (normalized) length ̅ = 0 and by extension a common normalized contraction velocity ̅ , 1146 as ̅ = ̅ . Third, to isolate ̅ = ̅ and rearrange EQ. 27 as a differential equation of the common 1147 normalized MU length ̅ , the MU FV relationship was simplified to be independent from the 1148 individual MU active state , which yielded EQ. 28. 1149 ̅ = ̅ = −1 ( ̅̅̅̅̅̅ − ̅̅̅̅̅̅ ∑ ,0 ̅̅̅̅̅̅ • ( ( ̅ ) • ( , ̅ )) ) EQ. 28 Although EQ. 28 can be solved numerically for ̅ , the computational cost and the complexity of the 1150 model were further reduced by first assuming the TA tendon to be rigid {4}. As 0 = 3.5 for the 1151 participant's TA muscle, this assumption was judged acceptable for low-to-medium TA isometric 1152 contractions, considering that 0 = 1 describes a 'very stiff' tendon in the literature (Zajac Felix, 1989). 1153 With this assumption, the tendon length is constant and equals the tendon slack length (Millard et al., 1154(Millard et al., 2013 = . As is also constant during isometric contractions, so is the common normalized 1155 MU length ̅ , in which case = 1. Second, because = + 1.16 • 0 for the segmented TA 1156 muscle in the conditions of the experiment (see Results Section for the subject-specific values 1157 calculated with EQ. 2), it was finally approximated {5} that the FGs worked at the common normalized 1158 length ̅ = 1.16, in which case ̅̅̅̅̅̅ = 0.02 ≈ 0, according to EQ. 29 with 1 = 5 and 2 = 0.6 for 1159 human young adults (Thelen, 2003). With these five key simplifications {1-5}, the PEE and SEE (in grey in Figure 1F) were neglected and the 1161 rheological structure was reduced to the cohort of in-parallel MU actuators (in black in Figure 1F), the 1162 dynamics of which reduced to EQ. 10 in the main text of the manuscript 1163 1164 Figure 10: Step-by-step simplification of the MN-driven model. After neglecting the pennation angle between 1165 the TA muscle and tendon, the MUs are first assumed to share a common optimal length, which is the whole