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

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.


Methods
The force   (in Volts) recorded by the load cell during the experiments was first converted into an ankle torque  (in N.m) using a constant gain derived from the load cell documentation and the geometry of the dynamometer.It was then assessed with EQ.S1, for the recordings made during the second session where bEMG electrodes were involved, the amount ∆() of ankle torque taken over the trapezoidal contraction by agonist EHL, EDL and antagonist SOL, GM, and GL muscles.In EQ.S1,   () is the ankle torque produced by the i th muscle,   ̅̅̅̅̅̅̅ () is the normalized EMG envelope for muscle i computed previously,  is the moment arm between the tendon of muscle i and the ankle joint previously measured from the subject-specific MSK model,  0,  is the maximum isometric force of muscle i estimated previously, and   is the scaling factor calculated with the FL relationship in EQ.11.
̅̅̅ is the normalized length of muscle i and was calculated from the subject-muscle-specific    ,  ,  , and  0,  values derived previously, assuming the tendons to be rigid.It was also assumed that the EDL and EHL muscles produced the same electrical activity.
From EQ. S1, the continuous relationship between the level of load sharing Δ() and the recorded ankle torque  was computed by fitting the exponential function in EQ.S2 to the (; Δ) cloud of points.

Δ𝑇(𝑇) = 𝑎 • 𝑒 𝑏•𝑇
EQ. S2 The fitted Δ() relationship in EQ.S2 was then applied to the torque () recorded during the first experimental session, for which the level of co-contraction Δ(()) was estimated.The experimental TA muscle force   () produced during the first experimental session was finally estimated with EQ.S3, where   is the moment arm the TA tendon makes with the ankle joint.

Results
During the second experimental session, voluntary ankle torques  (black dashed traces in Fig S1) of 5.2, 8.8, and 17.5 Nm were recorded experimentally over the dorsiflexion tasks up to 30%, 50% and 100% MVC, respectively.The results are of the same order of magnitude as the maximum voluntary torque of 22.5 Nm (+28%) calculated with published age-sex-dependent joint angle-torque relationships for ankle dorsiflexion [1].Using the subject-muscle-specific properties of maximum isometric forces  0  and moment arms  reported in Table 2, a maximum torque in ankle dorsiflexion of 41 Nm was theoretically reachable for the participant if the TA, EDL, and EHL muscle produced their maximum force (maximum myoelectric activity at optimal length, i.e.,  ̅̅̅̅̅̅̅ = 1 and  ̅ = 1 in EQ.S1) and no co-contraction of the SOL, GM, and GL muscles occurred.In the conditions of the experiment, the TA, EDL, EHL, SOL, GM, and GL worked at 1.18, 0.78, 1.20, 0.60, 0.68, and 0.80 normalized lengths assuming rigid tendons, respectively (Table 2).The myoelectric activity of the EDL/EHL, SOL, GM, and GL muscles reached 27%, 5%, 3%, 6%, and 50%, 16%, 8%, 16%, and 100%, 37%, 30%, 47% of their maximum recorded EMG amplitude, at 30% MVC, 50% MVC, and 100% MVC respectively.
Very similar relationships Δ(), calculated with EQ.S1 and EQ.S2 with the subject-muscle-specific properties reported in From EQ. S3, the force developed by the TA was estimated to reach 260 N and 656 N (torques displayed as solid red lines in Fig S1) at 30% and 50% MVC respectively, i.e. 25% and 63% of the TA's maximum isometric force  0  (Table 2), which seems reasonable for a participant who aimed to isolate the TA during those contractions.At 30% MVC and 50% MVC, the agonist EDL and EHL muscles collectively produced around 112 N and 215 N of force respectively, i.e., less than half the amount of force produced by the TA muscle despite having more than half the force-generating capacities, which is again consistent with a participant who aimed to isolate the TA during those contractions.and EHL and antagonist SOL, GM, and GL muscles (blue dotted traces) was calculated with EQ.S1 using the subject-muscle-specific properties reported in Table 2. Δ was negative during the contractions, i.e., the agonist and antagonist muscle altogether produced a torque in plantarflexion.The ankle torque developed by the TA muscle was calculated as  − Δ (green dash-dotted traces) and as  − Δ() (solid red traces), where Δ() given in EQ.S4 is the continuous relationship fitted to the (; Δ) cloud of points.

