An Electromyography-Assisted Musculoskeletal Simulation With Concurrent Optimization of Muscle Excitations and Knee Joint Kinematics Open Access

Numerous musculoskeletal (MSK) models of the knee have been developed and vary from simple planar to complex multidegree-of-freedom (DoF) representations of the joint. These models aim to discern knee biomechanics, such as knee kinematics and joint contact force (JCF) [1–4], as these mechanical quantities are difficult to measure noninvasively [5]. The MSK models have been used in different fields, including the prediction, prevention, and treatment of knee MSK disorders, such as knee osteoarthritis prediction and evaluation of surgical treatments. Hence, accurate estimates of knee biomechanics are essential to understanding the etiology, consequences, and treatments of neuromusculoskeletal disorders [6]. Nevertheless, MSK model estimates are sensitive to modeling features, such as the neural solution employed [7,8], knee model's DoF [1,8–11], and knee contact geometry and contact formulation [12–14]. Thus, it is important to understand which MSK model is required for valid and reliable estimation of the quantities of interest.

Within a majority of MSK models, the knee is modeled as a kinematically driven joint in which the estimated knee kinematics are independent of the knee kinetics, e.g., muscle forces. Among these, rigid hinge-type (1DoF) knee joints are more conventional [3,13,15], in which patellofemoral and tibiofemoral secondary kinematics (e.g., translations and internal/external and abduction/adduction rotations) are constrained as functions of knee flexion angle. These constraints are typically formulated based on passive flexion of cadaveric knees [3,15]. More complex 1DoF closed-kinematic-chain multilink knee mechanisms have also been developed with subject-specific articulating surfaces and ligaments, demonstrating great potential for accurate estimation of secondary knee kinematics [16,17]. However, these knee models are kinematically driven [16,17], and hence, the estimated kinematics are independent of knee kinetics, and some such models, when loaded, do not well estimate subject-specific secondary knee kinematics, especially the translational DoF [1,9,10,18]. Accurate estimation of secondary knee kinematics may be crucial, for instance, in treating MSK knee disorders, which may develop due to the altered knee kinematics, even in reduced knee JCFs [19]. Nevertheless, a 1DoF representation of the knee may be able to replicate in vivo knee JCF in a limited number of activities, namely, walking [1,10]. But, these simplified models often fail to estimate plausible muscle excitations and JCF in activities other than walking [11,18] or analysis of individuals with neuromusculoskeletal disorders [20–23]. Also, to our knowledge, the more complex multilink mechanisms knee models have, thus far, not been incorporated into the MSK models [16,17].

Musculoskeletal analyses often use optimization to solve muscle redundancy and estimate muscle forces. Most of these algorithms rely solely on the MSK model's state variables, e.g., position, velocity, and acceleration of the model's DoF, and mass and inertial system dynamics. Static optimization (SO) [24–26], computed muscle control [27], and force-dependent kinematics [28] can be categorized as such methods. However, these methods do not consider individualized muscle excitation patterns. Indeed, muscle excitation patterns are known to be abnormal in individuals with neuromusculoskeletal disorders (e.g., knee osteoarthritis) [22,29], and are likewise altered by the specifics of the task demands (e.g., position versus force control [30]), pain [31], and due to training [32].

Electromyography (EMG)-assisted MSK models with calibrated parameters for muscle–tendon contraction and muscle activation dynamics have been developed to account for subject-dependent variations in these parameters [33–35], as well as subject- and task-specific muscle recruitment strategies [2,33]. The calibration process refines the MSK model's muscle–tendon parameters, which often vary across subjects, to minimize the mismatch between predicted and experimental joint moments and muscle excitations [2,34]. Once calibrated, EMG-assisted MSK models have a demonstrably more accurate estimation of muscle excitations and JCF compared to SO-based and computed muscle control neural solutions [2,34,36–38]. Despite encouraging developments, currently available EMG-assisted MSK models have limited DoF knee joint models (e.g., 1DoF) and are thus subjected to the aforementioned limitations.

To overcome the drawbacks of kinematically driven knee joint models, multi-DoF (up to 24DoF) knee models that also respond to load have been developed for concurrent optimization of secondary knee kinematics and muscle excitations [1,25,26,28,39,40]. Although these models have shown excellent agreement with experimental measures [1,26,39], they are limited to SO-based neural solutions [1,25,26,28,39–41], thus potentially incapable of generating joint kinematics and kinetics to subject- and task-specific muscle excitation strategies [2,8]. Further, available MSK models with multi-DoF knee joints omit muscle contraction and activation dynamics from their analysis [1,25,26,39]. Efforts have also been made to develop multi-DoF EMG-assisted knee models [8,42], but these models are computationally demanding (up to several days), and the knee kinetics are estimated using a simplified 1DoF knee model that is decoupled from the effects of secondary knee kinematics. While continuous development in computational resources makes these heavy simulation pipeline more accessible, certain applications, such as gait retraining that involves (real-time) biofeedback, rely on rapid (milliseconds to minutes) simulation methods.

