Solving musculoskeletal biomechanics with machine learning

Musculoskeletal polynomial model for generation of training and testing datasets

We have previously developed the method of autogenerated polynomial models (Sobinov et al., 2020). In these polynomials, the composition of terms was expanded using objective information measurements, i.e., the corrected Akaike Information Criterion. In brief, the posture-dependent musculotendon actuator length and joint moment arms for each muscle in the upper-limb model (Gritsenko et al., 2016) were accurately approximated using the selection of up to 5th power polynomial terms, where muscle length and moment arms were connected through a partial derivative of the muscle length in local coordinates corresponding to limb posture. Overall, the 18 DOF model of the human arm and hand is actuated by 33 muscles, each spanning about 3 DOFs and up to 6 DOFs for thumb muscles. Thus, each actuator is represented by a set of one length- and about 3 moment arm-posture polynomials. The goal of this development was to bypass costly calculations of geometrical transformations with high-quality approximations. Previously, we have demonstrated high-fidelity of these approximations with kinematic errors below 1% (Sobinov et al., 2020). In biomechanics, the errors of 1–5° in joint angles are expected from flaws in the observations in motion capture, and errors of 2° and less are not meaningful in the clinical context (McGinley et al., 2009). Thus, the errors below 1% of joint range of motion are negligible. These polynomial models of muscle posture-dependent state were used to develop an ANN-based approximation method for the musculoskeletal dynamics in this study.

Training, validation and testing datasets

Training, validation and testing datasets for the assessment of model performance were generated by the musculoskeletal polynomial model, which was used as a reference in this study. Input-output relationships were extracted randomly with uniform distribution where inputs were 18 DOF vectors of joint angles and the outputs were 33 length vectors and 99 moment arm vectors. An average muscle crosses 3 DOFs and has, consequently, 3 moment arm relationships on average. We used the supervised learning approach for training ML models. The training dataset was used for two tasks, tuning the model hyper-parameters and model training, to maximize the model performance in replicating the desired outputs with given inputs. The testing dataset contained ~5% of all data (5 × 10⁴ samples). The remaining ~95% were divided into the training dataset (80%, 8 × 10⁵ samples) and the validation dataset (20%, 2 × 10⁵ samples). The validation dataset was used to prevent overfitting, i.e., higher performance on the training data as compared to that on the validation data. These datasets were similarly used for training ANN and LGB models, described below. Overall, the training time was about 15 times longer for ANN then for LGB models. The training of all ANN and LGB models on the standard hardware took about 3.5 days.

Metrics

The performance of the trained models was further evaluated with the testing dataset, which was not used during the training procedure. We expected to reach the same error tolerances as in our previous polynomial fitting method study (Sobinov et al., 2020). Consequently, we used the same normalization of lengths and moment arms as in our previous work. The RMSE values were calculated as the absolute difference between reference and predicted muscle length values. To normalize the results, we divided each reference and predicted length value by the muscle length range respectively:

$$RMS{E}_L={\displaystyle \frac{1}{m}\sum _{l=1}^m{\displaystyle \frac{1}{n}\sum _{i=1}^n\sqrt{{\left({\displaystyle \frac{x_{r,l,i}-x_{p,l,i}}{Lmax{}_l-Lmin{}_l}}\right)}^2}}}$$

where m is the number of muscles (m = 33), n is the number of test samples, $x_{r,l,i}$ and $x_{p,l,i}$ are reference and predicted length values, respectively, $Lma{x}_l$ and $Lmi{n}_l$ are the maximum and minimum values over the full range of lth muscle length.

Similarly, the RMSE of moment arms was calculated as the absolute difference between reference and predicted values, which were normalized to the moment arm maximum ( $Mma{x}_l$):

$$RMS{E}_{MA}={\displaystyle \frac{1}{m}\sum _{l=1}^m{\displaystyle \frac{1}{n}\sum _{i=1}^n\sqrt{{\left({\displaystyle \frac{x_{r,l,i}-x_{p,l,i}}{Mma{x}_l}}\right)}^2}}}$$

