Motoneuron-driven computational muscle modelling with motor unit resolution and subject-specific musculoskeletal anatomy

Abstract

The computational simulation of human voluntary muscle contraction is possible with EMG-driven Hill-type models of whole muscles. Despite impactful applications in numerous fields, the neuromechanical information and the physiological accuracy such models provide remain limited because of multiscale simplifications that limit comprehensive description of muscle internal dynamics during contraction. We addressed this limitation by developing a novel motoneuron-driven neuromuscular model, that describes the force-generating dynamics of a population of individual motor units, each of which was described with a Hill-type actuator and controlled by a dedicated experimentally derived motoneuronal control. In forward simulation of human voluntary muscle contraction, the model transforms a vector of motoneuron spike trains decoded from high-density EMG signals into a vector of motor unit forces that sum into the predicted whole muscle force. The motoneuronal control provides comprehensive and separate descriptions of the dynamics of motor unit recruitment and discharge and decodes the subject’s intention. The neuromuscular model is subject-specific, muscle-specific, includes an advanced and physiological description of motor unit activation dynamics, and is validated against an experimental muscle force. Accurate force predictions were obtained when the vector of experimental neural controls was representative of the discharge activity of the complete motor unit pool. This was achieved with large and dense grids of EMG electrodes during medium-force contractions or with computational methods that physiologically estimate the discharge activity of the motor units that were not identified experimentally. This neuromuscular model advances the state-of-the-art of neuromuscular modelling, bringing together the fields of motor control and musculoskeletal modelling, and finding applications in neuromuscular control and human-machine interfacing research.

Introduction

During voluntary skeletal muscle contractions, the pool of alpha-motoneurons (MNs) and the muscle fibres they innervate, constituting the muscle’s population of motor units (MUs), transform in a series of chemical-mechanical events the neural message from the spinal cord into molecular forces [1]. These forces sum across the multiscale architecture of the muscle to build the individual MU forces and the whole muscle force. Extensive experimental investigations in animals in vivo have advanced our understanding of the neuromechanical dynamics responsible for voluntary muscle contraction. For example, the force generated by skeletal muscles is modulated by the discharge frequencies of the innervating MNs [2], and by the number of discharging MUs, which are sequentially recruited in the increasing order of their size and the force they can generate [3,4]. Because invasive in vivo practices in humans are reduced to specific electrode insertions [5,6] or challenging surgeries [7–9], our understanding of human voluntary muscle contraction mainly relies on non-invasive in vivo techniques, such as bipolar (bEMG) and high-density (HDEMG) surface electromyography [10–12], imaging [13,14], and joint torque recordings. These non-invasive techniques are however challenging to apply during real-time contractions on multiple muscles, provide limited information that is usually not easily interpreted, or are commonly constrained to simple tasks like isometric contractions. When some neuromuscular properties and dynamics cannot be measured in humans, they are often extrapolated from animals, despite important limitations [9,15,16].

To support experimental investigations into understanding and simulating the interplay between human motor control and force generation, mathematical models of skeletal muscles were developed to mimic the dynamics of muscle contraction. Hill-type neuromuscular models are popular solutions which have been extensively reviewed recently [17]. They rely on few parameters, are computationally cheap and conceptually simple, while complex enough to describe the chain of neuromechanical events responsible for muscle force generation, which makes them accurate, flexible and capable of refinement [18]. Hill-type models have been used to explore voluntary neuromuscular control with EMG-driven forward simulations of muscle contraction [19–23]. With these methods, the recorded EMG signals are processed to provide a unique neural control to a single Hill-type actuator that phenomenologically describes whole muscle dynamics. This approach is well-established and implemented in automated tools [21] and has had important applications, for example in human-machine interfacing [24–27].

However, and despite attempts to address this limitation [28–32], this single-input approach lumps the recruitment and discharge dynamics of the muscle’s individual MUs into a unique phenomenological neural control for the whole MU population, the resulting force-generating events of which are hence lumped into single representative quantities for the whole MU pool. These representative quantities usually overlook the continuous distribution of the MUs’ neuromechanical properties across the MU pool and are difficult to calibrate and interpret [32]. Consequently, this single-actuator macroscopic approach does not provide an adequate structure to physiologically describe the independent and inter-related [33–35] dynamics of a muscle’s MU pool. These limitations hinder the investigation of human neuromuscular control with computational tools. In response to this limitation, other studies developed more physiological muscle models described as populations of Hill-type actuators, each of which simulated the dynamics of individual MUs [36–38]. Yet, those models received synthetic motoneuronal signals and were not tested, and most of them not validated in conditions of voluntary neuromuscular control. Recent advances in the recording and decomposition of surface HDEMG signals [39,40] allow the non-invasive identification of trains of action potentials for large samples of MUs discharging during human voluntary contraction. To date, those vectors of identified MU spike trains were systematically compiled into single neural controls to drive single Hill-type muscle actuators [22,23]. Consequently, the experimental spiking behaviour of the identified discharging MNs was never used to control Hill-type-like models of individual MUs in the simulation of human voluntary contraction.

To address this gap, in this study a novel MN-driven neuromuscular model was developed and validated. The model was described as a population of Hill-type actuators, each of which simulated the neuromechanical and force-generating dynamics of individual MUs. In forward simulations of human voluntary contraction, the MN-driven model transformed an experimental vector of MN spike trains decoded from HDEMG signals [40], that comprehensively described the dynamics of MU recruitment and rate coding, into a vector of simulated MU forces that summed into the predicted whole muscle force. Sampling the muscle into individual MUs provided an adequate structure for proposing robust multiscale simplifications, advanced models of the MU activation dynamics, and physiological muscle-specific distributions of the MU’s neuromechanical properties across the MU pool, that were scaled to subject-specific values derived from a subject-specific musculoskeletal (MSK) model. It was shown that accurate muscle force predictions were obtained when the experimental motoneuronal controls accurately described the real discharge activity of the MU pool. This was obtained with experimental samples of MN spike trains derived from large and dense grids of EMG electrodes, and, to some extent, when the discharge activity of the complete MU pool was reconstructed from experimental data with computational methods [41]. By proposing novel solutions for controlling and designing Hill-type-like neuromuscular models, this study advances the state-of-the-art of neuromuscular modelling, opens avenues for investigating the interplay between the central nervous system and the neuromuscular machinery during human voluntary contraction, reconciles the complementary fields of motor control and MSK modelling, and finds applications in numerous fields, including the investigation of the human neuromuscular dynamics and neural synergies [35] during voluntary contractions, and the design and control of neuroprosthetics [42,43]. The implementation of the method is publicly available at https://github.com/ArnaultCAILLET/MN-driven-Neuromuscular-Model-with-motor-unit-resolution. The segmented medical images and the subject-specific MSK model are publicly available at https://zenodo.org/records/10069266.

Methods

1. Overview of MN-driven neuromuscular modelling

We developed an experimental and computational method (Fig 1) for simulating, with a MN-driven neuromuscular model, the neuromechanical dynamics of a subject-specific population of individual MUs during human voluntary isometric muscle contraction. From surface HDEMG signals recorded on the Tibialis Anterior (TA) muscle (Fig 1A and 1B), we decoded the experimental discharge activity of the population of identified MUs (Fig 1C), which was extended to the complete MU pool (Fig 1D) using an open-source computational method [41], available at https://github.com/ArnaultCAILLET/Caillet-et-al-2022-PLOS_Comput_Biol. The accuracy of both the experimental and reconstructed populations of spike trains spi(t) in estimating the neural drive to muscle was validated (Fig 1E). We described the MN-driven neuromuscular model (Fig 1F) as a collection of in-parallel Hill-type Force Generators (FGs), each of which was controlled by a dedicated spike train spi(t). Each FG described the neuromechanics of a MU and included a Neuromechanical Element (NE) and a Contractile Element (CE) to describe the excitation-contraction coupling properties (MU excitation and activation dynamics) and mechanical properties (MU contraction dynamics) of the MU, respectively [17]. The FGs individually transformed the spike trains spi(t) into the MU forces that collectively generated the whole muscle force FM(t) (Fig 1G). We further scaled the MN-driven neuromuscular model to be subject-specific with muscle-tendon properties derived from a subject-specific MSK model obtained from segmented magnetic resonance images (MRIs) (Fig 1H and 1I). We finally validated the predicted force FM(t) against a reference muscle force FTA(t) (Fig 1G), that was estimated (Fig 1J) from the recorded ankle torque T(t) and the bipolar EMG (bEMG) activity of the TA’s agonist and antagonist muscles crossing the ankle joint.

[Figure omitted. See PDF.]

Fig 1. Description of the steps involved in the forward prediction and validation of the TA muscle force during isometric contractions using the developed subject-specific MN-driven neuromuscular model.

During trapezoidal isometric ankle dorsiflexions (A), HDEMG signals were recorded from the TA muscle with a grid of 256 EMG electrodes, built with four grids of 64 electrodes (B), which was artificially downsampled into two grids of lower electrode density. Using convolutive blind-source separation, the HDEMG signals were decomposed into Nr identified motor unit spike trains (C), from which the discharge activity of the complete MU pool was inferred and reconstructed (D) [41]. The experimental and reconstructed populations of spike trains spi(t) were validated for their accuracy in estimating the normalized effective neural drive to muscle (E). The MN-driven neuromuscular model developed in this study (F), consists of a collection of in-parallel Hill-type MU Force Generators (FGs), each including a Neuromechanical Element (NE) and a Contractile Element (CE), that collectively generate the whole TA muscle force FM. FM was predicted (G) from the experimental spi(t) trains in a blind approach (blue dotted lines), after location of the identified MUs into the MU pool (green dashed lines), 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-specific MSK model (I), obtained from segmented MRI scans (H). The predicted TA force , or , depending on the type of neural control, was finally validated (G) against the experimental TA force FTA (J). FTA was estimated from the recorded ankle torque T and the estimated co-activity of the TA’s agonist and antagonist muscles crossing the ankle joint.

https://doi.org/10.1371/journal.pcbi.1011606.g001

In the following, when subject-specific properties could not be measured in vivo, the coefficients of the mathematical equations describing the model were tuned with experimental data from the literature, considering in turn and decreasing order of preference experimental studies on individual human MUs or bundles of fibres, human fibres and sarcomeres, cat fibres, rodent fibres, rodent muscles, and amphibian fibres or muscles. Any experimental quantity measured at the sarcomere and/or fibre scales was directly assigned to the MU scale.

2. Experimental data and subject-specific neuromusculoskeletal quantities

2.1. Contraction tasks, force and EMG recordings, and subject-specific spiking activity.

The experimental approach represented in Fig 1A was described elsewhere [40] and is briefly summarized here. One healthy male participant (age: 26 yr; height: 168 cm; body weight: 52 kg) volunteered to participate in the experimental session of the study. The participant sat on a massage table with the hips flexed at 30°, 0° being the hip neutral position, and his knees fully extended and strapped to the massage table to avoid any knee and hip motion and co-contraction of the thigh muscles during the ankle dorsiflexion. The ankle torque signals T(t) were recorded with a load cell (CCT Transducer s.a.s, Turin, Italy) connected in-series to the pedal of a commercial dynamometer (OT Bioelettronica, Turin, Italy), to which the participant’s foot was fixed at 30° in the plantarflexion direction, 0° being the foot perpendicular to the shank. The effects of gravity, including the rig weight, were cancelled by removing the initial measured offset during the experiments. After registering the participant’s maximum voluntary contraction (MVC) in dorsiflexion, the participant performed two trapezoidal isometric contractions at 30% and 50% MVC with 120 s of rest in between, consisting of linear ramps up and down performed at 5%/s and a plateau maintained for 20 s and 15 s at 30% and 50% MVC, respectively.

To derive subject-specific motoneuronal controls to the neuromuscular model, HDEMG signals were recorded over the TA muscle during a first experimental session with a total of 256 EMG electrodes (Fig 1B), distributed between four rectangular grids of 64 electrodes, with 4 mm interelectrode distance covering 36 cm2 of the muscle surface (10 cm x 3.6 cm; OT Bioelettronica). HDEMG and force signals were recorded using the same acquisition system (EMG-Quattrocento; OT Bioelettronica) with a sampling frequency of 2,048 Hz. As previously described [40], the 256-electrode grid was downsampled by successively discarding rows and columns of electrodes and artificially generating two new grids of lower electrode density (Fig 1B). The two grids covered the same muscle area with interelectrode distance of 8 mm and 12 mm, involving 64 and 36 electrodes, respectively. After visual inspection and band-pass filtering, the HDEMG signals recorded with the three grids were decomposed into MUs spike trains using convolutive blind-source separation, as previously described [39]. After the automatic identification of the MUs duplicates were automatically removed according to a threshold of 30% of discharge times in common between identified spike trains [40], and all the MUs spike trains were visually checked for false positives and false negatives [11]. Only the MUs which exhibited a pulse-to-noise ratio > 28 dB were retained for further analysis. The final Nr spike trains spi(t) identified experimentally (Fig 1C) were stored as binary vectors of zeros and ones identifying the sample times when a discharge occurred. The ankle torques Tth at which the Nr MUs were recruited, also called torque recruitment thresholds in the following, were calculated as the average of the recorded T values over a 10 ms range centred around the MUs’ first identified discharge time. In the following, the Nr MUs identified experimentally were ranked in the ascending order of measured recruitment torques with the index i∈⟦1; Nr⟧. According to the Henneman’s size principle [3,4], these MUs were therefore also ranked in the same order of current recruitment threshold Ith, maximum isometric forces , and innervation ratios IR according to Eq 1.

(1)

In a second session where the isometric contractions were repeated following the same protocol, bEMG signals were recorded from the pair of agonist extensor digitorum longus (EDL) and extensor hallucis longus (EHL) muscles, and from the antagonist Gastrocnemius Medialis (GM) and Lateralis (GL), and Soleus (SOL) muscles with adhesive bipolar electrodes (OT Bioelettronica). During this second session, the participant additionally produced two MVCs in plantarflexion. The bEMG signals were band-pass filtered (10–450 Hz), full-wave rectified, and then low-pass filtered (2 Hz) [22], before the resulting envelopes were normalized with respect to the peak processed EMG values obtained in the MVC trials.