In this study, we developed and validated an EMG-assisted MSK modeling method that explicitly incorporated a 12DoF knee joint model with nonrigid subject-specific knee joint, subject- and task-specific muscle excitation patterns, muscle–tendon calibration, and concurrent optimization of muscle excitations and secondary knee kinematics. Secondarily, we performed similar MSK analyses but without calibration and EMG assistance, as well as using conventional MSK models with a 1DoF knee joint with and without EMG assistance, to determine whether the new modeling method outperformed other previously published methods. Model predictions were compared against measured EMG, joint moments from inverse dynamics, and in vivo tibiofemoral JCF from two grand challenge datasets [5].

The calibrated EMG-assisted MSK model with the 12DoF knee joint was created, and its performance was compared to four different models, with all five models using the same experimental data. The model with the 12DoF knee had a 6DoF tibiofemoral and a 6DoF patellofemoral joints constrained by knee ligaments. Knee contact was an elastic foundation model with contact geometries from subjects' knee implant.

The four comparison MSK models were (i) the EMG-assisted model with 12DoF knee joint but with uncalibrated muscle–tendon parameters, (ii) the SO-based model with 12DoF knee joint and uncalibrated muscle–tendon parameters, (iii) an EMG-assisted model with 1DoF knee joint and calibrated muscle–tendon parameters, and (iv) an SO-based model with 1DoF knee joint. The EMG-assisted and SO-based models with 12DoF or 1DoF knee models were identical except for the different muscle–tendon parameters and neural solutions. We also performed two additional analyses to further investigate the effect of muscle–tendon calibration and EMG assistance, which are presented in the Supplemental Materials on the ASME Digital Collection for interested readers (Sec. S1.5 available in the Supplemental Materials).

We used the experimental data from the fifth and sixth grand challenge datasets [5]. Two subjects were analyzed: (1) subject code PS with left knee implant (male, mass = 75 kg, height = 180 cm); and (2) subject code DM with right knee implant (male, mass = 70 kg, height = 172 cm) [5]. The experimental data consisted of motion capture data (i.e., marker trajectories and ground reaction forces), synchronized with EMGs from 15 muscles of the studied leg, and in vivo tibiofemoral JCF measured through the implant [5]. EMGs were measured from medial and lateral gastrocnemius, soleus, rectus femoris, vastus lateralis, vastus medialis, long head of biceps femoris, semimembranosus, peroneal, adductor magnus, gluteus maximus, gluteus medius, tensor fascia latae, tibialis anterior, and sartorius. To generate EMG envelopes, the recorded EMG signals were bandpass filtered (30–300 Hz), notch filtered (50 Hz, to remove powerline noise), full-wave rectified, zero-lag low-pass filtered (6 Hz), and then normalized to the single greatest value of similarly processed EMG data from maximum isometric voluntary exertion trials and the successful trials the subject undertook [43]. We used implant geometries provided in the dataset for modeling [5,44].

For subject PS, we analyzed 29 successful trials of the stance phase from overground normal gait (six trials), treadmill normal gait (20 trials), and right turn gait (three trials). For subject DM, we analyzed 98 successful trials of the stance phase from overground normal gait (six trials), overground crouch gait (five trials), overground bouncy gait (six trials), overground smooth gait (two trials), treadmill normal gait (14 trials), treadmill slow-speed gait (six trials), treadmill fast-speed gait (19 trials), treadmill multiple-speeds gait (29 trials), and treadmill crouch gait (ten trials). Preprocessing of data is explained in the Supplemental Materials (Sec. S1.1 available in the Supplemental Materials).

Importantly, the in vivo JCF was not used for MSK analyses; it was only used for comparison against the results.

The 12DoF MSK models (Fig. 1) had three computational modules: (1) creating the knee joint model, (2) calibrating the EMG-assisted MSK model's muscle–tendon parameters, and (3) executing the EMG-assisted and SO-based MSK analyses. These steps are detailed in Secs. 2.2.1–2.2.3.

