Musculoskeletal modeling is used to simulate neuromuscular control in various human movements (Delp and Loan, 2000, Dembia et al., 2020, Hatze, 1976; Falisse et al., 2019b; Heinrich et al., 2022, Navacchia et al., 2019, Uchida et al., 2016). This information is useful to understand the neuromuscular control strategy for healthy and pathological movements (Afschrift et al., 2018, Falisse et al., 2020, 2019b; Navacchia et al., 2019, Pitto et al., 2019). Commonly the kinematics and kinetics data for musculoskeletal modeling is obtained with a marker-based three dimensional motion capture system that includes force plates in a laboratory space (Sylvester et al., 2021). These devices are gold standard for the biomechanical measurements (Winter, 2009). For clinical assessment, however, not all the hospitals and clinics can operate the costly motion capture system with force plates. To address this problem, previous studies proposed approaches to predict ground reaction force (GRF) (Dijkstra and Gutierrez-Farewik, 2015, Fluit et al., 2014, Muller et al., 2020, Skals et al., 2017).

A previous study utilized a foot–ground contact model to predict GRF in forward dynamics based musculoskeletal simulations (Falisse et al., 2019b; Koelewijn and van den Bogert, 2020). The contact model applies a force when it penetrates the ground plane. The reported estimated muscle activations were in agreement with measured electromyography (EMG) (Falisse et al., 2019b). It indicated that muscle activations and forces could be successfully estimated without force plates. However, this approach has some difficulties in performing the simulation. It is an optimal control problem that includes a large number of states, control and constraints. Although the performance and robustness has been improved, it is time consuming (36 min for half gait cycle (Falisse et al., 2019b)) and may not converge to find the reasonable solution (Anderson and Pandy, 2001, Dembia et al., 2020). In contrast, inverse dynamics based simulation, especially static optimization, is relatively easy and stable to find optimal results (Erdemir et al., 2007).

Musculoskeletal simulations have been performed with inverse dynamics based approaches in previous studies (De Groote et al., 2016, Mokhtarzadeh et al., 2013, Ueno et al., 2020, Zargham et al., 2019). A challenge in the inverse dynamics based simulation is the errors between measured kinematics and force data (Erdemir et al., 2007), which is called residual forces and moments. These residuals affect inverse dynamics results and subsequently the static optimization results, e.g. muscle activations and forces (Faber et al., 2018). A previous study proposed an approach to predict GRF from motion with contact models and an optimization technique (Muller et al., 2020). In this contact force optimization step, the sum of squared contact forces was minimized while respecting the multibody dynamics. Residuals were mostly eliminated in this step. They reported that the correlation coefficients and root mean square error averaged over the task between L5/S1 joint moment computed with predicted and measured GRF were 0.95 and 14 Nm in flexion/extension direction during a manual material handling task. With its accuracy with eliminated residual forces and moments, it seemed that their approach can derive accurate results from the static optimization in musculoskeletal modeling. However, the validity of muscle activation estimated with predicted GRF has not been reported. For gait trials, the validity of predicted GRF and joint moments using Muller’s approach remain unclear as well. Predicted GRF in double leg stance phase of gait would be less accurate since objective function was simplified from Fluit et al. (2014) which addressed the indeterminacy problem of double leg stance phase. Therefore, the purpose of this study was to determine the validity of predicted GRF and joint moments as well as the muscle activation estimated with predicted GRF in the inverse dynamics based musculoskeletal simulation during gait. It was hypothesized that the predicted GRF and joint moments would be accurate and could derive estimated muscle activation comparable to that with measured GRF.

Methods.

An open-source motion capture dataset that contains gait data from 50 healthy subjects (24 women and 26 men, 37.0 ± 13.6 years, 1.74 ± 0.09 m, 71.0 ± 12.3 kg) was used as the experimental data (Schreiber and Moissenet, 2019). Briefly, participants performed walking at five different speeds: between 0 and 0.4 m.s−1, between 0.4 and 0.8 m.s−1, between 0.8 and 1.2 m.s−1, self-selected spontaneous speed and self-selected fast speed. 52 marker trajectories were recorded with 10 camera optoelectronic system sampled at 100 Hz (OQUS4, Qualisys, Sweden). Two force plates sampled at 1500 Hz (OR6-5, AMTI, USA) were used to record 3D ground reaction forces and moments. Following SENIAM recommendation, a wireless EMG system sampled at 1500 Hz (Desktop DTS, Noraxon, USA) was used to record 8 muscles of the right leg: gluteus maximus, gluteus medius, rectus femoris, vastus medialis, semitendinosus, gastrocnemius medialis, soleus, and tibialis anterior. All these systems were synchronized using the Qualisys Track Manager software (QTM 2.8.1065, Qualisys, Sweden). The data in which the subject failed to step on the force plates with each foot was excluded from the analysis.