Physiological and mathematical steps taken to simplify the neuromuscular model
The +   =   EQ.S5 This approach neglects both the geometrical orientation of the individual MUs and the lateral force transmission between MUs, both of which would decrease the total force developed by the MU pool along the tendon if considered.
In the following, all TA quantities are normalized with respect to the MU or whole muscle optimal length or maximum isometric force and are reported with a bar following EQ.S6.In EQ.S6,  ,0  ,  ,0  ,    , and   ̅̅̅ are the maximum isometric force, the optimal length, the length, and the contraction velocity of the FG of MU , respectively.EQ.S8

𝐹
EQ. S8 cannot be solved without knowing, for each MU, the value of its individual fibre architecture (optimal fibre length, fibre length), which implies knowing their spatial distribution, at all frames and active state.The muscle quantities are currently unfeasible to measure and therefore we implemented a number of simplifications displayed in Fig S2 and detailed in the following.First, it was assumed that all the MUs in the MU pool shared the same optimal MU length  0  , which was taken to be the subjectspecific muscle optimal length  0  derived previously in the main manuscript in EQ. 2. Second, expecting the MU optimal lengths  0,  to be related to the MU lengths    , all the MUs in the MU pool were and rearrange EQ.S8 as a differential equation of the common normalized MU length  ̅ , the MU FV relationship   was simplified to be independent from the individual MU active state   , which yielded EQ.S9.
Although EQ.S9 can be solved numerically for  ̅ , the computational cost and the complexity of the model were further reduced by first assuming the TA tendon to be rigid {4}.As     0  = 3.5 for the participant's TA muscle, this assumption was judged acceptable for low-to-medium TA isometric contractions, considering that     0  = 1 describes a 'very stiff' tendon in the literature [3].With this assumption, the tendon length is constant and equals the tendon slack length [4]   =    .As   is also constant during isometric contractions, so is the common normalized MU length  ̅ , in which case   = 1.Second, because   =    + 1.16 •  0  for the segmented TA muscle in the conditions of the experiment (see Results Section for the subject-specific values calculated with EQ. 2), it was finally approximated {5} that the FGs worked at the common normalized length  ̅ = 1.16, in which case   ̅̅̅̅̅̅ = 0.02 ≈ 0, according to EQ. S10 with  1 = 5 and  2 = 0.6 for human young adults [5].
With these five key simplifications {1-5}, the PEE and SEE (in grey in Fig 1F ) were neglected and the rheological structure was reduced to the cohort of in-parallel MU actuators (in black in Fig 1F ), the dynamics of which reduced to EQ. S10 in the main text of the manuscript Step-by-step simplification of the MN-driven model.After neglecting the pennation angle between the TA muscle and tendon, the MUs are first assumed to share a common optimal length, which is the whole muscle's optimal length  0  {1} and a common normalized length  ̅ {2}.Then, the MUs' FV relationships are simplified to be activation-dependent {3}, before the tendon is assumed rigid {4}.Consequently, the SEE is neglected, the muscle contracts isometrically, and the FV relationship is neglected.Finally, based on subjectspecific MSK measurements, it yields  ̅ = 1.16, in which case the PEE develops negligible force and is neglected {5}.