The 12DoF EMG-assisted MSK model was implemented in opensim v4.2 [45,46], based on a previously developed and validated SO-based MSK model [26,39]. The knee had 6DoF tibiofemoral and 6DoF patellofemoral joints constrained by 12 bundles of nonlinear elastic ligaments [39,47]. The anterior cruciate ligament was excluded from the model as it was resected during the knee replacement surgery. Knee contact [26,48] was a nonlinear elastic foundation model with Young's modulus = 463 MPa and Poisson's ratio = 0.46 representing the polyethylene inserts [49].

The knee model was integrated into an MSK model [3,26] with 44 muscles on each leg. The model included Hill-type muscle–tendon model with muscle force–length–velocity (FLV) relationship [50] and muscle excitation dynamics (muscle excitation–activation relationship) [43], the latter using a second-order discrete linear model [43]. The EMG-assisted and SO muscle–tendon models [2,51] were implemented with rigid tendons and without parallel elastic components (primarily due to the limitations of opensim solver). For each subject, the MSK model was linearly scaled using the upright double-stance trial. However, knee articulating surfaces were defined subject-specifically using three-dimensional implant geometries from the grand challenge dataset [5,44].

Muscle–tendon parameters, as well as coefficients describing muscle excitation and activation dynamics, often vary across individuals. Thus, accurate estimation of muscle excitations and JCFs becomes challenging. To address this issue, muscle–tendon calibration has been devised to mitigate the effects of muscle–tendon uncertainties on the estimation of muscle excitations and JCFs [33,36,52]. To this end, following the MSK model scaling, Calibrated EMG-Informed Neuromusculoskeletal Modeling Toolbox (ceinms) was used to perform a multi-DoF calibration to optimize muscle–tendon and muscle activation dynamics parameters of all the 44 muscles of the studied leg [2,51]. Optimized parameters for the muscle–tendon model were: maximum isometric force, tendon slack-length, and optimal fiber-length; and for the muscle activation dynamics: recursive filter-coefficients and nonlinear shape-factor [2,43]. Calibration used one trial from each gait style and 5DoF lower limb moments, i.e., hip flexion–extension, abduction–adduction, and internal–external rotation; knee flexion–extension; and ankle dorsi–plantar flexion [2,51].

To lessen the computational burden, an iterative calibration method was used (Fig. 1). This assumed that small changes in the secondary knee kinematics would have negligible effects on the muscle–tendon kinematics and, thus, calibration results. Accordingly, the current calibration iteration was performed using the secondary knee kinematics estimated by muscle–tendon parameters from the previous calibration iteration. Iterations continued until all the changes in the muscle–tendon parameters and muscle forces were less than 5% of the previous iteration.

where $Mn,d,texp$ and $Mn,d,test$ are experimental and estimated moments at time t of nth calibration trial, acting on the model's dth DoF, $max(JCFnest,med)$ and $max(JCFnest,lat)$ are the peak of estimated medial and lateral tibiofemoral JCF for each trial, and $φ$ and $ϕ$ are positive weighting factors. $φ$ and $ϕ$ were set to $10−3$ to scale the JCF to an order of magnitude comparable with joint moments. $Ntrials$ and $NDoF$ are the total number of trials and DoF included in calibration, and $NT$ is the total of time points within nth calibration trial. The cost function in Eq. (1) is shown to improve the accuracy of the estimated knee JCF [12].

where $qpi¨̂$ and $qsi¨̂$ are primary and secondary DoF accelerations produced by a unit force from the ith muscle, $Fimt,limt, and vimt$are muscle–tendon force, length, and velocity of the ith muscle, respectively, and $Vi$ is the ith muscle volume appropriate for the muscle's maximum isometric force [39,54], $q¨ other$ is the nonmuscular generated accelerations including centripetal and Coriolis effects, gravity, ligament forces, contact forces, and external ground reaction forces (applied to the feet). In addition, $eim$ or $ejm$ are estimated excitations of the ith or jth muscles, respectively, $ejEMG$ is the processed EMG of the jth muscle, $Ul$ is contact energy associated with each face of a contact mesh (i.e., within articulating surfaces), $anres$ is the activation level governing the generalized forces from reserve actuators acting on the MSK model's nth kinematic DoF, and $α$⁠, $β$⁠, $γ$⁠, and $ω$ are positive weighting factors that are provided in the Supplemental Materials on the ASME Digital Collection (Sec. S1.3 available in the Supplemental Materials). Rationales for inducing zero-acceleration of knee secondary kinematics (Eq. (2b)) are also provided in the Supplemental Materials (Sec. S1.2.3 available in the Supplemental Materials). In Eq. (2), $Nm$ is the total number of muscles (=44), $NEMG$ is the total number of muscles with measured EMG (=15), $Nfaces$ is the total number of elements within the contact surface, and $Nres act$ is the total number of reserve actuators. Reserve actuators, conventionally used in opensim MSK analyses [25], were added to the net joint forces or moments of each DoF when muscles were incapable of generating joint forces or moments.