C3D files that stored the pre-processed marker trajectories and force plate data distributed with the dataset were used for this study. According to Schreiber and Moissenet (2019), the marker trajectories were smoothed using a fourth order Butterworth low pass filter with a 6 Hz cut-off frequency, while the ground reaction forces and moments were smoothed using a second order Butterworth low pass filter with a 15 Hz cut-off frequency in their study. Raw EMG data that distributed with CSV files were band pass filtered at 10–400 Hz, full-wave rectified and low pass filtered at 6 Hz in this study according to a previous study (Navacchia et al., 2019). Since EMG signals during maximum voluntary contraction were not available from the dataset, EMG was normalized to the peak value of estimated muscle activations as similar to a previous study (Rajagopal et al., 2016). When the EMG was compared to PRED, it was normalized to the peak value of the estimated muscle activation from PRED in the trial, while EMG was normalized to the peak value of the estimated muscle activation from EXP when the EMG was compared to EXP.

A summary of the musculoskeletal simulation pipeline is shown in Fig. 1. CusToM, a MATLAB toolbox was used for the musculoskeletal simulations (Muller et al., 2019). Full body model that consists of 42 degrees of freedom, 82 lower extremity muscles and 28 contact points on the feet was used as a generic model. Subject models were calibrated by optimizing the segment length and marker placements of the generic model (Andersen et al., 2010). Then inverse kinematics, external force prediction, inverse dynamics and static optimization were performed. The external force prediction in CusToM was designed to perform the approach proposed by the previous study (Muller et al., 2020). It first detected the contact points that were in contact state with the ground. When the height and speed of a contact point were below 0.05 m and 0.8 m/s, respectively, the contact point was considered to be in contact with the ground. Then the contact points for each foot in contact state were optimized by minimizing sum of squared forces on the contact points. In the equality constraint of the optimization, the forces and moments around the pelvis that simply calculated from the contact force vectors and pelvis-contact points distance were set equal to the forces and moments at the pelvis calculated from inverse dynamics only with the kinematics without external forces. The maximum force for each contact point was limited to 40 % of the subject's body weight. The predicted GRF was used for the calculation of inverse dynamics. The joint kinematics and GRF data were processed with fourth order Butterworth lowpass filter at 6 Hz applied twice in both the forward and reverse directions before inverse dynamics. Inverse Dynamics derived the joint moments and residual forces and moments. Finally static optimization was performed to simulate muscle activations and forces. The sum of cubed muscle activations was minimized in static optimization as it was the default setting in CusToM. In this study, two sets of inverse dynamics and static optimization were performed, one used predicted GRF (PRED) and another used experimentally measured GRF (EXP).

To evaluate the validity, root mean square error (RMSE) and Pearson's correlation coefficient for all time frames of GRF and joint moments between PRED and EXP were evaluated. Force and moment values were normalized to the subject’s bodyweight (BW) and bodyweight times body height (BW*BH), respectively. Pearson’s correlations between EMG and estimated muscle activations for both PRED and EXP were calculated to evaluate the validity using all time frames as samples for the analysis as similar to a previous study (Zargham et al., 2019). To evaluate the differences between PRED and EXP, paired t-test was used for the trial-wise Pearson’s correlation coefficient of the estimated muscle activation, as well as the root mean square of the residual forces and moments. To investigate how residuals relate to computation of the joint moments, Pearson’s correlation coefficient between residuals and joint moments using all time frames were computed. The values for residuals and joint moments were calculated as the absolute value of the PRED-EXP difference. Stance phase for GRF, full right gait cycle for joint moments and muscle activation were included for statistics. With regard to the residuals, the loading response phase (first double leg stance phase) lacked the GRF measurement on the other foot and was excluded from the statistical analysis. Note that the second double leg stance phase (pre-swing phase) was included in the analysis. The figures show full gait cycle.