Derivation of the free Calcium dynamics
Despite some disagreements [6], it has been repeatedly observed in in vitro experiments on amphibian fibres from the semitendinosus, iliofibularis and TA muscles [7][8][9][10] that the peak value of the calcium concentration nonlineary varies with the sarcomere length   .The experimental data provided by these studies was manually digitized and normalized to the highest peak [ 2+ ] value retrieved in each study and to the optimal sarcomere length  0  ( ̅ =    0  ), which is typically 2.1 in these amphibian muscles.The piecewise linear trendline  1 ( ̅ ) given in EQ. 15 was fitted to the normalized data, which reveals that the peak [ 2+ ] value is obtained between 1.15 and 1.30  0  .
The  2 ( ̅ ) scaling factor described in EQ. 16 was obtained following a similar approach and models the linear increase in   with longer normalized fibre lengths  ̅ above 1.15 0  observed experimentally [10].
In the literature of Hill-type muscle models, the coefficients scaling the mathematical equations describing the calcium transients ( 1 ,  2 ,  3 in EQ. 14) are typically tuned to match dated experimental observations of calcium waves obtained in vitro from amphibian muscle fibres in low-temperature baths [11].Yet, the calcium dynamics in the sarcoplasm are strongly species-dependent, which transients twice longer in amphibians than in mammals [12], and temperature-dependent [7,10,[13][14][15][16][17][18][19][20][21][22][23].Moreover, these dated experimental observations were typically obtained with chemical indicators (murexide, aequorin and arsenazo III for example) that have a high affinity to components other than calcium ions in the myoplasm, such as Mg ions, and report calcium twitches of too low peak values and delayed transients.Because of these limitations, rather than considering the values proposed in Hatze's work [24] for the  1 ,  2 ,  3 coefficients, which were obtained in frog fibres at 9°C with the aequorin indicator, the  1 ,  2 ,  3 coefficients were here tuned to match the experimental data from two experimental studies [12,25], where the experimental [ 2+ ] waves were measured at 28-35°C in rodent fibres with the reliable [10,12,15,19,21,25] mag-fura-2 indicator.Lacking human data in the literature, it was believed that the calcium dynamics in rodents at 35°C were the most representative for the human species at body temperature.Because both studies [12,25] performed their experiments at the same sarcomere length   = 3.8 ≈ 1.6 0  [26], the length-dependent scaling factors were set to  1 (1.6) = 0.82 and  2 (1.6) = 0.82 when tuning the  1 ,  2 ,  3 coefficients.

Fig S1 :
Fig S1: Estimation of the ankle torque developed by the TA muscle (solid red traces) from the recorded ankle torque  (dashed black traces) during the trapezoidal contractions up to 30% (Left) and 50% MVC (right)performed during the second experimental session.The amount of ankle torque Δ produced by the agonist EDL and EHL and antagonist SOL, GM, and GL muscles (blue dotted traces) was calculated with EQ.S1 using the subject-muscle-specific properties reported in Table2.Δ was negative during the contractions, i.e., the agonist and antagonist muscle altogether produced a torque in plantarflexion.The ankle torque developed by the TA

Fig S2:
Fig S2:Step-by-step simplification of the MN-driven model.After neglecting the pennation angle between the TA muscle and tendon, the MUs are first assumed to share a common optimal length, which is the whole muscle's optimal length  0  {1} and a common normalized length  ̅ {2}.Then, the MUs' FV relationships are simplified to be activation-dependent {3}, before the tendon is assumed rigid {4}.Consequently, the SEE is neglected, the muscle contracts isometrically, and the FV relationship is neglected.Finally, based on subjectspecific MSK measurements, it yields  ̅ = 1.16, in which case the PEE develops negligible force and is neglected {5}.

Table 2 ,
were obtained for both trapezoidal contractions up to 30% and 50% MVC between the amplitude of measured ankle torque () (black dashed traces in Fig S1)and the amount of ankle torque Δ taken by the group of agnostic EHL and EDL and antagonist SOL, GM, and GL muscles (blue dotted traces in FigS1).The average of the two relationships yielded the relationship in EQ.S4 (95% confidence ranges: [-0.1235 ; -0.1170] and [0.3428 ; 0.3479], r 2 = 0.84).EQ.S4 suggests that contractile activities of the antagonist muscles SOL, GM, and GL overtook the ankle torque produced by the agonist EDL and EHL muscles following an exponential tendency with increasing recorded torques  with, for example, Δ = −1.6  at  = 5.2  (at 30% MVC) and Δ = −7.7  at  = 8.8  (at 50% MVC).
[2]driven neuromuscular model developed in this study(Fig 1F)consists in a population of  inparallel FGs, that describes the force-generating activity of the MU pool according to the cascading neuromuscular dynamics of the individual MUs in Fig 3.The FGs are placed in-parallel with a single muscle-scale passive elastic element (PEE) and in-series with a common in-series elastic element (SEE), i.e., the tendon.isthe length of the vector of available input spike trains controlling the model and takes the values   or  in this study.The FGs individually develop the MU active forces    ,  ∈ ⟦1; ⟧, which linearly sum with the PEE force   to yield the total muscle force   .Neglecting the low pennation angle between the muscle belly and the tendon of the TA (approximately 11°[2]), which may underestimate   by approximately 2%,   equals the tendon force   by equilibrium in EQ.S5.