The cost function in Eq. (2c) is based on previous studies that aimed to estimate in vivo knee JCF. This cost function estimates the muscle excitations (⁠$eim$⁠) that minimize the squared sum of muscle excitations [1,39,54], follow the measured EMGs [2], minimize knee contact force [12,26], and minimize the use of reserve actuators [26,46]. Notably, there are important differences between EMG-assisted and SO-based neural solutions for the 12DoF knee joint MSK models. First, $β$ = 0 in Eq. (2c) for SO-based solution excludes tracking of EMGs and has been previously introduced and validated for this MSK model [26,39]. Second, when calibrating the EMG-assisted MSK model, Eq. (2) was used at the beginning of each calibration iteration (Fig. 1 and Sec. 2.2.2) to estimate secondary knee kinematics for the current calibration iteration.

The 1DoF knee MSK model [55] was adapted from a previous study that well estimated the in vivo tibiofemoral JCF of the grand challenge dataset [12]. Like the 12DoF knee model, the 1DoF knee model was scaled using the identical double-stance trial. Following scaling, the secondary knee kinematics of the 1DoF knee models were defined as a function of knee flexion angle, which has been previously shown to improve predictions of knee JCF [1,12]. These secondary kinematics were obtained from a passive knee flexion analysis using the MSK model with the 12DoF knee.

The EMG-assisted 1DoF knee MSK model utilized elastic tendons and passive muscle forces. The muscle–tendon parameters were first calibrated, like that described in Sec. 2.2.3, but with just one iteration. The EMG-assisted execution mode of ceinms was then used to estimate muscle excitations and forces, and knee JCF (more details in Sec. S1.3 available in the Supplemental Materials). Execution of the SO-based 1DoF knee MSK models was performed in opensim SO tool v4.2 [45,46], utilizing rigid tendons and excluding passive muscle forces (due to the limitations of opensim solver). Muscle forces were estimated to match the joint external moments while minimizing the sum of squared muscle excitations.

The in vivo medial and lateral compressive tibiofemoral JCF were calculated using the regression equations provided with the fifth grand challenge dataset [5] (Sec. S1.4 available in the Supplemental Materials). The two subjects had different implant designs, while the provided regression equations have been validated only for the fifth grand challenge dataset. Thus, we reported the medial and lateral compressive tibiofemoral JCF only for the fifth grand challenge dataset, i.e., 29 trials of subject PS.

We analyzed 127 successful trials of the stance phase (see Sec. 2.1). Root-mean-square error (RMSE) and coefficient of determination (R²) between experimental and estimated muscle excitation, joint moments, and JCF were calculated for each MSK modeling approach, separately for each subject's trial (including all the time points). RMSE and R² were calculated for the first and second peaks of the JCF, as well as the whole stance phase. We chose this statistical analysis method, instead of, e.g., Statistical Parametric Mapping, to avoid the influence of temporal variations across trials and subjects on the comparisons. Statistical analyses were performed on the RMSE and R² values using one-way repeated-measures analysis of variance (ANOVA) followed by post-hoc test with Bonferroni correction. The significance level was p < 0.05, and statistical tests were performed in spss v27 (IBM Corporation).

Subjects' in vivo secondary knee kinematics were not available; thus, we compared the mean and standard deviation of estimated secondary knee kinematics with the literature.

In general, EMG-assisted models had higher R² and lower RMSE between estimated and measured muscle excitations (i.e., EMG) compared to SO-based models (Fig. 2 and Figs. S1 and S2 available in the Supplemental Materials). For 13 out of 15 muscles with measured EMG, the calibrated EMG-assisted model with 12DoF knee had significantly higher (p < 0.001) muscle excitation R² (estimated versus measured) than the SO-based models and the EMG-assisted 1DoF knee model (Fig. 2 and Fig. S1 available in the Supplemental Materials). For 11 out of 15 muscles with measured EMG, the calibrated EMG-assisted model with 12DoF knee had significantly lower (p < 0.001) muscle excitation RMSE (estimated versus measured) than both uncalibrated EMG-assisted and SO-based 12DoF knee models (Fig. S2 available in the Supplemental Materials). Vastus medialis was the only muscle for which the EMG-assisted models had significantly (p < 0.01) higher muscle excitation RMSE (but not higher R²) compared with SO-based neural solutions (Fig. 2 and Figs. S1 and S2 available in the Supplemental Materials on the ASME Digital Collection). Likewise, the uncalibrated EMG-assisted model with 12DoF knee had significantly lower (p < 0.001) muscle excitation RMSE in 11 of the 15 muscles compared to the same model but with SO-based neural solution (Fig. S2 available in the Supplemental Materials). When comparing calibrated EMG-assisted models, overall, the model with 12DoF knee had significantly higher R² and lower RMSE of muscle excitation in muscles acting on the hip than the model with 1DoF knee (Fig. 2 and Figs. S1–S14 available in the Supplemental Materials).