where $x_{r,l,i}$ and $x_{p,l,i}$ are reference and predicted values, m is the number of moment arms (m = 99), and n is the number of test samples.

Machine learning models

Two types of ML models were used to map the musculoskeletal input-output relationships. We used the light gradient boosting machine (LGB) models and artificial neural network (ANN) with two hidden layers. The models were trained and tested according to the workflow in Fig. 2A. Validation accompanied the training process to prevent overfitting.

Estimation of the training dataset size

The selection of the training dataset size for ANN and LGB models is a non-trivial step in the model development. Our source of data was expressed functionally allowing unlimited source of training data. However, the selection of an optimal dataset that captures the relationships without the tendency for overfitting was the initial goal of our development. We used RMSE metric for both length and moment arm models trained with several datasets of incremental size. The relationship between the metric and the dataset size are shown in Fig. 3. As the size of the dataset increased logarithmically (from 10³ to 10⁶ samples), the training accuracy also increased, with minor improvement in the range above 10⁵ samples. The improvements with the dataset size were not as pronounced showing 1.45% and 1.94% errors with the smallest datasets (10³ samples). The improvement curve of LGB is flat, showing no further improvement, after 10⁵ sample size. The performance of the relatively simple (Biceps Brachii Long Head) and complex (Extensor Pollicis Longus) muscles is illustrated in Figs. 3C and 3D. These two muscles are on the opposite extremes of complexity expressed as the number of terms required for accurate polynomial fit (Sobinov et al., 2020). The improvements are qualitatively similar for both of these muscles indicating no strong dependency on muscle path complexity with both types of models. We used the largest dataset size for all further model development described below.

Model accuracy

High accuracy was achieved with both LGB and ANN model types. The distribution of errors is shown in Fig. 4 with the histograms of RMSE values for the testing dataset (5 × 10⁴ samples). The highest achieved performance of ANN models was with 10⁷ with absolute errors at 0.08 ± 0.05% for muscle lengths and 0.53 ± 0.29% for moment arms. Similarly, LGB models generated accurate predictions with large training datasets, 0.18 ± 0.06% and 0.13 ± 0.07% errors, respectively. This is the expected error rate based on our previous analysis (Sobinov et al., 2020).

Overall, the error span did not exceed 0.6% for muscle lengths and 2% for muscle moment arms, shown in Figs. 4 and 5. We found no relationship between length and moment arm errors (p = 0.746, R² = 0.003). The level of errors was comparable for both simple and complex muscles (spanning more than 3 DOF), e.g., the errors of ECR_BR (a two DOF muscle) were comparable to those of EDM. The accuracy of LGB and ANN models was comparable. The interquartile ranges (IQR), corresponding to the distance between 25% and 75% level for the distribution of all length error values were 0.075% (ANN) and 0.216% (LGB). The 25-75% IQRs for moment arm errors in Fig. 4B were 0.464% (ANN) and 0.0782% (LGB). Median RMSE values for all models were less than 0.01%. To check if accuracy declines in the extreme postures (1st and 4th quartiles), we repeated testing and found a similar rate of errors, 0.144% and 0.587% for lengths and moment arms.

The distribution of absolute errors is shown for each muscle in Fig. 5. The muscles are sorted according to the number of DOFs they actuate with relatively simple muscles on the left and complex (thumb) muscles on the right. Overall, the majority of distributions were below 0.2% for 75% of all vales, as indicated by the top value of the interquartile range in the box plots. The largest length errors were observed in OP, which was also one of the most difficult muscles to design structurally (Boots et al., 2020). In this muscle the top value (Q3) of the interquartile range was about 0.3%, which corresponds to the error of 0.018 mm. In Fig. 5A, the normalized errors of muscle length are presented for each muscle. For all muscles, the most errors (up to 75% of the distributions) are below 0.2%. The prominent exception is OP with the highest normalized errors, which is explained by the minimal full physiological range of only 6 mm. The error of 0.3% in OP length corresponds to the absolute error of 0.018 mm. The errors did not increase with muscle structural complexity. The evaluation errors in muscle length where generally larger in the group of muscles spanning 2 DOF (blue, see Extensor Carpi Radialis, ERC_LO) and were comparable to the errors in complex muscles (e.g., Abductor Pollicis Brevis, APB).