2.2. MRI scans, subject-specific MSK quantities, and experimental force estimation.

To develop a subject-specific model of the MSK system to couple with the neuromuscular model, a 3D T1-weighted VIBE (volumetric interpolated breath-hold examination) sequence was used to acquire high resolution images (0.45 x 0.45 mm pixel size, 1 mm slice thickness and increment) from the knee joint to mid-foot of the participant’s dominant leg. The sequence consisted in three blocks with 50 mm overlap, during which the participant was asked not to move to facilitate a successful merging of the adjacent blocks.

The volumetric shapes of the tibia, fibula, and foot bones, and of the TA muscle belly and tendon were carefully identified by manual segmentation (Fig 1H) using ITK-SNAP [44] and an anatomical atlas as support [45]. Slices of the agonist EDL and EHL and antagonist SOL, GM, and GL muscle bellies were segmented at regular intervals along the tibial length so that their centroidal lines of action could be outlined. All the main tendons crossing the ankle joint could be clearly identified and were fully segmented. An expert radiologist confirmed the accuracy of the segmentation. The segmented shapes were manually adjusted on the transversal plane to avoid inconsistency of the muscle geometry in correspondence of the MRI sequence blocks using Netfabb (https://www.autodesk.com/products/netfabb/) and Meshlab [46]. The TA muscle volume VM computed from the segmentation was found consistent with the relationships between anthropometry and muscle volumes proposed in [47] from a cohort of segmented MRI scans. Therefore, the volume VM of the EDL, EHL, SOL, GM, GL muscle bellies was estimated with the volume-height-mass relationship provided in the supplementary material of [47].

A subject-specific skeletal model of the ankle joint compatible with OpenSim (version 4.4) [48,49] was automatically generated from the segmented bone volumes (Fig 1I) using the open-source automated toolbox STAPLE [50]. With this approach, the ankle joint was modelled as a hinge joint, the axis of rotation of which was personalised based on the shape of the articular surfaces of the segmented talus bone. Using another automated approach [51], the centroidal line of action of TA was automatically computed from the segmented muscle geometry. Similarly, the centroidal lines of action for EDL, EHL, SOL, GM, and GL were manually outlined from the segmented muscle slices. The muscle tendons paths were manually defined in NMSBuilder [52] with an approach previously proposed [53] that uses the segmented tendons as reference.

To reproduce the experimental conditions, the plantarflexion angle of the final subject-specific MSK model was set to 30°. At this angle, the subject-specific muscle-tendon length lMT of the six muscles and their moment arm L with respect to the ankle joint were computed using OpenSim. The subject-muscle-specific optimal length and the tendon slack length of the six muscles were estimated from the muscle-specific values and proposed in a generic published model [54], which were scaled with Eq 2 by maintaining their ratio with respect to the entire musculotendon lengths and lMT in the neutral ankle position.

(2)

Considering a specific tetanic tension of σ = 60 N/cm2 [54], the subject-specific maximum isometric force of the six muscles was obtained from the known muscle volumes VM with Eq 3.

(3)

Combining these subject-muscle-specific properties with the bEMG signals recorded from the co-contracting flexor muscles during the second experimental session (Fig 1A), the experimental force FTA(t) developed by the TA during the first experimental session was inferred (Fig 1J) from the measured experimental ankle torque T(t) with Eq 4. In Eq 4, LTA is the moment arm the TA tendon makes with the ankle joint, and ΔT(T) relates the level of measured ankle torque T(t) to the algebraic amount of ankle torque ΔT taken by the group of agonist EHL and EDL and antagonist SOL, GM, and GL muscles. The continuous ΔT(T) relationship was obtained by trendline fitting the (T; ΔT) cloud of points measured during the second experimental session. Details on the calculation of ΔT(T) are provided in S1 Text (Section 1).

(4)

3. Generic properties of the MU pool in the human TA muscle

From an extensive review of the experimental human TA literature performed elsewhere [41], it was assumed that a typical adult TA muscle included = 200,000 muscle fibres gathered into N = 400 MUs. This review also returned that the TA MU torque recruitment thresholds Tth, expressed in percentage of the maximum recorded torque (% MVC), followed the generic linear-exponential continuous distribution in Eq 5 and displayed in Fig 2A. In Eq 5, j is an integer identifying the jth MU in the entire pool ranked in the ascending order of MU torque recruitment thresholds.

(5)

[Figure omitted. See PDF.]

Fig 2.

In a generic human TA muscle, distributions in the pool of N = 400 MUs of the MU torque recruitment thresholds Tth (A), MU maximum isometric forces (B), MU twitch forces normalized to the muscle maximum isometric force and expressed in percentage (C), and MU innervation ratios IR (D), calculated with Eq 5 to Eq 9. (A) As identified by the dashed line, 231 TA MUs, i.e., 58% of the MU pool is recruited below 20% MVC, which is consistent with previous conclusions [2]. (B) The normalized distribution in Eq 8 was scaled to Newtons with the subject-specific maximum isometric force of the TA derived previously. From the subject-specific moment arm of the TA derived previously, TA MUs are expected to produce 20 to 230∙1−3 Nm twitch torques with this distribution, which is consistent with the measurements reported in [5]. (D) As identified by the dashed line, the 359 smallest TA MUs gather 72% of the fibre constituting the muscle and are assumed slow-type MUs. According to Eq 1, the MUs are ranked in the same ascending order of Fth, and IR.

https://doi.org/10.1371/journal.pcbi.1011606.g002

The Nr MUs identified in the experimental trials were ranked with the index i (i∈⟦1; Nr⟧, with Nr<N) in the ascending order of measured . Using the Tth(j) distribution in Eq 5, these experimental MUs were mapped into the complete Tth-ranked pool of 400 MUs, where they were assigned a new index Ni (Ni∈⟦1; N⟧), as shown in Eq 6. To do so, the equation was solved for the unknown Ni for the Nr experimental MUs.

(6)

Never measured in human muscles in vivo, the continuous distribution of the normalized MU maximum isometric forces in the TA was here approximated from a continuous distribution of MU twitch force normalized to the muscle maximum isometric force . Processing published experimental TA data [5,55,56] with the same method as used for Eq 5, Eq 7 was obtained. In Eq 7, a 16.3-fold range was taken for the f tw values, consistently with Fig 3A in [5] for a single subject, as opposed to a 250-fold range proposed in a recent review [57] that relied on the same figure, but considering unnormalized merged data from ten subjects.

(7)

The normalized distribution in Eq 8 was then estimated by scaling the coefficients in Eq 7 to obtain a 11-fold range for the values in the TA MU pool, to account for the different twitch-tetanus ratios obtained between slow and fast fibres in the rodent literature [58]. In Eq 8, . The distribution in Newtons, displayed in Fig 2B, is finally obtained by multiplying Eq 8 by the value obtained with Eq 3.

(8)

Human muscles also display a continuous distribution of slow-to-fast MU types in the MU pool. It is known that, in a human TA of N = 400 MUs, 72% of the = 200,000 TA fibres are Type 1 ‘slow’ fibres [59]. Considering that the innervation ratio IR of a MU is roughly proportional to the amplitude of the MU force twitch f tw [4], Eq 9 gives the continuous distribution of IR, displayed in Fig 2D.

(9)

Although it is debated [4] whether the MU type, currently determined by the MU twitch force f tw, is also correlated or not to the MU size, i.e., to the MU Force recruitment threshold in humans, it was here assumed that the small low-threshold MUs were of the ‘slow’ type, and the large high-threshold MUs of the ‘fast’ type. Therefore, from the cumulative summation of the MU IRs, the 359 lowest-threshold MUs that have the lowest innervation ratios were considered to account for 72% of and to be of the ‘slow’ type in the following.

4. Development of a MN-driven neuromuscular model described as a population of Hill-type MU actuators

4.1. Rheological description and preliminary simplifications.

The MN-driven neuromuscular model developed in this study (Fig 1F) consists of a population of n in-parallel FGs, that describe the force-generating activity of the MU pool with mathematical models of the individual MU’s force-generating dynamics (Fig 3). n is the length of the vector of available input spike trains controlling the model and takes the values Nr or N in this study. Using step-by-step simplifications justified in detail in S1 Text (Section 2), {1} the TA tendon was assumed rigid (tendon length set to the tendon slack length [60]) and the in-series elastic element (SEE, in grey in Fig 1F) and its passive force dynamics were neglected, {2} all the MU actuators were assigned the same optimal length, set as the muscle’s , and the same normalized constant length, set as according the subject-specific MSK measurements, {3} the Force-Velocity properties of the MUs were therefore neglected, and {4} the passive force of the CE was found negligible, so that the passive elastic element (PEE, in grey in Fig 1F) was neglected. With these assumptions, and after normalizing the force FM and MU length l state variables to the and parameters and reporting them with a bar, as described in S1 Text (Section 2), the whole muscle dynamics in isometric contractions were described with Eq 10. In Eq 10, at time t and for a population of n modelled MUs, FM is the whole muscle force in Newtons, is the normalized force developed by the kth MU (k∈⟦1; n⟧), and is the MU-specific normalized maximum isometric force that the kth MU can develop. In Eq 10, the MU active states and Force-Length (FL) scaling factors were derived from the input spike trains spk(t) and the common MU length with six cascading mathematical models of the MU’s excitation, activation, and contraction dynamics summarized in Fig 3. Finally, because the MUs build force upon the discrete activity of their fibres that receive the same MN action potential with a random time delay αc, the calculation of the normalized MU forces in Eq 10 was updated to to account for this non-simultaneous twitch activity, which acts as a low-pass filter.

(10)

[Figure omitted. See PDF.]

Fig 3. For a representative motor unit (MU) of the pool of Hill-type MUs in Fig 1F, description of the cascading dynamics that transform an input spike train sp(t) into a normalized MU force .

The time-history of is obtained with Eq 10 by product of the scaling factor fFL(t), the MU active state a(t), and the MU normalized isometric force . fFL is obtained with a normalized FL relationship (Fig 4A) that entirely describes the Contraction Dynamics of the Contractile Element (CE) in isometric conditions. The Neuromechanical Element (NE) is controlled by the binary spike train sp(t), which predicts during the MU’s Excitation Dynamics the trains of motoneuron (MN) and muscle Unit (mU) Action Potentials (APs) in Volts (Fig 4B), identified with the e(t) and u(t) state variables, respectively. The discharge activity u(t) triggers the MU’s activation dynamics by successively controlling, with the c(t), P(t) and a(t) state variables, the concentration in Mols of free calcium ions ([Ca2+]) in the sarcoplasm (Fig 4C–4F), the concentration in Mols of the Calcium-Troponin complex (CaTn) in the sarcomeres, and the MU active state a. The interplay between the MU’s NE and CE is modelled with a nonlinear dependency (dashed arrows) of the MU’s Activation and Contraction Dynamics to the MU’s length and activation a(t), respectively. The MU’s force-generating capacities also depend on some MU-specific properties (dotted lines), including its type (slow/fast) and its normalized isometric force .

https://doi.org/10.1371/journal.pcbi.1011606.g003

For improved readability, the subscript k was discarded in the following of the Methods section to refer to the dynamics of a representative MU of the pool.

4.2. Contractile element and contraction dynamics—MU force-length scaling factor fFL.

The MU fFL scaling factor in Eq 10 was defined with a normalized FL relationship and entirely describes the MU’s contraction dynamics (Fig 3). To the authors’ knowledge, no FL relationship was measured at the MU scale in humans. Despite limitations related to the multiscale approach [17], a human Force–Sarcomere Length relationship obtained from measurements in skinned human fibres [7] was here considered the most suitable option. The sarcomere-to-MU multiscale assumption is here acceptable as muscles modelled as a scaled sarcomere for their normalized isometric FL properties can provide accurate Hill-type model predictions for mammalian whole muscles [61]. Other available solutions were considered less appropriate, such as torque-length and torque-angle measurements in humans which are affected by muscle-tendon dynamics and muscle co-contraction, or FL measurements from fibre bundles in cat, rabbit and rat specimens, that have different optimal sarcomere lengths than humans [62]. Rather than a piecewise linear relationship consistent with the sliding filament theory at the sarcomere scale [63], a gaussian mathematical description in Eq 11 was fitted to the FL experimental data (red curve and black dots in Fig 4A) to account for the smoother relationship with larger plateau observed at larger scales [36,64]. For higher physiological credibility and improved prediction accuracy [65], the MU FL relationships were made nonlinearly dependent to the MU activation states a(t) [20] in Eq 11 at submaximal contraction levels (dashed lines in Fig 4A).

(11)

[Figure omitted. See PDF.]

Fig 4. Contraction, Excitation, and Activation Dynamics (Fig 3) of the slow and fast motor units (MUs) constituting the MN-driven neuromuscular model.

(A) Normalized FL relationship describing the scaling factor fFL. The experimental data points (dots), obtained from human skinned fibres [7], were fitted (solid red curve) with Eq 11. fFL is nonlinearly dependent on the MU active state a at submaximal levels (dashed lines) [20], where the apparent optimal MU length increases by up to 15% (dotted line). (B) Time-history of the MN (black curves) and fibre (blue curves) normalized APs simulated with Eq 12 and Eq 13, respectively. The solid lines were obtained after tuning the coefficients in those equations (Table 1) to fit cat experimental data, while the dotted lines, reported for comparison, were obtained from simulations using the original amphibian data as in Hatze’s work [28]. (C-F) Time-course of the free calcium concentration c in μM for slow (C, E) and fast (D, F) MUs. The calcium transients simulated after tuning the coefficients in Eq 14 to match experimental rodent data (Table 1) are displayed with solid black lines. The simulated transients were validated against other experimental [78,79] traces (green dotted lines) at 67Hz (C) and 125Hz (D) stimulation, respectively. In (E, F), the simulated transients were compared against those proposed in Hatze’s work [28] for amphibians (red dotted curves). (G, H) Relationship between the amplitude of the simulated active state a and the MN firing frequency in Hz. (G): for slow-type MUs, time-history of the active state a(t) with increasing MN firing frequency at 2 (blue), 10 (orange), 15 (green), 20 (red), 25 (purple), 30 (brown) and 50 Hz (pink) where the activation twitches fuse to reach 1.0 at steady-state. (H) Activation-Frequency relationship in the steady-state. Simulations for slow and fast-type MUs are displayed in blue triangles and red symbols respectively. Experimental torque-frequency relationships obtained from FDL (dotted trace, [87]), EDB (solid trace, [86]) and thenar (dashed trace, [85]) muscle fibres are superimposed. (I, J) The bottom row displays the time-course of the five state variables describing the dynamics of the Neuromechanical Element when the MU receives a unique nerve impulse (I) and a train of nerve impulses at 30Hz (J). Black and blue solid curves: trains of MN and muscle unit action potentials e(t) and u(t). Green dashed curve: free calcium concentration c(t). Red dotted curve: CaTn concentration P(t). Purple solid curve: MU active state a(t). Except for the MU active state, the quantities were normalized to their maximum physiological values for visual purposes.

https://doi.org/10.1371/journal.pcbi.1011606.g004

4.3. MU excitation dynamics.

Dynamics of MN AP elicitation. The time-history of the MU active state a(t) in Eq 10 was inferred by first transforming the input binary train of MN spikes sp(t) into a train e(t) of MN action potentials (APs) (Fig 3). The MN APs were phenomenologically modelled with Eq 12 as analytical sine waves [28] of half period and voltage amplitude Ve. The waves are triggered at the discharge times ti when sp(t = ti) = 1, while the MN membrane remains at equilibrium, i.e., e(t) = 0 for . The individual AP shapes recorded in the literature of anesthetized and dissected cats previously reviewed [17] were fitted with Eq 12, yielding Ve = 90mV and T = 1.4ms (Table 1). T is here shorter (Fig 4B) than the amphibian T = 2ms [66] used in a previous neuromuscular model [28].

(12)

[Figure omitted. See PDF.]

Table 1. Coefficients scaling Eq 12 to Eq 18 that describe the MUs’ Excitation and Activation dynamics for both slow-type and fast-type MUs.

Those coefficients were tuned to match contemporary experimental data on mammals at body temperature.

https://doi.org/10.1371/journal.pcbi.1011606.t001

Dynamics of fibre AP elicitation. The train of the fibre APs (u(t) in Fig 3) was modelled in Eq 13 as the 2nd-order response [28] to the train e(t) of MN APs. In Eq 13, when the fibre receives an isolated MN discharge, a1 determines the peak-to-peak amplitude Ap of the fibre AP, while a2 and a3 respectively determine its time-to-peak duration ttp and exponential decay time constant τd. The fibre APs simulated with the tuned ai coefficients in Table 1 matched typical experimental cat values for Ap (77 mV), ttp (0.6 ms) and τd (0.6 ms) [67–77], and were five times shorter (Fig 4B) than previously proposed for amphibians at low temperature [28].

(13)

4.4. MU activation dynamics.

Dynamics of free calcium concentration. The MU excitation state u(t) was modelled to control the level c(t) of calcium concentration [Ca2+] in the MU sarcoplasm (Fig 3) with the 2nd-order ordinary differential equation (ODE) in Eq 14, although alternative reviewed approaches were previously proposed [17]. In Eq 14, when the fibre receives an isolated fibre AP, b1 determines the peak-to-peak amplitude Ac of the [Ca2+] twitch, and b2 and b3 its time-to-peak duration ttp and decay time constant τd. The equation proposed in [28] was modified to introduce the length-dependent scaling factors and , defined in Eq 15 and Eq 16, that respectively make Ac and τd nonlinearly dependent to the normalized MU length . The procedure taken to tune the ci coefficients and derive the scaling factors to match available experimental data for amphibian and mammal slow and fast fibres in vitro is reported in S1 Text (Section 3).

(14)(15)(16)

As reported in Fig 4D for fast fibres, the calcium transients simulated with the tuned bi coefficients (Table 1) compared well with the experimental trace at 125 Hz stimulation from [78] once steady-state was reached. The calcium twitches at steady-state returned the expected ttp = 2.0ms, τd = 6.6ms, and Ac = 20μM with less than 5% difference. As Eq 14 does not account for some physiological mechanisms, such as receptor saturation and the CASQ-Traidin-RYR1 action, the typically larger Ac and τd of the first experimental [Ca2+] twitch [79] were underestimated. Experimental data for slow type rodent fibres were available in [79] at 16°C, but not at 35°C body temperature. In [78], Ac, ttp and τd respectively increased by 26%, remained constant, and decreased by 56% when increasing the temperature from 16°C to 35°C for fast-type fibres. The Ac, ttp and τd coefficients reported at 16°C [79] were scaled to 35°C using the same scaling factors. As shown in Fig 4C, the simulated calcium transients at 35°C for slow fibres under 67Hz stimulation compared well with the scaled experimental data with similar ttp, and τd of 19 and 44ms and an average underestimation of Ac of 20%. The calcium twitches simulated at optimal fibre length with Eq 14 and the parameter values in Table 1 have ten times shorter half-widths and 2.3 times lower amplitudes (Fig 4E and 4F) than the twitches obtained with Hatze’s model [28], which was calibrated on frog data at 9°C.

Dynamics of Calcium-Troponin concentration. The sarcomeric calcium-troponin (CaTn) binding-unbinding process was described with Eq 17 to predict from the free calcium concentration c the time-course of the concentration P of the Ca-Tn structures in the sarcomeres (Fig 3). The initial description [80] was extended to account for the nonlinear length-dependency of the CaTn transients with the scaling factors. In Eq 17, P0 is the total troponin concentration in the myoplasm, and c1, c2 are respectively the forward and backward reaction rates of the Ca-Tn binding-unbinding process. When triggered by a unique calcium twitch, scales the peak-to-peak amplitude At of the CaTn twitch, and and determine its time-to-peak ttp and decay constant τd. The P0, c1, c2 coefficients, originally calibrated to match experimental force data with a neuromuscular model [80], were here tuned to match contemporary experimental data [79,81] following a method described in S1 Text (Section 4).

(17)

The tuned P0, c1, c2 coefficients reported in Table 1, with which the experimental At, ttp, and τd quantities could be reproduced with Eq 17 with less than 5% error at all experimentally measured lengths and for both fibre types, are strongly consistent with the typical coefficient values reported in the cited literature c1 = 3.6∙1012 M−2∙s−1, c2 = 100 s−1, and P0 = 2.2∙10−4 M.

Dynamics of MU activation. The MU active state can be defined as the normalized instantaneous ratio of formed cross-bridges in the force-generating state in a population of myosin heads from a representative sarcomere [17]. Although semi-phenomenological models of multi-state cross-bridge attachments [82] are available, a phenomenological model [83] was used to infer the MU active state a from the CaTn concentration P (Fig 3) to reduce the computational cost, the number of parameters to tune, and the overall model complexity. In Eq 18, d1, d2 and d3 respectively control the amplitude Aa of the activation twitch, its time-to-peak ttp and its time of half-relaxation t0.5, and were tuned to match, lacking more suitable experimental measurements in the literature, the normalized experimental twitch torque or force recorded in individual human MUs in the TA muscle [5,55,56] following a method described in S1 Text (Section 5).

(18)

Because of the MU type-specificity of the [Ca2+] and CaTn dynamics modelled previously, a unique set of values was obtained for the tuned d1, d2 and d3 coefficients (Table 1) for both slow and fast MUs, with less than 7% error in reproducing the experimental At, ttp and t0.5 quantities when simulating single activation twitches for both slow and fast MUs.

When receiving artificial excitatory spike trains at 2, 10, 15, 20, 25, 30, and 50 Hz, the NE of a slow-type MU at optimal length produced the activation traces in Fig 4G. As supported by previous findings in mammals [2], the simulated activation twitches fused for MN firing frequencies above 2Hz and the initial rate of increase of the MU active state was independent from discharge frequency, while the 0-to-0.8 rate of increase increased with higher firing frequencies from 3.5/s at 15Hz and 5.8/s at 50Hz. Consistent with previous findings [84], every new elicited activation twitch increased the MU activation level at steady-state by 100% at 5Hz, and only by 1% at 50Hz, following a negative exponential tendency with increasing MN discharge frequency.

Fig 4H displays the simulated steady-state activation-frequency relationship for both slow-type MUs (blue triangles) and fast-type MUs (red symbols). Despite obvious limitations but lacking experimental measurements of human activation states in the literature, those results were compared to the literature of experimental human muscle force and torque twitches. In Fig 4H, simulated twitch-tetanus ratios of 0.30 and 0.19 were obtained at steady-state for the simulated slow and fast-type MU active states, respectively. These ratios are of the same order of magnitude as those experimentally obtained for human hand muscles (range 0.1–0.3, [85–87]). The activation-frequency results for the simulated slow-type MUs compared relatively well with the findings from these experimental studies with a maximum difference of 30% at 15Hz with the extensor digitorum brevis (EDB) muscle data [86] and less than 10% at all physiological frequencies above 5Hz with the data obtained from thenar muscles [85]. Yet, the simulated rate of increase in a with firing frequency largely underestimated the available experimental data for the fast-type MUs with up to 200% difference. The lack of available experimental human data of fast MU twitches prevents for now from deriving more advanced conclusions on the accuracy of the simulated MU active states.

4.5. Solving the dynamics of the MU pool.

Working at fixed MU length for each of the n modelled MUs, the cascading ODEs in Eq 12 to Eq 18, that describe the MU-specific excitation and activation dynamics, were first numerically solved at each time step, yielding a vector of MU active state time-histories ai(t). Fig 4I and 4J displays the time-histories of the five state variables that describe the NE’s activity at 1 and 30 Hz, respectively. An estimation of the MU normalized forces and of the whole muscle force was then obtained with Eq 10.

5. Neural control to the MN-driven neuromuscular model

5.1. Neural controls to the neuromuscular model.

In this study, we assessed the MN-driven neuromuscular model (Fig 1F) with two different types of neural controls, i.e., the experimental samples of Nr identified spike trains (Fig 1C) and the complete population of N MU spike trains (Fig 1D) obtained with a MU pool reconstruction method [41]. Both approaches were tested with the sets of spike trains identified with the three electrode densities (Fig 1B) at both contraction intensities up to 30% and 50% MVC.

When directly inputting the Nr experimental spike trains to the neuromuscular model, the Nr MUs were assigned with Eq 19 representative maximum isometric forces , necessary parameter in Eq 10, to account for the force-generating capacities of their neighbouring MUs, in the sense of their recruitment thresholds, that were not identified experimentally. The muscle force (blue dotted line in Fig 1G) was first predicted when the Nr identified MUs were assumed to be homogeneously spread across the Tth-ranked MU pool, in which case the distribution in Eq 8 was evenly split in Nr domains along the x-axis (Fig 2B). A second TA force (green dashed line in Fig 1G) was predicted after assigning to the Nr MUs representative maximum forces with Eq 19 after a physiological mapping of these MUs to the complete Tth-ranked MU pool using Eq 6. For example, when Nr = 30 MUs are identified at 30% MVC, where Na = 300 MUs are supposed to be recruited (Eq 5), the 9th, 10th, and 11th identified MUs in the Tth-ranked experimental sample are mapped into the complete Tth-ranked pool to the 90th, 100th, and 110th MUs with the first method (considering even steps of MUs) and to the 90th, 120th, and 180th MUs with the second method (using Eq 6). The 10th identified MU represents 10 MUs and produces the combined force with the first method, while it represents 45 surrounding MUs and generates with the second method.

(19)

The TA force (red solid line in Fig 1G) was predicted from a comprehensive cohort of N = 400 simulated spike trains, that described the discharge activity of a completely reconstructed MU population (red solid line in Fig 1E). The complete population of N MUs was obtained from a validated computational method [41] that maps with Eq 6 the Nr MUs into the MU pool and infers from the Nr identified spike trains a continuous distribution of the MN electrophysiological properties across the MU pool to physiologically scale a cohort of N = 400 models of motoneurons that predicts the discharge activity of the N−Nr MUs that were not identified experimentally. The implementation of the method is publicly available at https://github.com/ArnaultCAILLET/Caillet-et-al-2022-PLOS_Comput_Biol. Here, the distribution of MU maximum forces in Eq 8 was directly applied to the reconstructed population of N MUs.

5.2. Neural drive to muscle and quality of the experimental neural control.

Cumulative spike trains were computed as the temporal binary summation of the pools of spike trains and were low-pass filtered in the bandwidth [0–4] Hz relevant for force generation [88] to approximate the effective neural drive to muscle D(t). As suggested for isometric contractions [33], the normalized effective neural drive was compared for validation against the normalized experimental TA force (black trace in Fig 1E) with calculation of three metrics: the error in seconds in identifying the onset of the neural drive Δ1, the normalized root-mean-square error (nRMSE) and the coefficient of determination r2. Of note, was obtained by normalization to the average of the FTA(t) values over the plateau of contraction and not to the maximum registered FTA value to minimize the impact of possible discharge artefacts onto the interpretation of the nRMSE metric. To discuss the importance of the quality of the experimental data in deriving a reliable neural control to the MN-driven neuromuscular model developed in this study, the validation between the and traces was performed for both the pools of Nr experimental and N reconstructed spike trains (blue dotted and plain red traces in Fig 1E, respectively) derived from the three grids of increasing electrode density (Fig 1B).

6. Validation of the predicted forces

The TA forces FM(t) predicted with the MN-driven neuromuscular model (Fig 1F) were validated against the experimental TA force FTA(t) (solid black trace in Fig 1G), that was estimated from the recorded ankle torques T(t) with Eq 4. The FM(t) versus FTA(t) validation was performed with calculation of {1} the error in seconds in identifying the onset Δ1 of generated force, i.e. when 2% of the max generated force was produced, {2} the experimental force developed after the Δ1 delay, {3} the maximum error (ME) in Newtons, {4} the nRMSE in percentage, and {5} r2. The nRMSE was calculated for the complete contraction, the ascending ramp (nRMSEr1), the plateau (nRMSEp), and the descending ramp (nRMSEr2). Δ1 and assess the delay in predicting the onset of force and only depends on the first discharge time identified experimentally [41], ME identifies the maximum error in predicting the whole muscle force and is related to the Δ1 and metrics in the regions of low generated force, the nRMSE metric evaluates the accuracy in predicting the amplitude of the predicted whole muscle force, which mainly relies on both the distribution of MU tetanic forces and the accurate prediction of the MU active states, and r2 evaluates the time-course of the variation of the muscle force, which mainly relies on the distribution of the MU neuromechanical properties across the MU pool.

Results

1. Subject-specific MSK quantities and experimental TA force

The subject-specific MSK model (Fig 1I) yielded the subject-muscle-specific properties in Table 2 for the TA muscle and for its agonist EDL and EHL, and antagonist muscles SOL, GM, and GL in ankle dorsiflexion. By reapplying with Eq 4 the exponential ΔT(T(t)) relationship derived from the bEMG and torque recordings obtained during the second experimental session to the torque traces T(t) recorded during the first experimental session, the experimental TA force FTA(t) (solid red traces in Fig 5) was estimated (see S1 Text (Section 1) for details of the calculation). In the normalized space (bottom row in Fig 5), compared well to the muscle force that would be computed if the TA alone produced the total recorded ankle torque T (black dashed traces in Fig 5), with a maximum error ≤ 15%, nRMSE ≤ 7%, and r2 ≥ 0.98. Yet, ignoring the co-contracting activities of the EDL, EHL, SOL, GM, and GL muscles would underestimate the TA force FTA by up to 49 N and 283 N, i.e., 18% and 60% of the maximum FTA value, at 30% and 50% MVC respectively (upper row in Fig 5).

[Figure omitted. See PDF.]

Fig 5.

Estimation of the force developed by the TA muscle during the trapezoidal dorsiflexions up to 30% (left) and 50% MVC (right) performed during the first experimental session. The TA forces (first row: in Newtons, second row: normalized), displayed in solid red traces, were estimated by correcting the recorded torque T for co-contraction with the relationship in Eq 4 (See S1 Text (Section 1) for details about the calculation). The black dashed traces display the TA force calculated without considering the co-contraction of surrounding muscles (as if the TA alone produced the entire recorded ankle torque T).

https://doi.org/10.1371/journal.pcbi.1011606.g005

[Figure omitted. See PDF.]

Table 2. Summary of the subject-muscle-specific MSK properties of the six flexor muscles when the ankle joint is set to a 30° plantarflexion angle.

The moment arm L made by the tendons with the ankle joint and the muscle-tendon lengths lMT were identified by MRI segmentation and measurements in NMSBuilder [52] and OpenSim [48,49]. The optimal length and the tendon slack length were scaled from the generic values reported in [54] with the ratio of generic to subject-specific muscle-tendon lengths lMT at default 0° plantarflexion angle. The muscle maximum isometric force was estimated from the muscle volume VM, which was obtained from MRI segmentation for the TA, and from anatomical relationships [47] for the other muscles. The normalized muscle lengths , at 30° angle in ankle plantarflexion and assuming a rigid tendon, were calculated with Eq S1 in S1 Text (Section 1).

https://doi.org/10.1371/journal.pcbi.1011606.t002

2. Experimental neural control

For the contractions up to 30% and 50% MVC, 81, 32, and 14 MUs, and 55, 15, and 3 MUs (Fig 1C) were respectively identified [40] with the EMG grids of 4, 8, and 12 mm interelectrode distance (Fig 1B). The dataset of three identified MUs was disregarded in the following results because too few and only high-threshold MUs were identified in this dataset, which therefore provided a highly inaccurate description of the full spectrum of MU properties and discharge activity of the real MU pool [41]. The histogram distribution of the identified MUs across the MU pool according to their recruitment threshold Tth in %MVC is displayed in Fig 6. Importantly, the denser the EMG grid and the lower the generated force, the more homogeneous the distribution of the identified MUs across the MU pool [40]. For example, 12% to 20% of the MUs decoded with the densest grid at 30% MVC were systematically identified in each of the 5% ranges sampling the recruitment range (histogram in Fig 6C). Higher heterogeneity and lower representation of low-threshold MUs were obtained with shallower grids and at 50% MVC. Classically, the decomposition algorithms converge towards the large MUs that contribute most to the EMG signals, while action potentials of the small MUs of lower energy are masked by the potential of larger units, more of which are recruited at 50% MVC than at 30% MVC. At 30% MVC, increasing the electrode density better samples the action potential profiles of the low-threshold unmasked MUs across multiple electrodes, enabling their identification. The normalized effective neural drive to muscle (Fig 1E) was computed with these cohorts of spike trains and was validated against the experimental normalized TA force (black solid traces in Fig 6) with the three validation metrics reported in Table 3.

[Figure omitted. See PDF.]

Fig 6.

For the three grid configurations of increasing electrode density (in columns) and for both contractions up to 30% MVC (first row: A, B, C) and 50% MVC (second row: D, E), validation against the experimental normalized force trace (solid dark trace) of the normalized effective neural drive to muscle estimated from the completely reconstructed population of N = 400 MUs (solid red traces) and from the experimental samples of Nr identified MUs (dotted blue traces). The TA force FTA(t) (Fig 5) was estimated from the recorded ankle torque T with Eq 4. The histogram distributions of the identified MUs across the recruitment range Tth in %MVC were normalized to the number of MUs identified experimentally.

https://doi.org/10.1371/journal.pcbi.1011606.g006

[Figure omitted. See PDF.]

Table 3. Validation of the normalized effective neural drive to muscle for both contractions up to 30% and 50% MVC, computed from the spike trains obtained with the three grid configurations of 4, 8, and 12 mm interelectrode distance that involve 256, 64, and 36 EMG electrodes respectively.

The validation was performed for both the experimental samples of Nr identified MUs, and for the completely reconstructed pool of N MUs.

https://doi.org/10.1371/journal.pcbi.1011606.t003

The normalized effective neural drive estimated from the experimental samples of Nr MUs (blue dotted lines in Fig 6) was always better approximated when the electrode density increased, especially in the regions of low generated forces. For example, at 30% MVC, it was obtained Δ1 = 0.07 s, nRMSE = 8.0%, and r2 = 0.97 from the 81 discharging MUs identified with 4 mm interelectrode distance, while the 14 MUs identified with 12 mm interelectrode distance returned Δ1 = 0.14 s, nRMSE = 17.5%, and r2 = 0.89 (Table 3). The effective neural drive estimated from the experimental MU samples was always better approximated at 30% MVC than at 50% MVC (blue dotted traces in upper row versus lower row in Fig 6), with two to three times lower nRMSE values, and strongly lower Δ1 and higher r2 values, respectively (Table 3).

Complete populations of N = 400 MNs (Fig 1D) were also reconstructed from the discharge activity of the Nr identified MNs by estimating the discharge activity of the MUs not identified experimentally [41]. In this case, the estimated at 30% MVC (solid red lines, upper row, in Fig 6) was accurately predicted and was unrelated to the quality of the input experimental samples, with nRMSE<7%, and r2 = 0.98 for the three EMG grid configurations (Table 3). At 50% MVC, the derived from the N MNs (solid red lines, bottom row, in Fig 6) was better estimated when the MU pool was reconstructed from denser grids of electrodes, with lower nRMSE and higher r2 values (Table 3). As for the experimental MU samples, estimated from the N = 400 MUs was better approximated at 30% MVC than at 50% MVC (solid red traces in upper row versus lower row in Fig 6). At 50% MVC, the reconstructed populations of discharging MUs returned inaccurate estimations of the onset and ascending ramp of force, with Δ1>1s, nRMSE>12%, and r2≤0.95. Overall, the normalized effective neural drive was systematically better approximated by the reconstructed pool of N = 400 MUs than by the experimental samples of Nr identified MUs (solid red versus blue dotted traces in Fig 6) with lower nRMSE and higher r2 values (Table 3).

Therefore, was best approximated by samples of MU spike trains of identified MUs that spanned across the whole recruitment range and followed a homogeneous distribution across the MU pool. Obtaining such representative description of the discharge behaviour of the real MU pool was possible with EMG signals recorded at 30% MVC with high electrode density (4 mm interelectrode distance) or, to some extent, with the computational reconstruction of the complete MU pool, both of which aimed to correct for the systematic bias of HDEMG decomposition towards identifying the higher-threshold MUs (histograms in Fig 6) of the discharging MU pool [40].

Activation and Contraction dynamics of the MU pool and TA Force validation

The neuromuscular model was built as populations of Nr or N in-parallel FGs, whether it received as neural control the experimental samples of Nr identified spike trains (Fig 1C) or the completely reconstructed populations of N MU spike trains (Fig 1D). According to Fig 2, around 298 and 355 MUs were recruited over the trapezoidal contractions up to 30% and 50% MVC, respectively, in which case all the modelled discharging MUs were assigned slow-type properties for their activation dynamics (Fig 3).

Fig 7A displays the distribution of the maximum values of the modelled MUs’ activation states over the contraction up to 30% MVC, that were estimated from the cohorts of Nr = 81 (green triangles) and N = 400 (red dots) spike trains with Eq 12 to Eq 18. Consistent with the onion skin theory [89] where MN discharge rate and recruitment threshold are negatively related, the estimated MU activation was overall negatively related to recruitment threshold, with low-threshold MUs reaching maximum activation states up to a = 0.80 and higher-threshold MUs as low as a = 0.18. Because of trendline fitting during the reconstruction process of the MU pool [41], the monotonically decreasing distribution of active states of the N-population did not account for the physiological variability in MN discharge rate with recruitment threshold (green triangles in Fig 7A). The same conclusions were obtained with the neural controls obtained with the grids of lower EMG electrode densities and at 50% MVC.

[Figure omitted. See PDF.]

Fig 7.

Distribution across the Tth-ranked MU population of the maximum values at 30% MVC of the MU activation states ak (A) and MU forces (B-D) simulated with the MN-driven neuromuscular model. The model either received as neural input the Nr = 81 experimental spike trains (green triangles in A for ak; plots B and C for ) or the discharge activity of the complete pool of N = 400 MUs (red dots in A for ak; plot D for ). In B-D, the black dots represent the distribution of the MUs’ maximum forces , which is continuous for the completely reconstructed pool (D). In B and C, the Nr identified MUs were assigned representative maximum forces (black dots) to account for the force-generating properties of the MUs not identified experimentally. To do so, in B, the Nr identified MUs were assumed to be homogeneously distributed across the Tth-ranked MU pool, while they were accurately mapped to the real pool with Eq 6 in C. In B, the non-continuous distribution of the values is explained by the fact that is not an integer value in Eq 19.

https://doi.org/10.1371/journal.pcbi.1011606.g007

Fig 7B–7D displays the distribution of the maximum forces (blue crosses, green triangles, red dots) developed by the modelled MUs over the contraction up to 30% MVC, that were estimated with Eq 10 from the MUs’ active states a(t), the common MU length , and the MU-specific maximum isometric forces (black dots), for cohorts of Nr = 81 (Fig 7B and 7C) and N = 400 (Fig 7D) MUs. Depending on the modelling approach, the modelled MUs were assigned maximum MU forces (black dots in Fig 7B–7D) by either blindly assuming the Nr MUs to be evenly distributed across the MU pool (Fig 7B), mapping the Nr MUs to the real MU pool with Eq 6 (Fig 7C), or directly applying the distribution in Eq 8 to the complete MU pool (Fig 7D). Depending on the approach, the individual MUs produced maximum forces over the contraction up to 30% MVC in the range 1.8–6.6 N, 0.4–16.5 N, and 0.5–1.5N (blue crosses in Fig 7B, green triangles in Fig 7C, and red dots in Fig 7D, respectively), and in the range 4.2–16.2 N, 1.5–136.7 N, and 0.7–2.2 N at 50% MVC. Only the reconstructed pool of N = 400 MUs (Fig 7D) provided a window onto the continuous distribution of the dynamics of the MU pool with physiological generated MU forces.

The simulated MU forces were processed and linearly summed with Eq 10 to yield the three predicted whole muscle forces in Fig 8A–8E (experimental MU population with blind distribution—blue dotted traces), (experimental MU population with MU mapping—green dashed traces), and (reconstructed MU population—solid red traces), which were validated against the experimental TA force FTA(t) (solid black traces) with the validation metrics in Table 4.

[Figure omitted. See PDF.]

Fig 8.

For the three grid configurations of increasing electrode density (in columns) and for both contractions up to 30% MVC (first row: A, B, C) and 50% MVC (second row: D, E), validation against the experimental force trace FTA(t) (solid dark trace) of the TA force FM(t) estimated with the MN-driven neuromuscular model in Fig 1F. FM(t) was estimated from the completely reconstructed population of N = 400 MUs (, solid red traces) and from the experimental samples of Nr identified MUs when ‘blind’ (, blue dotted traces) and Tth-informed (, green dashed traces) distributions of the MU maximum isometric forces were assigned to the MU samples. The experimental TA force FTA(t) was estimated from the recorded ankle torque T with Eq 4.

https://doi.org/10.1371/journal.pcbi.1011606.g008

[Figure omitted. See PDF.]

Table 4. Validation against the experimental force trace FTA(t) of the predicted TA force FM(t) at 30% and 50% MVC.

FM(t) was predicted from the spike trains obtained at 4-, 8-, and 12-mm electrode density. FM(t) was estimated from the completely reconstructed population of N = 400 MUs () and from the experimental samples of Nr identified MUs when ‘blind’ () and Tth-informed () distributions of the MU maximum isometric forces were assigned to the MU samples. nRMSE, nRMSEr1, nRMSEp, nRMSEr2 were calculated for the whole contraction, and over the ascending ramp, plateau, and descending ramp, respectively.

https://doi.org/10.1371/journal.pcbi.1011606.t004

When high numbers of MUs were homogeneously identified over the full range or recruitment (Nr = 81 MUs, 30% MVC, 4 mm interelectrode distance), i.e., when was accurately estimated from the experimental samples of spike trains (Fig 6C), FTA was predicted by , and with the same level of accuracy (Fig 8C) for the eight validation metrics in Table 4 with r2≥0.98, all nRMSE values below 10%, and ME in 36–66%, irrespective from the three choices of assignment (black dots in Fig 7B–7D).

In all other cases, FTA(t) was the most accurately predicted by both when the complete MU pool was modelled (solid red traces in Fig 8A–8E, nRMSE in range 8–13% and r2 in range 0.97–0.99 at 30% MVC, and nRMSE in range 15–16% and r2 in range 0.90–0.95 at 50% MVC) and when the population sample of Nr MUs was assigned representative values according to a MU mapping (green dashed traces in Fig 8A–8E, nRMSE in 6–11% and r2 in 0.97–0.99 at 30% MVC, and nRMSE in 14–20% and r2 in 0.92–0.95 at 50% MVC). The least accurate predictions were systematically obtained with when the population sample of Nr MUs received experimental neural inputs and was blindly assigned values (blue dotted traces in Fig 8A–8E, nRMSE in 7–18% and r2 in 0.89–0.98 at 30% MVC, and nRMSE in 27–32% and r2 in 0.86–0.92 at 50% MVC). and performed better than by especially over the ramps of force with around 1.5–3 times lower ME, nRMSEr1, and nRMSEr2 values (nRMSEr1 and nRMSEr2 metrics calculated over the ramps of force—Table 4), because the blind distribution of values in (Fig 7B) did not correct, contrary to the two other approaches, for the experimental bias towards identifying relatively more high-threshold MUs with decomposed HDEMG signals and overestimated their force-generating activities in the MU pool. As long as the Nr decoded MUs were homogeneously spread across the identified portion of the MU pool, returned slightly more accurate predictions of FTA than , with slightly lower nRMSE and maximum error values over the whole contraction. Mainly, returned higher nRMSEp values (Table 4) calculated over the plateau of force because of the noisy force-generating activity of the highest-threshold recruited MUs modelled with the reconstruction method [41]. With experimental samples of lesser quality (50% MVC, 8 mm interelectrode distance, Fig 6), performed better than .

For all three modelling approaches at both contraction levels, FTA was more accurately predicted according to all the validation metrics when the electrode density increased, especially during the ramps of force. For example, at 30% MVC, it was obtained Δ1 = 0.1 s, N, nRMSE in 6–8%, and r2 in 0.98–0.99 from the 81 discharging MUs identified with 4 mm interelectrode distance, while the 14 MUs identified with 16 mm interelectrode distance returned Δ1 = 0.3 s, N, nRMSE in 11–18%, and r2 in 0.89–0.98 (Table 4).

As expected from the neural drive estimation in Fig 6, for all the modelling approaches and EMG grid configurations, FTA was more accurately predicted at 30% MVC than at 50% MVC (upper row versus lower row in Fig 6), with two to three times lower nRMSE values, and strongly lower Δ1 and and higher r2 values, respectively (Table 4).

Discussion

Summary of the work

This study reports a novel approach, summarized in Fig 1, to develop a MN-driven neuromuscular model, that is controlled at the level of the individual MU with experimental motoneuronal data. This model advances the recently reviewed [17] state-of-the-art of neuromuscular modelling on multiple aspects, and finds later-discussed applications in the complementary fields of motor control and MSK modelling. The population of Hill-type MU actuators (Fig 1F) transforms, with advanced models of the MU’s excitation, activation, and contraction dynamics (Figs 3 and 4), a vector of MN spike trains into a vector of transient MU forces that sum across the modelled MU population to yield the whole muscle force (Fig 8). The neuromuscular model and the continuous distribution of the MU’s recruitment and force-generating properties across the MU pool were respectively made subject-specific using a subject-specific MSK model derived from medical images (Fig 1H and 1I), and muscle-specific, using results from the literature (Fig 2). The neural control to the MN-driven model is also subject-specific as the vector of input MN spike trains is derived from HDEMG signals recorded during the participant’s voluntary task (Fig 1A and 1B). The experimental vector of Nr identified spike trains resulting from HDEMG signal processing was either directly input to the neuromuscular model of Nr MUs (Fig 1C) or first extrapolated to an estimation of the discharge activity of the complete population of N MNs [41], which was then input to the neuromuscular model of N MUs (Fig 1D). The accuracy in estimating the neural drive to muscle of this novel MN-driven approach was assessed using experimental datasets of varying quality (Fig 6), with the control vectors of Nr and N MN spike trains (Fig 1E). Finally, the whole muscle force predicted by the MN-driven neuromuscular model was validated against an experimental muscle force (Figs 1G and 8), that was estimated from measured joint torque and bEMG signals recorded from agonist and antagonist muscles (Figs 1J and 5).

Advancing the state-of-the-art in neuromuscular modelling

The neuromuscular model developed in this study was built as a collection of in-parallel MU actuators, which proved beneficial for advancing the recently reviewed [17] state-of-the-art of neuromuscular modelling on several aspects.

For the first time in Hill-type modelling, we controlled the individual MUs of a modelled MU population with a vector of dedicated experimental motoneuronal controls. The vector of MN spike trains obtained from HDEMG signals provided a comprehensive description of both the dynamics of MU recruitment and MU rate coding for the modelled population of MUs. Conversely in EMG-driven models of whole muscle actuators, the recruitment and discharge dynamics of the MU pool are lumped into a single phenomenological macroscopic neural control where they become indistinguishable, that is either obtained from bEMG envelopes [19,20,90] or from the temporal summation of filtered cumulative spike trains and signal residuals from HDEMG signals [22,23,91]. Besides, in multiscale models of single representative MUs, the dynamics of MU recruitment were either overlooked or accounted for with phenomenological models [28–32,36]. A few studies also assigned motoneuronal controls to modelled populations of MUs [36–38,92–94] but used synthetic data and/or phenomenological models, such as Fuglevand’s formalism [95], to describe the discharge and recruitment dynamics of the MU pool. Consequently, these models of MU populations were not used in conditions of voluntary muscle contraction and the force they predicted was indirectly validated against results from other models and not against synchronously recorded experimental data like performed in this study. A recent study [96] implemented a neuromechanical model of a population of twitch-type models of MUs controlled by experimental vectors of MU spike trains. This model was tested with limited samples of few experimental MN spike trains (i.e., 2–22 identified MUs per contraction per muscle); such samples are typically not representative of the real neural drive to muscle, which is a key challenge in MN-driven neuromuscular modelling, as discussed in the present study (Fig 6). In this perspective, it is worth highlighting that our study did not perform a calibration step of the muscle model parameters, contrary to [96], to avoid any kind of error cancellation that would non-physiologically correct for this experimental limitation when validating the model. Then, [96] proposes a linear summation of the fusing twitches in a MU, that contradicts some experimental evidence discussed in the present study. The physiological nonlinear summation of the MU twitches during muscle contraction, classically modelled in twitch-type models with a nonlinear gain that is MU-specific and depends on the instantaneous firing rate [95], is in our study predicted by the detailed models of the MU activation dynamics described in Eqs 14–18 (Fig 4G and 4H). Finally, the inter-relation in a MU pool proposed in [96] between MU firing properties, MU twitch contraction time, and MU type, as well as their bimodal distribution in the MU pool, remain unproven in human muscles with contradictory findings [5,55,84,97], as reviewed [98] and discussed [4] previously. For these reasons, the MN-driven approach developed in the present study is the first that can control with experimental vectors of inputs, that are representative of the real neural drive to muscle, a population of physiologically accurate models of the individual MUs’ force-generating dynamics in the forward prediction of human voluntary muscle contraction. Moreover, the proposed MN-driven model while also relies on a detailed description of the excitation-contraction dynamics of the active muscle tissue and on a set of experimental spike trains that is more easily interpretable than the phenomenological neural controls used in single-input multiscale models of whole muscle actuators in the literature.

Discretely sampling the active muscle tissue into individual MUs provided a convenient framework for the precise distribution of the muscle’s properties across the MU pool, including the novel experimental distributions in Fig 2 of MU recruitment threshold, maximum isometric forces, and type, some of which were never included in neuromuscular models. This description of the muscle as a cohort of MUs also provided a convenient structure for deriving physiological assumptions and simplifications to the rheological structure of the model (see S1 Text (Section 2) for details).

This MU-scale approach was ideal for developing advanced models of the MU’s excitation, activation, and contraction dynamics (Fig 4). Experimental data from the literature was used in this study to derive adequate mathematical equations (Eq 11 to Eq 18), identify fitting physiological parameters (Table 1), and validate the modelled dynamics (Fig 4). Importantly, the identification of the parameter values in Eq 11 to Eq 18 and the validation of those mathematical descriptions were performed using different sets of experimental data from the literature that were measured at different discharge rates (see Appendices A.1 and A.1 for details). Furthermore, we limited the amount of multiscale and inter-species scaling inaccuracy, which is an important limitation of single-actuator Hill-type models [17], by considering source experimental data, in turn and decreasing order of preference, from studies on individual human MUs or bundles of fibres, human fibres and sarcomeres, cat fibres, rodent fibres, rodent muscles, and amphibian fibres or muscles. In doing so, the simulated MN APs, fibre APs, and free calcium twitches were more physiological and consistent with mammalian dynamics than previous approaches based on experimental amphibian data at low temperature (Fig 4B, 4E and 4F). For the first time in Hill-type-like modelling, the dynamics of calcium-troponin concentration in the sarcoplasm were described according to experimental data (Eq 17) (see S1 Text (Section 4) for details) and the free calcium transients were made nonlinearly dependent to both MU length and MU type based on experimental measurements (Eq 14 to Eq 16). Hence, this model further reconciled the phenomenological Hill-type modelling approach with the real physiological mechanisms responsible for muscle force generation it describes.

While this study presents novel techniques to develop state-of-the-art neuromuscular models as MN-driven pools of MU Hill-type actuators, it remains to assess the performance of these more physiological and complex models compared to the more common phenomenological single-input single-actuator approaches, like the standard bEMG-driven Hill-type models [21]. Almost no studies that proposed models of MU pools performed this comparative assessment, which was therefore never reviewed [99]. Only three studies compared, for limited contraction tasks, the force outputs of a model of MU pool and of a single muscle-scale model, using Hill-type [100] and twitch-type [91,96] approaches. [91,100] reached the same conclusions. The models of MU pool and the single-actuator muscle-scale models can yield equally accurate predictions of muscle force in the time domain, although the MU pool model always performs better in the frequency domain. Contrary to single-actuator models, models of MU pool can predict the higher frequency force fluctuations, crucial in the analysis of steadiness for example. This is explained by the superiority as control signals of vectors of MU spike trains, immune to waveform cancellation, over rectified and smoothed bEMG signals. It is worth noting that the single-actuator Hill-type model in [100] is significantly more complex and elaborate in its phenomenological description of the MU dynamics than traditional EMG-driven Hill-type models, as reviewed [17], while the single-actuator twitch-type model in [91] was fully calibrated to match the force traces. It is therefore possible that MN-driven models of MU pool behave more accurately than more standard single-actuator bEMG-driven approaches, especially between contractions where the neural strategy changes, e.g., between different contraction intensities or when fatigue occurs, considering that calibrated bEMG-driven models are known to behave poorly at other activation levels [101]. In this respect, cumulative spike train driven single-actuator models [22], that describes the MU pool recruitment and firing strategies, were shown to more accurately predict muscle forces than classic bEMG-driven Hill-type approaches.

Consequently, the field requires a systematic assessment of the comparative performance of standard single-actuator neuromuscular models and models of MU pools in a variety of contraction tasks and in generic, subject-specific, and calibrated approaches. Finally, beyond the prediction accuracy, MN-driven models of MU pools provide detailed insights into the muscle’s internal dynamics and finds practical applications that standard single-actuator bEMG-driven models cannot provide, as later discussed.

Modelling limitations

Although the approach presented in Fig 1 provides a state-of-the-art approach for investigating and modelling the dynamics of the MU pool in forward simulations of human voluntary muscle contraction, it suffers from the following modelling limitations. First, the muscle model was simplified to account for a rigid tendon. This approach was acceptable for the TA and decreased the computational load. Yet, the force equilibrium between the elastic tendon and the active muscle tissue should be considered with Eq S9 in S1 Text (Section 2) for muscles with more compliant tendons or higher ratios, in which case would vary and the PEE would not be neglected anymore. Second, physiologically describing the TA muscle with a MU resolution required assigning to the MU pool physiological distributions of neuromechanical parameters (Fig 2) obtained from TA-specific experimental measurements from the literature. Although multiple studies also investigated these parameters in the human thenar, ankle extensor, and masseter muscles, those properties remain mostly unknown for the other human muscles, as recently reviewed [57]. Third, the pipeline in Fig 1 was only applied to one subject in this study. Yet, the approach is general and can be reapplied to other subjects following the same method. Fourth, we had to correct a limitation of our reconstruction method [41] for the simulated MUs recruited during the plateau of constant recorded torque. When the discharge frequency of such simulated MU was below 7 Hz in average over the plateau, the predicted MU force was set to zero to avoid adding random noise to the predicted whole muscle force at the highest developed force levels.

Experimental neural control and prediction accuracy

The accuracy of the whole muscle force predictions in Fig 8 was mainly related to the accuracy of the experimental neural control in approximating the neural drive to muscle (Fig 6). It is currently impossible to identify the discharge activity of the complete MU pool of a muscle with surface EMG because of the filtering effect of the volume conductor. It is however possible to increase both the number of identified MUs and the ratio of identified low-threshold MUs by increasing the size and the electrode density of HDEMG surface grids to obtain experimental samples of identified MUs that are representative of the discharge activity of the complete MU pool [40].

The experimental sample of Nr = 81 MUs, that was obtained at 30% MVC with a dense and large grid of 256 electrodes with 4 mm interelectrode distance, was representative of the discharge activity of the complete MU pool because the identified MUs both spanned across the entire recruitment range, with the lowest-threshold identified MU recruited at 0.6% MVC, and were homogeneously distributed across the recruited MU pool, with 12% to 20% of the identified MUs in each of the 5% ranges sampling the recruitment range, as displayed with the histogram in Fig 6C. Consequently, the neural drive to muscle was accurately estimated with the discharge activity of the Nr = 81 MUs (Fig 6C and Table 3), which produced the most accurate force estimations in this study when input to the neuromuscular model (Fig 8C and Table 4). Consistent with previous conclusions [40], the experimental samples obtained at 30% MVC with lower electrode density (8 and 12 mm interelectrode distance) provided a less accurate description of the discharge activity of the complete MU pool. At this force level, increasing the electrode density better samples the action potential profiles of the low-threshold MUs across multiple electrodes, enabling their identification. In these datasets, although the identified MUs also spanned across the entire recruitment range, fewer MUs were identified and their distribution across the MU pool shifted towards relatively larger sub-population of high-threshold MUs (histograms in Fig 6A and 6B). Consequently, the neural drive to muscle (Fig 6A and 6B and Table 3) and the predicted muscle force (Fig 8A and 8B, higher nRMSEr1 and nRMSEr2 values in Table 4) were underestimated in the regions of low-force generation when low-threshold MUs are recruited. This experimental under-representativity of the low-threshold MUs was corrected using a computational method [41] by deriving the continuous distribution of the MN’s electrophysiological properties across the entire MN pool and reconstructing the discharge activity of the MUs that were not identified experimentally (Fig 1D). The reconstruction method populated the two experimental samples with simulated populations of low-threshold MUs, in which case the accuracy in predicting both the neural drive (blue dotted versus solid red traces in Fig 6A and 6B) and the muscle force (Fig 8A and 8B) systematically increased (Tables 3 and 4). As expected, this reconstruction step did not improve the accuracy of the predictions for the experimental sample of Nr = 81 MUs, which was already representative of the discharge activity of the MU pool.

The discharge activity of the low-threshold MUs was not identified during high-force contractions up to 50% MVC, even with high electrode density. Classically, the decomposition algorithms converge towards the large MUs that contribute most to the EMG signals, while action potentials of the small MUs of lower energy are masked by the potential of larger units, more of which are recruited at 50% MVC than at 30% MVC [40]. Beyond shifting the MU identification towards the highest-threshold MUs and yielding imbalanced MU distributions (histograms in Fig 6D and 6E), the two experimental samples of identified MUs at 50% MVC did not span across the entire recruitment range, the identified MUs being first recruited above 5% and 20% MVC. Consequently, those experimental samples inaccurately predicted the onsets of the neural drive (Fig 6D and 6E) and of the whole muscle (Fig 8D and 8E) with 1.2 to 5.2 seconds delay and initial underestimations of the muscle force up to 268 N (Δ1 and in Table 4). Because the reconstruction method relies on the neural drive estimated from the experimental sample (Fig 6D and 6E), it could not correct for this source of inaccuracy (solid red traces in Figs 6 and 8), although it improved the predictions by better distributing the MU’s discharge activity across the identified MU pool, as discussed. It is worth noting that low-threshold MUs are usually not identified with grids of low electrode density at 30% MVC either, as hinted with the histograms in Fig 6A and 6B, and the resulting experimental samples of MUs hence identified usually do not span across the entire recruitment range, yielding similar prediction inaccuracies. To address this limitation when controlling single whole muscle actuators with HDEMG signals, previous studies [22,23,91] summed the experimental cumulative spike train with the EMG signal residual that was not explained by the identified MU spike trains to approximate the neural control to muscle. Consequently, the accuracy of the predictions made by the neuromuscular model developed in this study are sensitive to the number, the span, and the distribution across the MU pool of the experimental sample of identified MUs, which can, to some extent, be extrapolated to the discharge activity of the complete MU pool using our published reconstruction method [41].

As displayed in Fig 8A–8E, the accuracy of the predicted forces FM(t) also depended on the assignment of the maximum MU isometric forces to the modelled MU pool. The values were obtained from the muscle-specific distribution in Eq 8, which was directly mapped to the N modelled MUs when the complete MU pool was reconstructed. This physiological approach implemented the continuous distribution of the MU forces across the MU pool (Fig 7D) and almost systematically returned the most accurate FM(t) predictions (Fig 8A–8E and Table 4). When the muscle was described as populations of Nr MUs controlled by the experimental vector of MN spike trains, representative values were derived from the distribution and assigned to the experimental MUs to account for the force-generating properties of the MUs not identified experimentally. In such case, the distribution of simulated MU forces (Fig 7B and 7C) was not physiological and interpretable. We showed that locating the Nr MUs into the real Tth-ranked MU pool with Eq 6 before assigning representative values (Fig 7C) corrected for the non-homogeneous distribution of the experimental samples across the MU pool, which was previously discussed, and allowed predictions accuracies of FM(t) close to those obtained with the complete population of N MUs (green dashed versus solid red traces in Fig 8 and Table 4). For example, the few low-threshold MUs identified with grids of low electrode density were assigned high representative values to represent the force-generating properties of the large population of small MUs that were not identified experimentally. When the representative values were assigned under the assumption of a homogeneous distribution of the Nr MUs across the MU pool (Fig 7B), the predictions of the whole muscle force were systematically less accurate, again under-evaluating the force-generating activity of the low-threshold MU population.

Interfacing the fields of subject-specific motor control, neuromuscular modelling, and MSK modelling

The neuromuscular model developed in this study brings together state-of-the-art experimental and modelling techniques from the three complementary fields of motor control, neuromuscular modelling, and MSK modelling.

The subject-specific neural control was obtained from recently developed experimental and processing techniques (Fig 1A–1C) that yielded the identification of much larger samples of discharging MUs from recorded HDEMG signals than commonly obtained in the literature [11], i.e., up to 81 MUs in the TA muscle, and high ratios of low-threshold MUs [40].

The neuromuscular model, the modelling novelty of which was discussed previously, was scaled with the muscle-specific distribution of MU maximum isometric forces derived from the literature (Fig 2) and from the muscle architectural parameters (Table 2), that were obtained from a subject-specific MSK model (Fig 1H and 1I) built from the segmentation of MRI scans using state-of-the-art automated tools [50,51]. MSK model predictions can be sensitive to the uncertainties in parameter identification and MSK architecture [102,103]. The subject-specific properties of the TA muscle used in this study were compared to those from a generic published model [54], and the same model scaled to the anthropometry of the subject with the Opensim built-in tools [48,49]. The values of the optimal length and tendon slack length did not vary between the three models (<3% variation). The subject-specific muscle-tendon length lMT at 30° plantarflexion was 3% shorter and 12% longer than in the scaled and generic models, respectively. Considering the simplification of rigid tendon, these differences in length linearly propagate to the estimation of the normalized MU length and nonlinearly to the FL and length-dependent activation dynamics of the MUs (Fig 4). The subject-specific MRI-based maximum isometric force was 15% lower than proposed in the generic model. This difference linearly propagates to the linear scaling of the distribution of MU isometric forces in Eq 8 (Fig 7B–7D) and to the amplitude in N of the predicted whole muscle force (Fig 8). The highest difference was obtained for the muscle moment arm with the ankle joint in dorsiflexion, where the subject-specific quantity was 33% and 40% lower than the scaled and generic ones. This difference linearly propagates to the estimation of experimental muscle force from experimental joint torque (Fig 5) that was used for model validation (Fig 8).

The model proposed in this study and the predictions displayed in Fig 8 did not rely on any parameter calibration that would minimize external cost functions, like commonly proposed in EMG-driven approaches where some subject-muscle-specific properties are, for example, calibrated by minimizing the difference between predicted and experimental joint torques [20,21].

Beyond the MSK parameters directly measured with the subject-specific MSK model, it is worth noting that the remaining parameters in the equations that define the model’s force-generating properties (Eq 8 and Eq 11 to Eq 18) are fitting experimental data from the literature, that are independent from the EMG and torque signals that were used in this study to control and validate the model.

With this subject-specific and MN-driven approach, the accurate prediction of the whole muscle force amplitude at 30% MVC in Fig 8E–8G were consequently obtained without parameter calibration owing to a physiological and comprehensive description of the neural drive (Fig 6), which the bEMG envelopes normalized by MVC signals cannot achieve, an adequate distribution of the MU’s maximum isometric forces across the modelled MU population (Fig 7), and an estimation of muscle co-contraction during ankle dorsiflexion (Fig 5), as discussed.

In cases for which subject-specific measurements are not possible or muscle-specific data are not as available in the literature as for the TA, which is the case for most human muscles [57], missing muscle architectural parameters or the coefficients in Eq 8 could be included in an optimization routine aiming to minimize the difference between experimental and predicted joint torques, for example.

The only calibrated parameter in the pipeline in Fig 1 was the size of the experimentally identified MNs, that was required for the specific test case of complete reconstruction of the MN pool (Figs 6 and 8, red traces) to derive the continuous distribution of electrophysiological properties across the population of MN LIF models. Again, this parameter identification was independent from the experimental measurements of joint torque used in this study to validate the model. The two other approaches taken to predict the muscle force solely relied on the experimental samples of identified MNs (Fig 8, blue and traces) and did not require the calibration of the MN size parameter.

Applications and future perspectives

The MN-driven multi-MU model proposed in this study comprehensively describes the recruitment and firing strategy of the MU pool during human voluntary muscle contractions. For these reasons, this detailed model finds applications in answering scientific questions where the lump dynamics of the classic single-actuator bEMG-driven Hill-type models provide limited help.

For example, current research [42,43] aims for a new generation of high-performance multi-articulating prosthetic limbs, which, besides relying on direct skeletal attachment via osseointegration, muscle reinnervation, and implanted sensors, also require advanced MU-based algorithms, like the one proposed in this study, that implements the mechanistic relationship between the individual motoneuronal activity and motor function.

In the field of biomechanics, the presented MU-based model can actuate the recently developed volumetric representations of muscles in MSK models [51] by mapping the modelled population of MUs to the volumetric population of lines of action, that are consistent with the segmented muscle geometry. Volumetric representations of Hill-type-actuated muscle models have only rarely been proposed in the literature [36,104]. Beside addressing the current limitations of modelling muscles as single rectilinear segments, this volumetric mapping would also provide a solution to the indeterminacy in in S1 Text (Section 2) by assigning MU-specific lengths to the modelled MU actuators and would shift the current approach towards a more physiological nonlinear summation of the MU forces based on muscle architecture. Modelling the interaction between individual MU lines of action with the tendon and skeletal structures would further bridge the gap in modelling the interplay between motor control and resulting human motion.

In the study of the human muscle architecture, by combining the aforementioned mapping of the MU pool in volumetric muscle representations with recently developed ultrasound measurement techniques for tracking MU twitches and mechanical properties [34,105], one could gain insights into the distribution of MU territories within human muscles, which remains an open question. The volumetric mapping of MUs also provides a convenient framework for modelling the transversal mechanical interactions between MUs [106], and the resulting force-varying load path within the muscle tissue, that results in a nonlinear summation of the individual MU forces.

In the field of neurophysiology, the MN-driven model of MU population proposed in this study is suited for integrating MN synergies in simulations. Contrary to single-input Hill-type models that lump the dynamics of synaptic input, MN recruitment, and MN discharge into a single phenomenological signal, the MU pool modelled in Fig 1 can be divided into functional clusters and receive different common inputs to reduce the dimensionality of the control [35]. For example, in multi-muscle MSK models with MU-actuated volumetric representations of muscles, the model proposed in Fig 1 is also suited for the investigation of neural synergies between muscles, where MN clusters span across muscles [35], with strategies specific to the muscle groups [107]. In the field of neurophysiology, the model of MU population proposed in this study is also applicable for the computational generation of surface EMG signals during voluntary movement [108].

In the study of human muscle contraction dynamics, the model developed in this study provides a credible window onto the distribution of the excitation-contraction dynamics of individual MUs across the MU pool (Fig 4) that cannot be measured in human in vivo, and advances our understanding of the muscle-specific neuromechanical strategies for muscle force generation. For example, Fig 4 suggests that, at 30% MVC, the highest-threshold recruited MUs with the largest force-generating capacities are not those producing the highest forces within the recruited pool due to low activation states explained by the onion skin theory (Fig 8A and 8D). As the model integrates realistic descriptions of the MU-specific firing, recruitment, excitation, and activation dynamics responsible for muscle force generation, it also provides a convenient framework for integrating and investigating MU-related mechanisms, such as the effect of fatigue on the MU firing and recruitment strategies [109] or other fatigue-induced changes in MU activation dynamics as already proposed in single-MU models [110,111].

However, some currently investigated limitations can prevent the proposed MN-driven MU-based model from being readily applicable to the aforementioned scientific questions. First, although the pipeline in Fig 1 is fully suited for dynamic contractions, the study proposed here was constrained to isometric tasks due to current experimental limitations in the decomposition of HDEMG signals into MU spike trains during motion. New computational approaches are emerging to identify individual spike trains in quasi-static [112] and dynamic [113] tasks, where 7–20 MUs can be currently identified. Second, data acquisition remains challenging and time-consuming. As opposed to straightforward bEMG recording and filtering, the pipeline in Fig 1 requires decomposing HDEMG signals and manually editing the identified spike trains editing, while taking precautions to identify the full spectrum of discharging MUs for accurate predictions (Fig 6), as discussed. These challenges should be addressed with the rapid emergence of guidelines for manufacturers on HDEMG grid design [40], open-source tools for automatic HDEMG decomposition and spike trains edition [114], automated spike train identification approaches based on machine learning [115] and blind-source separation methods [116], and MU pool reconstruction methods [41].

Finally, with four ODEs to solve at each time step for each of the modelled MUs, the model proposed in this study is currently computationally too expensive for real-time applications, as opposed to the classic single-input Hill-type actuators, which usually include zero to two ODEs. For the simulation of a 40s-long muscle contraction using a standard laptop (RAM: 12 GB, CPU: one Intel Core i7-1165G7 2.80 GHz), the pipeline in Fig 1 implemented in Python ran in 18 minutes for a MU pool of Nr = 16 MUs (Fig 8A, dotted trace), in 63 minutes for Nr = 81 MUs (Fig 8C, dotted trace), and in 240 to 300 minutes for the complete pool of N = 400 MUs (Fig 8, red traces), that additionally required preliminary MU pool reconstruction. For comparison, a single-actuator rigid-tendon Hill-type model [60] ran in less than a minute on the same machine for the same 40s-long trapezoidal isometric contraction where the HDEMG signals were averaged and normalized to MVC signals. Besides implementations of the model in a faster compiled programming language (for example C++) and using better performing machines, the computational speed could be drastically increased by using parallel computing to solve the independent MU dynamics, by considering smaller reconstructed populations of 50 MUs without loss of accuracy in the estimated neural drive to muscle, as previously discussed [41], or by, for example, simplifying the current model of MU excitation and activation dynamics with analytical descriptions of the fibre APs and calcium transients instead of the two 2nd-order ODEs in Eq 13 and Eq 14.

Conclusion

We developed the first MN-driven neuromuscular model of a population of individual Hill-type MUs controlled by a vector of dedicated experimental motoneuronal controls (Fig 1). The model distinguishes the dynamics of MU recruitment from rate coding and produces the whole muscle force as the summation of the forces generated by the individual modelled MUs. The model is subject-specific (Table 2), muscle-specific (Fig 2), and includes an advanced and physiological model of the MUs’ activation dynamics (Figs 3 and 4 and Table 1). The motoneuronal controls, derived from HDEMG signals, are experimental and decode the subject’s intention, which makes the neuromuscular model applicable to the simulation and investigation of human voluntary muscle contraction. The model’s predictions of the whole muscle force are sensitive to the quality of the experimental neural control. Accurate force predictions were obtained when the effective neural drive to muscle was accurately estimated from the decoded MN spike trains (Figs 6C and 8G), i.e., when the experimental samples of identified MUs were representative of the discharge activity of the complete MU pool. This was obtained when the muscle’s myoelectric activity was recorded with large and dense grids of EMG electrodes during medium-force contractions, in which case the identified MUs span across the complete range of recruitment and are homogeneously distributed across the MU pool (Figs 6 and 8). Otherwise, the discharge activity of the low-threshold MUs is typically not identified, especially during high-force contractions, and the force predictions are inaccurate in the regions of low-force generation (Table 4). Inferring with a computational method the discharge activity of those MUs that were not identified experimentally improves the results to some extent and provides a window onto the continuous distribution of the MUs’ force-generating dynamics across the MU pool. The accuracy of the force predictions also relies on a physiological assignment of the MU-specific force-generating properties to the modelled population of Hill-type MUs (Fig 8B–8D). This MN-driven model advances the state-of-the-art of neuromuscular modelling, brings together the interfacing fields of motor control and MSK modelling, and finds applications in numerous fields, including the investigation of the human neuromuscular dynamics during voluntary contractions, neural synergies, and human-machine interfacing. The implementation of the method is publicly available at https://github.com/ArnaultCAILLET/MN-driven-Neuromuscular-Model-with-motor-unit-resolution. The segmented medical images and the subject-specific MSK model are publicly available at https://zenodo.org/records/10069266.

Supporting information

S1 Text. Supporting methods and information for the derivation and the validation of the neuromuscular model.

https://doi.org/10.1371/journal.pcbi.1011606.s001

(PDF)

Citation: Caillet AH, Phillips ATM, Farina D, Modenese L (2023) Motoneuron-driven computational muscle modelling with motor unit resolution and subject-specific musculoskeletal anatomy. PLoS Comput Biol 19(12): e1011606. https://doi.org/10.1371/journal.pcbi.1011606

[/RAW_REF_TEXT]

References

1. Huxley N. Structural changes in muscle during contraction: Interference microscopy of living muscle fibres. Nature 1954, 173:971–973. pmid:13165697

2. Heckman CJ, Enoka RM. Motor unit. Compr. Physiol. 2012, 2:2629–2682. pmid:23720261

3. Henneman E. Recruitment of motoneurons: The size principle. Motor Unit Types, Recruitment and Plasticity in Health and Disease 1981, 26–60.

4. Caillet AH, Phillips AT, Farina D, Modenese L. Mathematical relationships between spinal motoneuron properties. eLife 2022, 11:e76489. pmid:35848819

5. Van Cutsem M, Feiereisen P, Duchateau J, Hainaut K. Mechanical properties and behaviour of motor units in the tibialis anterior during voluntary contractions. Can. J. Appl. Physiol. 1997, 22:585–597.

6. Muceli S, Poppendieck W, Holobar A, Gandevia S, Liebetanz D, Farina D. Blind identification of the spinal cord output in humans with high-density electrode arrays implanted in muscles. Science advances, 8(46), eabo5040. pmid:36383647

7. Gollapudi SK, Lin DC (2009) Experimental determination of sarcomere force–length relationship in type-I human skeletal muscle fibers. J.Biomech. 2022, 42:2011–2016.

8. Ferguson RA, Shur NF. Skeletal muscle biopsy: Techniques and applications. In Anonymous Sport and Exercise Physiology Testing Guidelines 2022, 205–211.

9. Binder-Markey BI, Persad LS, Shin AY, Litchy WJ, Kaufman KR, Lieber RL. Direct intraoperative measurement of isometric contractile properties in living human muscle. J.Physiol. 2023, 601(10), 1817–1830 pmid:36905200

10. De Luca CJ. The use of surface electromyography in biomechanics. Journal of applied biomechanics 1997, 13:135–163

11. Del Vecchio A, Holobar A, Falla D, Felici F, Enoka RM, Farina D. Tutorial: Analysis of motor unit discharge characteristics from high-density surface EMG signals. Journal of Electromyography and Kinesiology 2020, 53:102426. pmid:32438235

12. Waasdorp R, Mugge W, Vos HJ, De Groot JH, Verweij MD, De Jong N, Schouten AC, Daeichin V. Combining ultrafast ultrasound and high-density EMG to assess local electromechanical muscle dynamics: A feasibility study. IEEE Access 2021, 9:1.

13. Pincheira PA, Boswell MA, Franchi MV, Delp SL, Lichtwark GA. Biceps femoris long head sarcomere and fascicle length adaptations after 3 weeks of eccentric exercise training. Journal of Sport and Health Science 2022, 11:43–49. pmid:34509714

14. Dick TJ, Hug F. Advances in imaging for assessing the design and mechanics of skeletal muscle in vivo. J.Biomech 2023, 111640. pmid:37244210

15. Manuel M, Chardon M, Tysseling V, Heckman CJ. Scaling of motor output, from mouse to humans. Physiology 2019, 34:5–13. pmid:30540233

16. Highlander MM,Allen JM,Elbasiouny SM. Meta-analysis of biological variables’ impact on spinal motoneuron electrophysiology data. Journal of Neurophysiology 2020, 123, 1380. pmid:32073942

17. Caillet AH, Phillips AT, Carty C, Farina D, Modenese L. Hill-type computational models of muscle-tendon actuators: A systematic review. bioRxiv 2022, https://doi.org/10.1101/2022.10.14.512218.

18. Schmitt S, Günther M, Häufle DFB. The dynamics of the skeletal muscle: A systems biophysics perspective on muscle modeling with the focus on hill-type muscle models. GAMM-Mitteilungen 2019, 42:n/a.

19. Hof AL, Van den Berg J. EMG to force processing I: An electrical analogue of the hill muscle model. Journal of biomechanics 1981, 14:747,755–753,758. pmid:7334035

20. Lloyd DG, Besier TF. An EMG-driven musculoskeletal model to estimate muscle forces and knee joint moments in vivo. Journal of Biomechanics 2003, 36:765–776. pmid:12742444

21. Pizzolato C, Lloyd DG, Sartori M, Ceseracciu E, Besier TF, Fregly BJ, Reggiani M. CEINMS: A toolbox to investigate the influence of different neural control solutions on the prediction of muscle excitation and joint moments during dynamic motor tasks. J.Biomech. 2015, 48:3929–3936. pmid:26522621

22. Sartori M, Yavuz U, Farina D. In vivo neuromechanics: Decoding causal motor neuron behavior with resulting musculoskeletal function. Scientific reports 2017, 7:13465–14. pmid:29044165

23. Kapelner T, Sartori M, Negro F, Farina D. Neuro-musculoskeletal mapping for man-machine interfacing. Scientific reports 2020, 10:1–10.

24. Sartori M, Reggiani M, Farina D, Lloyd DG. EMG-driven forward-dynamic estimation of muscle force and joint moment about multiple degrees of freedom in the human lower extremity. PloS one 2012, 7:e52618. pmid:23300725

25. Ao D, Song R, Gao J. Movement performance of human-robot cooperation control based on EMG-driven hill-type and proportional models for an ankle power-assist exoskeleton robot. TNSRE 2017, 25:1125–1134. pmid:27337719

26. Caggiano V, Wang H, Durandau G, Sartori M, Kumar V. MyoSuite—A contact-rich simulation suite for musculoskeletal motor control. arXiv 2022, 2205.13600.

27. Cimolato A, Driessen JJ, Mattos LS, De Momi E, Laffranchi M, De Michieli L. EMG-driven control in lower limb prostheses: A topic-based systematic review. Journal of NeuroEngineering and Rehabilitation 2022, 19:1–26.

28. Hatze H. A myocybernetic control model of skeletal muscle. Biological cybernetics 1977, 25:103–119. pmid:836914

29. Dorgan SJ, O’Malley MJ. A nonlinear mathematical model of electrically stimulated skeletal muscle. T-RE 1997, 5:179–194. pmid:9184904

30. Song D, Raphael G, Lan N, Loeb GE. Computationally efficient models of neuromuscular recruitment and mechanics. Journal of neural engineering 2008, 5:175–184. pmid:18441419

31. Dick TJM, Biewener AA, Wakeling JM. Comparison of human gastrocnemius forces predicted by hill-type muscle models and estimated from ultrasound images. Journal of experimental biology 2017, 220:1643–1653. pmid:28202584

32. Hussein M, Shebl S, Elnemr R, Elkaranshawy H. A new muscle activation dynamics model, that simulates the calcium kinetics and incorporates the role of store-operated calcium entry channels, to enhance the electromyography-driven hill-type models. J.Biomech.Eng. 2021, 144:011002.

33. Farina D, Negro F. Common synaptic input to motor neurons, motor unit synchronization, and force control. Exercise and sport sciences reviews 2015, 43:23–33. pmid:25390298

34. Lubel E, Grandi Sgambato B, Rohlen R, Ibanez J, Barsakcioglu D, Tang M, Farina D. Non-linearity in motor unit velocity twitch dynamics: Implications for ultrafast ultrasound source separation. bioRxiv 2023, 24.533983. pmid:37703141

35. Hug F, Avrillon S, Ibáñez J, Farina D. Common synaptic input, synergies and size principle: Control of spinal motor neurons for movement generation. J.Physiol. 2023, 601:11–20. pmid:36353890

36. Hatze H. Myocybernetic Control Models of Skeletal Muscle—Characteristics and Applications. University of South Africa Press, Pretoria, 1980.

37. Hamouda A, Kenney L, Howard D. Dealing with time-varying recruitment and length in hill-type muscle models. Journal of biomechanics 2016, 49:3375–3380. pmid:27614612

38. Carriou V, Boudaoud S, Laforet J, Mendes A, Canon F, Guiraud D. Multiscale hill-type modeling of the mechanical muscle behavior driven by the neural drive in isometric conditions. Computers in biology and medicine 2019, 115:103480. pmid:31629271

39. Negro F, Muceli S, Castronovo AM, Holobar A, Farina D. Multi-channel intramuscular and surface EMG decomposition by convolutive blind source separation. JNE 2016, 13:026027. pmid:26924829

40. Caillet AH, Avrillon S, Kundu A, Yu T, Phillips Andrew T.M., Modenese L, Farina D. Larger and denser: An optimal design for surface grids of EMG electrodes to identify greater and more representative samples of motor units. eNeuro 2023, ENEURO.0064-23.2023. pmid:37657923

41. Caillet AH, Phillips AT, Farina D, Modenese L. Estimation of the firing behaviour of a complete motoneuron pool by combining electromyography signal decomposition and realistic motoneuron modelling. PLOS Computational Biology 2022, 18:e1010556. pmid:36174126

42. Farina D, Vujaklija I, Sartori M, Kapelner T, Negro F, Jiang N, Bergmeister K, Andalib A, Principe J, Aszmann OC. Man/machine interface based on the discharge timings of spinal motor neurons after targeted muscle reinnervation. Nature biomedical engineering 2017, 1:1–12.

43. Farina D, Vujaklija I, Brånemark R, Bull AMJ, Dietl H, Graimann B, Hargrove LJ, Hoffmann K, Huang H, Ingvarsson T, Janusson HB, Kristjánsson K, Kuiken T, Micera S, Stieglitz T, Sturma A, Tyler D, Weir RFf, Aszmann OC. Toward higher-performance bionic limbs for wider clinical use. Nature biomedical engineering 2021, 1–13.

44. Yushkevich PA, Piven J, Hazlett HC, Smith RG, Ho S, Gee JC, Gerig G. User-guided 3D active contour segmentation of anatomical structures: Significantly improved efficiency and reliability. Neuroimage 2006, 31:1116–1128. pmid:16545965

45. Möller TB, Reif E. Pocket Atlas of Sectional Anatomy. 2007.

46. Cignoni P, Callieri M,Corsini M, Dellepiane M, Ganovelli F, Ranzuglia G. Meshlab: An open-source mesh processing tool. Eurographics Italian chapter conference 2008, 129–136.

47. Handsfield GG, Meyer CH, Hart JM, Abel MF, Blemker SS. Relationships of 35 lower limb muscles to height and body mass quantified using MRI. J.Biomech. 2014, 47:631–638. pmid:24368144

48. Delp SL, Anderson FC, Arnold AS, Loan P, Habib A, John CT, Guendelman E, Thelen DG. OpenSim: Open-source software to create and analyze dynamic simulations of movement. TBME 2007, 54:1940–1950. pmid:18018689

49. Seth A, Hicks JL, Uchida TK, Habib A, Dembia CL, Dunne JJ, Ong CF, Demers MS, Rajagopal A, Millard M, Hamner SR, Arnold EM. OpenSim: Simulating musculoskeletal dynamics and neuromuscular control to study human and animal movement. PLoS computational biology 2018, 14:e1006223. pmid:30048444

50. Modenese L, Renault J. Automatic generation of personalised skeletal models of the lower limb from three-dimensional bone geometries. J.Biomech. 2021, 116:110186. pmid:33515872

51. Modenese L, Kohout J. Automated generation of three-dimensional complex muscle geometries for use in personalised musculoskeletal models. Ann.Biomed.Eng. 2020, 48:1793–1804. pmid:32185569

52. Valente G, Crimi G, Vanella N, Schileo E, Taddei F. nmsBuilder: Freeware to create subject-specific musculoskeletal models for OpenSim. Comput.Methods Programs Biomed. 2017, 152:85–92.

53. Modenese L, Montefiori E, Wang A, Wesarg S, Viceconti M, Mazzà C. Investigation of the dependence of joint contact forces on musculotendon parameters using a codified workflow for image-based modelling. Journal of Biomechanics 2018, 73:108–118. pmid:29673935

54. Rajagopal A, Dembia CL, DeMers MS, Delp DD, Hicks JL, Delp SL. Full-body musculoskeletal model for muscle-driven simulation of human gait. TBME 2016, 63:2068–2079. pmid:27392337

55. Andreassen S, Arendt-Nielsen L. Muscle fibre conduction velocity in motor units of the human anterior tibial muscle: A new size principle parameter. J. Physiol. 1987, 391:561–571. pmid:3443957

56. Van Cutsem M, Duchateau J, Hainaut K. Changes in single motor unit behaviour contribute to the increase in contraction speed after dynamic training in humans. J.Physiol. 1998, 513:295–305. pmid:9782179

57. Duchateau J, Enoka RM. Distribution of motor unit properties across human muscles. J.Appl.Physiol. 2022, 132:1–13. pmid:34709066

58. Celichowski J, Grottel K. Twitch/tetanus ratio and its relation to other properties of motor units. Neuroreport 1993, 5:201–204. pmid:8298075

59. Henriksson-Larsén KB, Lexell J, Sjöström M. Distribution of different fibre types in human skeletal muscles. I. method for the preparation and analysis of cross-sections of whole tibialis anterior. Histochem.J. 1983, 15:167–178. pmid:6343306

60. Millard M, Uchida T, Seth A, Delp SL. Flexing computational muscle: Modeling and simulation of musculotendon dynamics. Journal of biomechanical engineering 2013, 135:021005. pmid:23445050

61. Winters TM, Takahashi M, Lieber RL, Ward SR. Whole muscle length-tension relationships are accurately modeled as scaled sarcomeres in rabbit hindlimb muscles. J.Biomech. 2011, 44:109–115. pmid:20889156

62. Burkholder TJ, Lieber RL. Sarcomere length operating range of vertebrate muscles during movement. J.Exp.Biol. 2001, 204:1529–1536. pmid:11296141

63. Gordon Huxley, Julian. The variation in isometric tension with sarcomere length in vertebrate muscle fibres. The Journal of physiology 1966, 184:170–192.

64. Moo EK, Leonard TR, Herzog W. The sarcomere force–length relationship in an intact muscle–tendon unit. J.Exp.Biol. 2020, 223:jeb215020. pmid:32098882

65. Sun L, Sun Y, Huang Z, Hou J, Wu J. Improved hill-type musculotendon models with activation-force-length coupling. Technology and Health Care 2018, 26:909–920. pmid:29914041

66. Ebashi S, Endo M. Calcium and muscle contraction. Prog.Biophys.Mol.Biol. 1968, 18:123–183.

67. Farmer TW, Buchthal F, Rosenfalck P. Refractory period of human muscle after the passage of a propagated action potential. Electroencephalogr.Clin.Neurophysiol. 1960, 12:455–466. pmid:13821598

68. Brooks JE, Hongdalarom T. Intracellular electromyography: Resting and action potentials in normal human muscle. Arch.Neurol. 1968, 18:291–300.

69. Ludin HP. Microelectrode study of normal human skeletal muscle. Eur.Neurol. 1969, 2:340–347. pmid:5808475

70. Kopec J, Delbeke J, McComas AJ. Refractory period studies in a human neuromuscular preparation. Journal of Neurology, Neurosurgery & Psychiatry 1978, 41:54–64.

71. Wallinga-De Jonge W, Gielen FL, Wirtz P, De Jong P, Broenink J. The different intracellular action potentials of fast and slow muscle fibres. Electroencephalogr.Clin.Neurophysiol. 1985, 60:539–547. pmid:2408854

72. Wolters H, Wallinga W, Ypey DL, Boom HB. Ionic currents during action potentials in mammalian skeletal muscle fibers analyzed with loose patch clamp. American Journal of Physiology-Cell Physiology 1994, 267:C1699–C1706. pmid:7528975

73. Balog EM, Thompson LV, Fitts RH. Role of sarcolemma action potentials and excitability in muscle fatigue. J.Appl.Physiol. 1994, 76:2157–2162. pmid:8063681

74. Dillmann U, Hopf HC, Lüder G, Schimrigk K. Isometric muscle contractions after double pulse stimulation. comparison of healthy subjects and patients with myotonic dystrophy. Eur.J.Appl.Physiol.Occup.Physiol. 1996, 74:219–226. pmid:8897028

75. Dumitru D, King JC. Motor unit action potential duration and muscle length. Muscle Nerve 1999, 22:1188–1195. pmid:10454713

76. Arabadzhiev TI, Dimitrov GV, Dimitrova NA. Intracellular action potential generation and extinction strongly affect the sensitivity of M-wave characteristic frequencies to changes in the peripheral parameters with muscle fatigue. Journal of Electromyography and kinesiology 2005, 15:159–169. pmid:15664146

77. Kim G, Ferdjallah MM, McKenzie FD. An empirical muscle intracellular action potential model with multiple erlang probability density functions based on a modified newton method. Biomedical Engineering and Computational Biology 2013, 5:BECB. S11646. pmid:25288900

78. Hollingworth S, Zhao M, Baylor SM. The amplitude and time course of the myoplasmic free [Ca2] transient in fast-twitch fibers of mouse muscle. J.Gen.Physiol. 1996, 108:455–469. pmid:8923269

79. Baylor SM, Hollingworth S. Sarcoplasmic reticulum calcium release compared in slow-twitch and fast-twitch fibres of mouse muscle. J.Physiol. 2003, 551:125–138. pmid:12813151

80. Wexler AS, Ding Jun, Binder-Macleod SA. A mathematical model that predicts skeletal muscle force. TBME 1997, 44:337–348. pmid:9125818

81. Hollingworth S, Zeiger U, Baylor SM. Comparison of the myoplasmic calcium transient elicited by an action potential in intact fibres of mdx and normal mice. J.Physiol. 2008, 586:5063–5075. pmid:18772198

82. Zot HG, Hasbun JE. Modeling Ca2 -bound troponin in excitation contraction coupling. Frontiers in physiology 2016, 7:406. pmid:27708586

83. Ding J, Wexler AS, Binder-Macleod SA. Development of a mathematical model that predicts optimal muscle activation patterns by using brief trains. Journal of Applied Physiology 2000, 88:917–925. pmid:10710386

84. Milner-Brown HS, Stein RB, Yemm R. The orderly recruitment of human motor units during voluntary isometric contractions. J. Physiol. 1973, 230:359–370. pmid:4350770

85. Thomas CK, Bigland-Richie B, Johansson RS. Force-frequency relationships of human thenar motor units. J. Neurophysiol. 1991, 65:1509–1516. pmid:1875259

86. Macefield VG, Fuglevand AJ, Bigland-Ritchie B. Contractile properties of single motor units in human toe extensors assessed by intraneural motor axon stimulation. J Neurophysiol 1996, 75:2509–2519. pmid:8793760

87. Fuglevand AJ, Macefield VG, Bigland-Ritchie B. Force-frequency and fatigue properties of motor units in muscles that control digits of the human hand. J Neurophysiol 1999, 81:1718–1729.

88. Negro F, Yavuz U, Farina D. The human motor neuron pools receive a dominant slow-varying common synaptic input. J. Physiol. 2016, 594:5491–5505. pmid:27151459

89. De Luca CJ, Hostage EC. Relationship between firing rate and recruitment threshold of motoneurons in voluntary isometric contractions. J.Neurophysiol. 2010, 104:1034–1046. pmid:20554838

90. Zajac Felix E. Muscle and tendon: Properties, models, scaling, and application to biomechanics and motor control. Critical reviews in biomedical engineering 1989, 17:359–411. pmid:2676342

91. Thompson CK, Negro F, Johnson MD, Holmes MR, McPherson LM, Powers RK, Farina D, Heckman CJ. Robust and accurate decoding of motoneuron behaviour and prediction of the resulting force output. J.Physiol. 2018, 596:2643–2659. pmid:29726002

92. Legreneur P, Morlon B, Van Hoecke J. Simulation of in situ soleus isometric force output as a function of neural excitation. Journal of biomechanics 1996, 29:1455–1462. pmid:8894926

93. Callahan DM, Umberger BR, Kent-Braun JA. A computational model of torque generation: Neural, contractile, metabolic and musculoskeletal components. PloS one 2013, 8:e56013. pmid:23405245

94. Raikova R, Celichowski J, Angelova S, Krutki P. A model of the rat medial gastrocnemius muscle based on inputs to motoneurons and on an algorithm for prediction of the motor unit force. Journal of neurophysiology 2018, 120:1973–1987. pmid:30020845

95. Fuglevand AJ, Winter DA, Patla AE. Models of recruitment and rate coding organization in motor-unit pools. Journal of neurophysiology 1993, 70:2470–2488. pmid:8120594

96. Gogeascoechea A, Ornelas-Kobayashi R, Yavuz US, & Sartori M. Characterization of Motor Unit Firing and Twitch Properties for Decoding Musculoskeletal Force in the Human Ankle Joint In Vivo. IEEE Transactions on Neural Systems and Rehabilitation Engineering 2023. pmid:37756177

97. Elek JM, Kossev A, Dengler R, Schubert M, Wohlfahrt K, Wolf W. Parameters of human motor unit twitches obtained by intramuscular microstimulation. Neuromuscul. Disord. 1992, 2:261–267. pmid:1483052

98. Bigland-Ritchie B, Fuglevand AJ, Thomas CK. Contractile properties of human motor units: Is man a cat? The Neuroscientist 1998, 4:240–249.

99. Raikova R, Krutki P, Celichowski J. Skeletal muscle models composed of motor units: A review. Journal of Electromyography and Kinesiology 2023, 102774. pmid:37099899

100. Hatze H. A general myocybernetic control model of skeletal muscle. Biological Cybernetics 1978, 28 143–157. pmid:630009

101. Perreault EJ, Heckman CJ, Sandercock TG. Hill muscle model errors during movement are greatest within the physiologically relevant range of motor unit firing rates. J.Biomech. 2003, 36:211–218. pmid:12547358

102. Valente G,Pitto L,Testi D,Seth A,Delp SL,Stagni R,Viceconti M,Taddei F. Are subject-specific musculoskeletal models robust to the uncertainties in parameter identification? PloS one 2014, 9, e112625. https://www.ncbi.nlm.nih.gov/pubmed/25390896. pmid:25390896

103. Navacchia A, Myers CA, Rullkoetter PJ, Shelburne KB. Prediction of in vivo knee joint loads using a global probabilistic analysis. J.Biomech.Eng. 2016, 138:031002. pmid:26720096

104. Woittiez RD, Huijing PA, Boom H, Rozendal RH. A three-dimensional muscle model: A quantified relation between form and function of skeletal muscles. J.Morphol. 1984, 182:95–113. pmid:6492171

105. Lubel E, Sgambato BG, Barsakcioglu DY, Ibáñez J, Tang M, Farina D. Kinematics of individual muscle units in natural contractions measured in vivo using ultrafast ultrasound. Journal of Neural Engineering 2022, 19:056005.

106. Herzog W. The problem with skeletal muscle series elasticity. BMC biomedical engineering 2019, 1:28. pmid:32903293

107. Del Vecchio A, Germer CM, Kinfe TM, Nuccio S, Hug F, Eskofier B, Farina D, Enoka RM. The forces generated by agonist muscles during isometric contractions arise from motor unit synergies. Journal of Neuroscience 2023. pmid:36922028

108. Ma S, Guerra IM, Caillet AH, Zhao J, Clarke AK, Maksymenko K, et al. NeuroMotion: Open-source Simulator with Neuromechanical and Deep Network Models to Generate Surface EMG signals during Voluntary Movement. bioRxiv 2023

109. Taylor JL, Amann M, Duchateau J, Meeusen R, Rice CL. Neural contributions to muscle fatigue: From the brain to the muscle and back again. Med.Sci.Sports Exerc. 2016, 48:2294. pmid:27003703

110. Callahan DM, Umberger BR, Kent JA. Mechanisms of in vivo muscle fatigue in humans: investigating age-related fatigue resistance with a computational model. The Journal of physiology 2016, 594(12), 3407–3421. pmid:26824934

111. Shorten PR, O’Callaghan P, Davidson JB, Soboleva TK. A mathematical model of fatigue in skeletal muscle force contraction. Journal of muscle research and cell motility 2007, 28, 293–313. pmid:18080210

112. Oliveira AS, Negro F. Neural control of matched motor units during muscle shortening and lengthening at increasing velocities. J.Appl.Physiol. 2021, 130:1798–1813. pmid:33955258

113. Glaser V, Holobar A. Motor unit identification from high-density surface electromyograms in repeated dynamic muscle contractions. IEEE Transactions on Neural Systems and Rehabilitation Engineering 2018, 27:66–75. pmid:30571641

114. Avrillon S, Hug F, Gibbs C, Farina D. Tutorial on MUedit: An open-source software for identifying and analysing the discharge timing of motor units from electromyographic signals. bioRxiv 2023. 13.548568.

115. Clarke AK, Atashzar SF, Del Vecchio A, Barsakcioglu D, Muceli S, Bentley P, Urh F, Holobar A, Farina D. Deep learning for robust decomposition of high-density surface EMG signals. IEEE Transactions on Biomedical Engineering 2020, 68:526–534.

116. Rossato J, Hug F, Tucker K, Lacourpaille L, Farina D, Avrillon S. I-spin live: An open-source software based on blind-source separation for decoding the activity of spinal alpha motor neurons in real-time. bioRxiv 2023.04. 14.536933.