The reserve actuator moments were near zero within the calibrated 12DoF knee EMG-assisted MSK model except for knee flexion/extension DoF, which presented relatively higher RMSE and lower R² between estimated and experimental moments (Table S1 available in the Supplemental Materials). However, the calibrated 12DoF knee EMG-assisted MSK model had significantly and considerably lower RMSE (p < 0.001) and higher R² in all joint moments compared to 1DoF knee EMG-assisted MSK model (Table S1 available in the Supplemental Materials).

The RMSE between the estimated and in vivo tibiofemoral JCF, for the whole stance or the first and second JCF peaks, was either comparable or significantly lower (p < 0.01) in the calibrated 12DoF knee EMG-assisted MSK model compared to other models (Figs. 3–5 and Fig. S15 available in the Supplemental Materials). In the 12DoF knee models, the RMSE between the estimated and in vivo total, medial, and lateral compressive tibiofemoral JCF over all stance phase were significantly lower (p < 0.001) in the calibrated EMG-assisted compared to the uncalibrated EMG-assisted and SO-based MSK models (Figs. 3–5 and Fig. S15 available in the Supplemental Materials). In the 12DoF knee models compared to other models, the uncalibrated EMG-assisted model had either significantly higher (p < 0.05) or similar RMSE between the estimated and in vivo total, medial, and lateral tibiofemoral JCF (Fig. 5 and Fig. S15 available in the Supplemental Materials).

The RMSE between the estimated and in vivo anteroposterior JCF were significantly lower (p < 0.01) in the MSK models with 12DoF knee joint compared to the 1DoF knee models. However, no significant differences were observed within the MSK models with 12DoF knee in either anteroposterior or mediolateral (i.e., transverse) JCF (Fig. 5 and Fig. S15 available in the Supplemental Materials on the ASME Digital Collection, bottom row). The estimated secondary knee kinematics were slightly different between the two subjects, and minor changes were observed between secondary knee kinematics estimated by the three MSK models with 12DoF knee joint (Fig. 6).

We developed and validated a simulation method comprising a calibrated EMG-assisted MSK model with a 12DoF knee joint using concurrent calculation of secondary knee kinematics and muscle excitations. The knee joint had subject-specific articulating surfaces, elastic foundation contact, and nonlinear elastic knee ligaments [26,39]. Two key features of the developed method were the implementation of a muscle–tendon calibration with a multi-DoF knee joint model, followed by the implementation of assisted estimation of the model's neural solutions based on measured muscle excitations (i.e., EMGs). We compared the estimated muscle excitations and knee JCF against in vivo measurements [5], as well as those estimated by four different conventional MSK modeling approaches (different neural solutions and 1DoF knee joint). Overall, the calibrated EMG-assisted model with 12DoF knee improved the estimation of muscle excitations and transverse tibiofemoral JCF to a greater extent than compressive tibiofemoral JCF.

Muscle recruitment strategy is often noticeably affected by MSK disorders, altered gait patterns, knee pain, muscle weakness, etc. [37,56–60]. This variation was also observed within subjects of our study (Figs. S3–S14 available in the Supplemental Materials, e.g., Fig. S3 versus S6 available in the Supplemental Materials, etc.). Importantly, previous studies have reported SO-based MSK models typically fail to account for these subject- and task-related changes in muscle excitation patterns, especially the muscle co-contractions [2,8,61–64]. Our results (Fig. 2 and Figs. S1–S14 available in the Supplemental Materials) also showed that in more than 70% of the muscles with measured EMG, the SO-based models (with either 12 or 1DoF knee) had significantly (p < 0.05) lower R² and higher RMSE between the measured and estimated muscle excitations, compared to EMG-assisted MSK models.

In patients who undergo total knee arthroplasty, as in our study subjects, quadriceps and hamstring muscles often co-contract and hence are more activated during gait compared to healthy individuals [65]. As illustrated in Figs. S3–S14 available in the Supplemental Materials, the SO-based MSK models, in general, underestimated muscle excitations, including quadriceps and hamstrings, compared to the EMG-assisted models of the study. This may be attributed to the limitations of SO-based neural solutions when estimating muscle co-contractions [8,61,63,64].