Training and evaluation time

The execution times were compared for ANN and LGB models (1.4 GHz Quad-Core 8th-generation Intel Core i5) by measuring the duration of 1000 evaluations (using method time from the standard time library in Python 3.7). For a given posture, ANN models evaluated both muscle length and moment arms with the combined latency of 1.1 ± 0.6 ms, as compared to 43.1 ± 8.3 ms for LGB models, which were about 39 times slower.

Motor intent decoding

Estimating limb posture from EMG in real-time applications remains to be a challenge in human-machine interfaces due to: (1) the difficulty in the theoretical description and (2) the lack of experimental data to validate these models (Alber et al., 2019). In general, a statistical mapping between posture and recorded activity from descending pathways, nerves, and muscles has been used as the transformation to predict motor intent (Wang & Buchanan, 2002; Ting et al., 2019; Dantas et al., 2019), to investigate interplay of mechanical and neural components in pathologies (Sartori et al., 2017) or to control powered prosthetic limbs or exoskeletal devices (Kiguchi & Hayashi, 2012; Collinger et al., 2013; Zhang et al., 2017; Furui et al., 2019; George et al., 2020). However, the accuracy of decoding realistic movements remains to be a challenge especially for movements that require dexterous object manipulation (Downey et al., 2017). Many current decoding methods in brain-computer interfaces assume that neural activity is related to limb end-point position and/or velocity and lack the description of movement kinetics generated by muscle forces. The resulting movements of prosthetics are typically slower and less robust than natural movements. Up to five (Wendelken et al., 2017) to six (George et al., 2020) arm and hand DOFs can be simultaneously controlled using nerve signals recorded with penetrating electrodes. Accurate but slow movements can be generated for high-dimensional artificial hands with wrist and digits; however, the accuracy is challenged by changing limb posture and orientation. One potential solution is the use of closed-loop control systems that provide not only the forward control of prosthetic, but also incorporate the sensory feedback within neuroprosthetics (Hughes et al., 2020; Ganzer et al., 2020; Charkhkar, Christie & Triolo, 2020). The closed-loop control system that takes into account muscle forces would critically depend on the accurate representation of musculoskeletal actions that we describe with our novel approach.

The generalizable control solutions based on the biomechanical transformations have potential advantages over the nontransparent statistical approaches. Mechanistic musculoskeletal models of legs (Sartori et al., 2017) and arms (Crouch & Huang, 2016; Boots et al., 2017; Mansouri et al., 2017; Sartori et al., 2018) have been previously used for motor intent decoding. The major advantage of biomechanical modeling over statistical methods is in the explicit representation of kinematic and kinetic dependencies in the generated motor command signals during multisegmented limb movements. Thus, the transformation from the recorded biological signals to the proportional control of limbs should be intuitive, given that the underlying computations are sufficiently accurate and without extensive processing delays.

Together, these possible solutions motivated our exploration and rationale for developing accurate ML methods of approximating the musculoskeletal transformations with physics-informed neural networks. We used model driven training and testing of ML algorithms to approximate posture-dependent changes in muscle lengths and moment arms of distal arm and hand muscles. The previously developed polynomial model provided a functional representation of data across all possible limb postures. Since any volume of data could be generated, we tested the extent of data required to train ANN and LGB models. Figure 3 shows the expected inverse relationship between the approximation errors and the training dataset size. Using only 10⁶ samples for training ANN and LGB models resulted in kinematic errors that were less than 0.5%. This number of samples is about one order of magnitude lower than the number of samples required for our previous implementation using an information theory-based algorithm for approximating the same kinematic variables. This may indicate that ANN and LGB methods may further generalize the polynomial relationships that required a larger volume of samples to generate functional muscle-posture representations limited to polynomial power terms.

Computational delays

The transformation of EMG signals into movement rarely accounts for the musculoskeletal anatomy and physiology. This is partly due to the extreme complexity of muscle organization (Gritsenko et al., 2016; Murphy et al., 2018) and nonlinear intrinsic muscle properties that include independent force-length and force-velocity relationships (Zajac, 1989) and less popular in modeling, short-range stiffness, which is a hysteretic force-length property (Cui et al., 2008). The task of simulating muscle force generation requires adequate structural information about muscle paths and posture-dependent changes in moment arms. The development of complex musculoskeletal models has been recently simplified by the dedicated simulation tools for editing and simulating segmental dynamics—OpenSim (Seth et al., 2018), MuJoCo (Todorov, Erez & Tassa, 2012), Simscape (MathWorks, Inc.). The challenge remains in collating sufficient datasets of musculoskeletal measurements for creating complex musculoskeletal models and then in testing and validating these models across the full-range of motion to ensure their use in a wide range of applications.

The computational delays of solving the equations of motion governing limb dynamics have previously impeded the application of model-based prosthetic controllers. In particular, the evaluation of muscle force-velocity characteristic may require sub-millisecond latencies to decode accurately rapid movements computed by a physical engine, which requires additional time to execute (~1 ms in Mansouri et al., 2017). Our implementations demonstrated a clear speed advantage of the ANN over the LGB model (about 39× faster) and required about 1 ms on standard hardware. This performance was still about 15 times slower than our polynomial approximation (about 60 µs); yet, it provides a milestone for the future improvements due to the rapid development of both software and hardware solutions for ANNs. Further improvements in the performance of the ANN may be possible approaching the latencies appropriate not only for the feedforward computations, but also for predictive inverse computations inspired by theoretical neural transformations (Wolpert & Ghahramani, 2000). Another biomimetic feature of the ANN model is its potential solution for the necessity to increase computational complexity to accommodate the increase in the size of described structure, typically termed as “the curse of dimensionality”. Sobinov et al. (2020) demonstrated that the typical exponential increase could be replaced with the linear increase in the required number of terms within the polynomial approximations. Here, the same structure of the ANN model accommodated accurate calculations for a subset and for the full set of 33 muscles. We hypothesize that as long as the additional simulated muscles are relatively similar in anatomical complexity to the muscles represented in our current model, the same structure of the ANN will be able to embed their dynamics without the increase in the number of nodes in each layer.

Limitations

Our approach described the kinematic transformations only and improves the calculation of instantaneous muscle forces, yet, it lacks the description of equations of motion. The full transformation to the intended joint and segment kinematics requires the forward simulation of segmental dynamics by a physics engine, for example, MuJoCo (Todorov, Erez & Tassa, 2012). Another minor limitation is the use of prior simplified model of muscle kinematics to supervise the training in our ML models. Small errors of about 1% are expected to propagate to this implementation. Given that models of this complexity often include hundreds of poorly validated or guessed parameters, the expected performance is within the “good enough” qualitative range for this model (McGinley et al., 2009). For example, in our development of biomechanical models with validated moment arm kinematics, we have discovered large discrepancies between models and experimental data (Boots et al., 2020). It is also expected that the individual morphological differences may exceed the numerical errors in our model. In the future work, we will examine the use of ML models in the description of subject-specific musculoskeletal transformations, where the current generic model can be used as the starting pre-trained model. Another current limitation is the demanding training dataset and also the slow execution of the ANN implementation, which is appropriate for real-time applications but is still slower than our previous polynomial approximation.