Muscle–tendon parameters, such as maximum isometric force, differ across subjects, especially in those with MSK disorders [59,66]. To attenuate the uncertainty of the muscle–tendon parameters, calibration methods have been developed to optimize the subject's muscle–tendon parameters [2,33,37,51,67,68]. In our results, and consistent with previous studies [2,33,37,51,67,68], calibrating the MSK model with 12DoF knee significantly increased R² and decreased the RMSE between the estimated and measured muscle excitations as well as JCF (Figs. 2 and 5 and Figs. S1–S15 available in the Supplemental Materials). In our tests, increasing $β$ in Eq. (2c) decreased the muscle excitation RMSE of the uncalibrated model, but it substantially increased the RMSE between estimated and in vivo JCF despite adjusting the $ω$ weighting factor (Eq. (2c)). To further assess the role of calibration, we performed two additional analyses, i.e., muscle–tendon adjustment without EMGs and EMG-assisted calibration followed by SO-based neural solutions. Details are provided in Secs. S1.5 and S2 available in the Supplemental Materials.

Quadriceps, followed by gastrocnemius, had the highest muscle excitation RMSE for the calibrated EMG-assisted model with 12DoF knee. Yet, the R² values were comparable or significantly (p < 0.001) higher compared to other MSK models of the study (Fig. 2 and Fig. S1 available in the Supplemental Materials). This higher RMSE may be due to the reduced maximum isometric force of these muscles during calibration. The reduced maximum isometric forces were then compensated by more activation of quadriceps in the early stance and gastrocnemius in the late stance to provide the required knee flexion/extension moment (Figs. S3–S14 available in the Supplemental Materials), and consequently, increased the RMSE between the estimated and measured muscle excitations (Fig. S2 available in the Supplemental Materials). Also, we observed a relatively large RMSE between the estimated and experimental knee flexion/extension moment by the EMG-assisted model with 1DoF knee compared to the other joint and MSK models (Table S1 available in the Supplemental Materials). Higher knee moment RMSE means that muscles could not generate the required moment. Thus, increasing the lower bounds of the maximum isometric force of quadriceps and gastrocnemius muscles and decreasing JCF weighting factors (⁠$φ$ and $ϕ$ in Eq. (1)) in calibration may reduce the RMSE between the estimated and measured muscle excitations as well as knee flexion–extension moments.

Excluding muscle activation dynamics increased (<0.1 bodyweight increase, and p > 0.05) the RMSE between the estimated and in vivo lateral tibiofemoral JCF and, to a lesser degree, the medial tibiofemoral JCF (Fig. S16 available in the Supplemental Materials). Including muscle activation dynamics significantly (p < 0.05) increased the RMSE between the estimated and measured muscle excitations in all the muscles with measured EMG (Fig. S17 available in the Supplemental Materials). This increase may be due to the relationship between muscle excitation and muscle activation, in which muscle activations depend on the past values of the neural excitation [43,69], narrowing the solution space (note that tendons were rigid, and the parallel elastic element of muscles was omitted). Likewise, when including or excluding the muscle FLV relationship, we observed minor changes in the estimated knee JCF (Fig. S18 available in the Supplemental Materials) but noticeable changes in the estimated muscle excitations (Fig. S19 available in the Supplemental Materials). Higher RMSE of muscle excitations when including muscle FLV relationship may be attributed to the altered force generation capacity of the muscles at different fiber lengths and contraction velocities, especially given that passive muscle forces were omitted. We also saw convergence difficulties when including the muscle FLV relationship, possibly due to the nonlinear interconnections between muscle forces (⁠$Vi$ in Eq. (2c)) and joint accelerations. Anderson and Pandy [70] also reported equivalent muscle forces during the gait with and without muscle FLV relationship and activation dynamics. To conclude, our analysis showed minor effects of muscle FLV relationship and activation dynamics on the estimated knee JCF in activities like walking.

The RMSE between the estimated and measured tibiofemoral JCF by the calibrated EMG-assisted MSK model with 12DoF knee was 117±33 N and 283±34 N for medial and lateral tibiofemoral JCF, respectively (Fig. 5 and Fig. S15 available in the Supplemental Materials). These are comparable or relatively lower than those of the MSK models reported in the literature (unblinded studies), i.e., ranging from 120 N to 312 N for medial and 220 N to 326 N for lateral tibiofemoral JCF [11,71–73]. It is worth noting that in those studies, only one or a few trials (up to 5) and gait styles (up to 2) of a single subject have been analyzed, while here, we analyzed a total of 127 trials of two subjects.

The calibrated EMG-assisted MSK model with 12DoF knee had significantly (p < 0.01) lower JCF RMSE (estimated versus measured) compared to the SO-based and uncalibrated EMG-assisted MSK models, especially at the second peak of the knee JCF (Figs. 3–5 and Fig. S15 available in the Supplemental Materials). This was noticeable to a greater extent during the treadmill walking of subject PS (Fig. 3). However, the muscles crossing the knee were equally or even more activated by the calibrated EMG-assisted model than the SO-based or uncalibrated models with 12DoF knee (Figs. S3–S14 available in the Supplemental Materials). This contrast can be explained by the calibration results, in which muscle strengths were adjusted to track the muscle EMG and joint moments, leading to a plausible estimation of knee JCF. Specifically, when the model underwent calibration, the maximum isometric force was decreased in gastrocnemius, biceps femoris, and quadriceps but increased in soleus. Hence, decreasing strength of knee crossing muscles led to lower JCF despite greater excitation of these muscles, compared to other MSK models of the study (Fig. 3). Also, increasing soleus strength provided more capacity to counterbalance ankle plantarflexion moment, compensating for the decrease in gastrocnemius strength. Equivalently, the $α$ weighting factor (Eq. (2c)) within the SO-based and uncalibrated MSK models with 12DoF knee was set to encourage activation of the hip (gluteal muscles) and ankle (soleus) crossing muscles but discourage activation of the knee crossing muscles (rectus femoris, gastrocnemius, and hamstrings). This adjustment of $α$ has been previously suggested to decrease the RMSE between the estimated and in vivo tibiofemoral JCF in walking [7,39,74]. In other words, the calibration step tuned the muscle–tendon parameters, implicitly accounting for the optimum adjustment of the $α$ weighting factor in a subject-specific manner (i.e., tuning $Vi$ rather than $αi$in Eq. (2c)).

Assisting the 12DoF knee model with muscle excitations (i.e., EMG) resulted in minor changes in secondary knee kinematics compared to those of the previously validated SO-based model of the study (Fig. 6). This supports the assumption that minor changes in secondary knee kinematics will not affect the muscle–tendon calibration results. Of note, the calibration used simulated annealing [2,75], which approximates the global optimum. Thus, the calibrated muscle–tendon parameters would be minimally affected by (physiologically reasonable) perturbations in the initial guess. Consequently, the calibration took only one iteration (for both study subjects), while the second calibration iteration had <1% changes in the optimized muscle–tendon parameters and muscle forces compared to the first iteration. The implant, due to its design, may have also restrained secondary knee kinematics, although this effect seems to be minor as the transverse forces were not significantly different between the MSK models with the 12DoF knee joint (Fig. 5 and Fig. S15 available in the Supplemental Materials on the ASME Digital Collection). Nevertheless, secondary knee kinematics estimated by MSK models with 12DoF varied considerably across the subjects and compared to those estimated from passive knee flexion (i.e., used for models with 1DoF knee), yet were comparable with the literature [76–79] (Fig. 6). Conversely, increasing the knee DoF in SO-based models did not substantially improve the accuracy of the estimated compressive tibiofemoral JCF, but it did improve the estimated transverse tibiofemoral JCF (Figs. 3–5 and Fig. S15 available in the Supplemental Materials). Likewise, Imani Nejad et al. have reported high errors in the estimated anteroposterior tibiofemoral JCF when using a 1DoF knee model [11], and Marra et al. have reported that using a 12DoF SO-based knee model, compared to a subject-specific 1DoF knee model, did not substantially improve the prediction of compressive knee JCF but enabled personalized prediction of secondary knee kinematics [1].

Overall, if one seeks only to estimate muscle forces and the magnitude of compressive knee JCF (i.e., during walking), the EMG-assisted MSK model with a subject-specific 1DoF knee may suffice [1,8,12]. Nonetheless, our results showed that the EMG-assisted MSK model with 12DoF knee joint provided a more accurate simultaneous estimation of secondary knee kinematics, joint moments, muscle excitations, and compressive as well as transverse knee JCF. Although both EMG-assisted MSK models with either 1 or 12DoF knee joints showed great potential to accommodate patient-specific muscle excitation patterns, the use of the MSK model with 12DoF knee joint showed a capacity to simulate the interaction of ligaments with knee kinematics and kinetics.

Muscle–tendon paths were not subject-specific in the MSK models. Extracting these properties was not possible according to the available dataset, and more importantly, we aimed to develop the model based on conventional and practical inputs to facilitate its application. Also, it has been reported that muscle moment arms may be only slightly affected when using subject-specific geometries for the whole MSK model, as compared to a linearly scaled MSK model (a conventional method used in our study, as well) [80], although this needs further verification [81]. Extracting the subject's muscle–tendon parameters is infeasible in many cases. Instead, we calibrated muscle–tendon parameters to attenuate the effect of the muscle–tendon uncertainties on the estimated knee kinematics and kinetics [33]. Notably, validation of the calibrated muscle–tendon parameters requires extensive efforts and specific data, which was not provided within the grand challenge dataset. Notably, the MSK models of our study used different muscle–tendon models, i.e., the use of rigid/compliant tendons and passive muscle forces. This was due to our aim, i.e., to compare the MSK models' best output according to their previous developments and not to focus on a side-by-side comparison of the models' elements. It is worth mentioning that altering the muscle–tendon model often requires extensive model adjustment. For instance, adding passive muscle forces to the model can lead to compensatory muscle activation and, thus, excessive muscle activity and overestimation of joint contact force [82]. Compliant tendons can substantially increase computation times (up to ∼10 times more than rigid tendons); however, the discrepancies between the estimated muscle force profiles may be negligible [83].

Due to the lack of suitable medical images within the datasets, subject-specific ligament insertion points could not be extracted. Additionally, no information was available on the material properties (e.g., laxity and stiffness) of the ligaments; thus, material properties were used from the literature and a previously validated MSK model [26,39,47,84].

The study data consisted of two subjects performing several activities where muscle excitations and knee JCF had relatively different magnitudes and patterns. However, complementary studies may further evaluate the developed MSK model using a larger cohort and various activities.

Our previous studies on the MSK finite element (FE) analysis of the knee [8,9] showed that supplying a multi-DoF FE model of the knee with an EMG-assisted or SO-based MSK model with 1DoF knee often causes excessive and physiologically unrealistic secondary knee kinematics and hence compromising the simulations. This is indeed due to the overestimated transverse JCFs (Fig. 5 and Fig. S15 available in the Supplemental Materials on the ASME Digital Collection). However, we observed that using the MSK model with 12DoF knee eliminated this issue in MSK FE models [9]. Therefore, we believe the simulation method developed in the current study would substantially facilitate multiscale simulations.

As further developments, we plan to incorporate automated tuning of personalized muscle–tendon paths into the model, and also integrate the model with our previously developed atlas-based MSK FE modeling pipeline [85]. This potentially enables a more accurate estimation of subject-specific tissue-level joint biomechanics, governing tissue adaptation and degradation responses such as collagen network damage and cell death within the knee cartilage [6]. Ultimately, the multiscale analysis workflow aims to estimate knee cartilage degradation, allowing for personalized rehabilitation and gait modification to delay or slow knee osteoarthritis progression [8,86].

In this study, we developed and validated an EMG-assisted MSK simulation method with a 12DoF subject-specific knee joint and concurrent optimization of muscle excitations and secondary knee kinematics. The developed model significantly improved the accuracy of the estimated muscle excitations and knee kinetics compared to previously developed MSK models. Overall, our results encouraged using models with EMG-assisted neural solutions, as well as the calibration of muscle–tendon parameters. The developed model showcased a great potential to estimate knee kinetics of individuals with MSK disorders, which may be altered due to muscle excitation strategy or knee load-bearing tissue (such as ligaments and articulating surfaces).

• Australian National Health and Medical Research Council (NHMRC) (Ideas Grant Award No. APP2001734; Funder ID: 10.13039/501100000925).

• European Union's Horizon 2020 research and innovation program (Marie Skłodowska-Curie Grant Agreement No. 101108335; Funder ID: 10.13039/501100000780).

• Saastamoinen Foundation.

• Finnish Cultural Foundation (Grant No. 00230091; Funder ID: 10.13039/501100003125).

• Novo Nordisk Foundation (Grant No. NNF21OC0065373; Funder ID: 10.13039/501100004191).

• Research Council of Finland (Grant No. 363459; Funder ID: 10.13039/501100002341).

• Sigrid Juselius Foundation (Grant No. 250109; Funder ID: 10.13039/501100006306).

• Australian Research Council Discovery Early Career Research Award (DECRA) (No. DE220101249).

Data provided by a third party listed in Acknowledgements.

More information on the study is provided in the Supplemental Materials. The MSK modeling pipeline developed in this study, integrated into opensim, is available at this link.²