EMG-Driven Musculoskeletal Model Calibration With Wrapping Surface Personalization

Muscle forces and joint moments estimated by electromyography (EMG)-driven musculoskeletal models are sensitive to the wrapping surface geometry defining muscle-tendon lengths and moment arms. Despite this sensitivity, wrapping surface properties are typically not personalized to subject movement data. This study developed a novel method for personalizing OpenSim cylindrical wrapping surfaces during EMG-driven model calibration. To avoid the high computational cost of repeated OpenSim muscle analyses, the method uses two-level polynomial surrogate models. Outer-level models specify time-varying muscle-tendon lengths and moment arms as functions of joint angles, while inner-level models specify time-invariant outer-level polynomial coefficients as functions of wrapping surface parameters. To evaluate the method, we used walking data collected from two individuals post-stroke and performed four variations of EMG-driven lower extremity model calibration: 1) no calibration of scaled generic wrapping surfaces (NGA), 2) calibration of outer-level polynomial coefficients for all muscles (SGA), 3) calibration of outer-level polynomial coefficients only for muscles with wrapping surfaces (LSGA), and 4) calibration of cylindrical wrapping surface parameters for muscles with wrapping surfaces (PGA). On average compared to NGA, SGA reduced lower extremity joint moment matching errors by 31%, LSGA by 24%, and PGA by 12%, with the largest reductions occurring at the hip. Furthermore, PGA reduced peak hip joint contact force by 47% bodyweight, which was the most consistent with published in vivo measurements. The proposed method for EMG-driven model calibration with wrapping surface personalization produces physically realistic OpenSim models that reduce joint moment matching errors while improving prediction of hip joint contact force.


I. INTRODUCTION
C OMPUTATIONAL neuromusculoskeletal models have emerged as inexpensive and objective tools for estimating important biomechanical quantities that are difficult or impossible to measure experimentally, which may eventually facilitate the design of more effective interventions for neuromusculoskeletal disorders [1], [2].Electromyography (EMG)-driven neuromusculoskeletal models in particular have been developed to estimate muscle forces and joint moments computationally using EMG and kinematic data, thereby resolving the muscle redundancy problem while simultaneously calibrating musculotendon properties (e.g., optimal muscle fiber length, tendon slack length) [3], [4], [5], [6].In the majority of EMG-driven modeling methods, muscle-tendon lengths and moment arms (henceforth called "geometries") are determined using geometric models of muscle-tendon pathways implemented within musculoskeletal modeling software such as OpenSim [7], [8] and Anybody [9].These pathways are typically defined using generic geometric representations of muscle origin and insertion points along with muscle wrapping surfaces that specify how muscles wrap around neighboring muscles and/or bones.However, musculoskeletal geometry influences the estimation of muscle forces [6], [10], [11] and joint contact forces [12], [13].Unlike personalization of muscle origin and insertion points, which are relatively easy to identify on computed tomography (CT)derived bone geometry [14], [15], [16], personalization of muscle wrapping surfaces has received limited attention.
Though the most direct way of personalizing muscle wrapping surfaces involves the use of medical imaging data collected from individual subjects, only a small number of studies have followed this approach due to the challenges involved.First, it is often not possible to obtain CT or magnetic resonance imaging (MRI) data from a specific subject.Second, even when it is, segmentation of muscle volumes and subsequent determination of muscle lines of action can be difficult and time consuming [17].Third, even if muscle lines of action can be found from imaging data, the positions, orientations, and dimensions of three-dimensional analytical wrapping surfaces (typically cylinders, spheres, or ellipsoids) must then be determined via complex computational methods [18].Fourth, even if accurate analytical wrapping surfaces can be determined, the geometric calculations involved in finding muscle-tendon lengths and moment arms are too expensive computationally to be used within iterative numerical optimizations that calibrate model parameter values or predict model motion [7], [19].Despite these challenges, several studies have developed methods to customize muscle wrapping surfaces using MRI data [13], [14], [15], [18], [21], [22].Generally, these studies determined muscle wrapping surfaces for only a limited number of muscles using custom computational algorithms.Furthermore, in some studies [17], [21], [22], the resulting muscle wrapping surfaces and associated computational algorithms were not consistent with widely used musculoskeletal modeling software such as OpenSim and Anybody, which further limits usage of built-in analysis tools (e.g., joint reaction force analysis and motion prediction) available in those platforms.
Two studies have used data-driven optimization methods to define subject-specific wrapping surfaces implicitly, where satisfactory agreement with literature-reported experimental musculoskeletal geometry (e.g., moment arms) was achieved [23], [24].However, the optimization methods used in these two studies performed repeated OpenSim Muscle Analyses, which are computationally expensive due to complex anatomical path constraints and iterative muscle wrapping algorithms [7], [19].Consequently, the computational cost of solving these optimization problems was high, limiting their application to a small number of muscles.Additionally, the optimization cost functions sought to match muscle-tendon lengths and moment arms determined a priori from extensive medical imaging datasets, which are difficult to obtain for most experimental situations and are not specific to the individual subjects being modeled.Ideally, muscle wrapping surface parameters could be calibrated using data commonly collected in a gait lab, such as joint kinematics and joint moments, without requiring repeated expensive geometric analyses (e.g., OpenSim Muscle Analyses).
As an alternative to imaging-based methods, EMG-driven model calibration can potentially provide an innovative approach for personalizing musculoskeletal geometries using experimental movement data.EMG-driven modeling uses experimental EMG, kinematic, and external force data along with a geometric musculoskeletal model to estimate muscle forces and joint moments while simultaneously personalizing muscle-tendon model parameter values [3], [6], [25].To date, the majority of EMG-driven modeling studies have used scaled generic musculoskeletal geometry that does not explicitly represent individual subject anatomy, which hinders the accuracy of predicted muscle forces and joint moments [11], [26].A previous study adjusted musculoskeletal geometries within an EMG-driven model calibration process by optimizing polynomial coefficients defining surrogate musculoskeletal geometric models [6].Use of the surrogate geometric models substantially reduced errors between model-predicted and inverse dynamic joint moments.However, the predicted muscle-tendon lengths and moment arms became physically unrealistic for some muscles and could not be converted back into a physical geometric musculoskeletal model, which is important for applications where knowledge of muscle force directions is required.Thus, it would be valuable to have a computationally efficient EMG-driven modeling method that can personalize muscle wrapping surfaces in geometric musculoskeletal models.
This paper presents a novel EMG-driven modeling method that can personalize muscle-tendon model properties and OpenSim cylindrical muscle wrapping surfaces simultaneously.The method was developed and tested using gait data collected from two subjects post-stroke performing treadmill walking at two speeds.Computationally efficient surrogate geometric musculoskeletal models were constructed that accurately predicted muscle-tendon lengths and moment arms for both subjects given joint angles and wrapping surface parameter values as inputs.Wrapping surface parameter values were then optimized as part of the EMG-driven model calibration process.For comparison, EMG-driven model calibration was also performed with no geometric adjustments and using two simpler surrogate geometric models.Each geometric adjustment approach within the EMG-driven model calibration process was evaluated by comparing predicted and inverse dynamic joint moments as well as predicted and scaled generic musculoskeletal geometries.Hip joint contact forces were also compared using models with scaled generic and personalized muscle wrapping surfaces.

A. Experimental Data Collection
Previously published walking datasets collected from two hemiparetic subjects post-stroke were used in this study (Table S1) [6], [27].Motion capture (Vicon Corp., Oxford, United Kingdom), ground reaction (Bertec Corp., Columbus, OH, United States), and EMG (Motion Lab Systems, Baton Rouge, LA, United States) data were recorded simultaneously when subjects walked on a split-belt instrumented treadmill (Bertec Corp., Columbus, OH, United States) with belts tied at their self-selected and fastest-comfortable speeds.All experimental protocols were approved by the University of Florida Health Science Center Institutional Review Board (IRB-01), and both subjects provided written informed consent before performing the experiments.For each subject, 16 channels of EMG data were recorded from each leg using a combination of surface and fine wire electrodes (Table S2).Data from ten gait cycles (five cycles per speed) were randomly selected for analysis of each leg.Raw motion capture and ground reaction data were low-pass filtered at 7/tf Hz, where tf is the period of the gait cycle being processed, while raw EMG signals were high-pass filtered at 40 Hz, demeaned, full-wave rectified, low-pass filtered at 3.5/tf Hz, and finally normalized to maximum value over all experimental gait cycles [6].After being processed, each cycle of data was time-normalized by resampling to 101 time frames from heel-strike (0%) to subsequent heel-strike (100%) of the same leg.

B. Musculoskeletal Model Refinement
A generic full-body musculoskeletal model [28] was employed for development of the approach in OpenSim 4.0 [11], [12].The model possessed six degrees of freedom (DOFs) in each leg, which included three hip DOFs, one knee DOF, and two ankle DOFs.Four sequential steps were performed on the generic model using the OpenSim application Fig. 1.Workflow for personalizing muscle wrapping surfaces while simultaneously personalizing EMG-driven model properties.Briefly, a two-level surrogate muscle-tendon geometric model was developed for each muscle with a cylindrical wrapping surface in OpenSim.The outer-level surrogate model fitted muscle-tendon lengths and moment arms sampled from OpenSim as a function of relevant joint angles using a quintic polynomial, while the inner-level surrogate model fitted the outer-level polynomial coefficients as a function of cylindrical wrapping surface parameter values using a cubic polynomial.Each two-level surrogate model provided fast and reliable estimates of OpenSim muscle-tendon lengths and moment arms.During EMG-driven model calibration, wrapping surface and EMG-driven model parameter values were adjusted concurrently such that errors between model-predicted and inverse dynamics joint moments were minimized.A binary support vector machine (SVM) classification model and a posterior SVM score-to-probability model were developed to identify combinations of wrapping surface parameter values that led to non-physiological muscletendon geometry in OpenSim.These SVM models informed the optimization algorithm when a guessed combination of wrapping surface parameter values was non-physiological.The final step investigated how calculated joint contact force was impacted by wrapping surface personalization.
programming interface (API) with MATLAB (MathWorks, USA) to personalize the model and prepare for EMG-driven model calibration: 1) model scaling to match the individual anthropometry; 2) calibration of lower-extremity joint positions and orientations by tracking the surface marker data [29]; 3) inverse kinematic (IK) analysis to calculate joint angle time histories from marker data; 4) inverse dynamic analysis to derive experimental joint moments using processed ground reaction data and time histories of joint kinematics from IK analysis.After these steps, 34 and 33 muscles per leg were available for analysis for subject S1 and S2, respectively, where 20 muscles in each leg possessed a cylindrical wrapping surface (Table S2).

C. Surrogate Geometric Model Construction
1) Definition of Surrogate Model Structure: Surrogate musculoskeletal geometric models were designed as two-level multivariable polynomials (Fig. 1), which enabled actual OpenSim musculoskeletal geometries to be updated reliably when joint kinematics and OpenSim muscle wrapping surface parameters were changed.Within the outer-level surrogate model, timevarying muscle-tendon length (L mt ), muscle-tendon velocity (v mt ), and muscle moment arms (M A) calculated by OpenSim were approximated using a set of fifth-degree polynomial functions of joint angles (θ 1 , θ 2 , . . ., θ n ) and angular velocities ( θ1 , θ2 . . .θn ) that shared the same polynomial coefficients (c 0 , c 1 . . .c k ): Similarly, muscle moment arms were estimated as partial derivatives of muscle-tendon length with respect to relevant musculoskeletal model generalized coordinates [30]: where n denotes the number of generalized coordinates spanned by the muscle and k denotes the number of terms in the outer-level polynomials, which was different across muscles.The outer-level surrogate models were constructed using all polynomial terms up to degree five.The negative sign in (3) was implemented for consistency with the OpenSim modeling environment.
Where b 0 , b 1 . . .b m are the inner-level polynomial coefficients and m denotes the number of terms in the inner-level polynomials.Similarly, the inner-level surrogate models were Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.
constructed using all polynomial terms up to degree three.
For each inner-level polynomial function, the number of terms was identical since the muscle wrapping surfaces were controlled consistently by six variables.However, the number of polynomial functions in the inner-level surrogate model was muscle-specific, which was dependent on the composition of the outer-level polynomials shown in Eqs. ( 1), (2), and (3).
2) Sampling of Model Poses: During the training phase for constructing outer-level surrogate models, 500 combinations of lower-extremity joint angles for 6 leg DOFs were sampled (named "training poses").As shown in the Fig. S1, the first 400 combinations were randomly selected within the blue belts defined by the minimum and maximum joint angles at each time frame across all experimental walking trials such that the sampled joint poses were physiologically realistic.More specifically, for each of the 100 time frames (representing 1%-100% of the gait cycle) defining each gait cycle after time normalization, 4 combinations of joint angles were randomly sampled within the range of motion for that time frame (named "random poses").In addition, 100 combinations of joint angles from a fastest-speed walking cycle (named "smooth poses") were concatenated behind the "random poses," which facilitated removal of unrealistic wrapping surfaces using a smoothness criterion (Table S6).
3) Sampling of Wrapping Surface Parameters: To account for anatomical constraints and to avoid anatomical obstacles, we employed analytical cylindrical wrapping surfaces to simulate muscle wrapping behavior over joint capsules, bony prominences, deeper muscle layers, and neighboring soft tissues within the curved-line segments [7], [13], [19].
In the generic OpenSim model used for this study, the local reference frame fixed to a given wrapping surface was named the "nominal frame" (Fig. 1).For each wrapping surface, six parameters were introduced to represent adjustments with respect to the nominal frame, including: 1) X -, Y-, and Z-translational offsets (transX, transY, and transZ); 2) Xand Y-rotational offsets following a body-fixed X-Y Euler rotation sequence (rotX and rotY); 3) radius scale factor (Rr).For each wrapping surface, 1000 combinations of parameters were sampled in a six-dimensional hypercube space using a Hammersley quasirandom sequence [31] (named "training samples").The upper and lower bounds for rotX (−20∼20 deg), rotY (−20∼20 deg), and Rr (0.8∼1.2) were consistent across wrapping surfaces, while the bounds for transX, transY and transZ were customized to account for the subject's anatomy.Further details can be found in Table S8 of the supplementary materials.
4) Analysis of Muscle Geometry: Muscle analyses were performed to calculate muscle-tendon lengths and moment arms for all 500 sampled poses using all 1000 muscle wrapping surface parameter combinations, resulting in a total of 500,000 muscle analyses performed through the OpenSim-MATLAB API.For muscles that wrap around a cylindrical surface, musculoskeletal geometries were computed for each of the 1000 sampled muscle wrapping surface configurations, while for muscles that do not wrap around any cylindrical surfaces, muscle-tendon geometries were computed only once.

5) Removal of Non-Physiological Wrapping Geometries:
During simulated motions with the OpenSim model, nonphysiological wrapping behaviors can occur due to unrealistic wrapping surface configurations.For example, if a sampled muscle wrapping surface is too far away from or too close to the associated muscle path, the muscle path will not wrap around the cylindrical surface or will penetrate the wrapping surface through the entire simulated motion or at certain poses.In addition, muscle wrapping errors sometimes occur in an OpenSim model whereby the muscle path wraps around the circumference of the cylinder due to unrealistic relative locations of intermediate via points, which results in unrealistic overestimates of muscle-tendon length [24].Therefore, three criteria were defined to determine whether or not a sampled wrapping configuration was physiologically valid, which depended on measurements of deviation, variability, and smoothness derived from the collected musculoskeletal geometries (Table S6).For each wrapping surface, all sampled combinations of wrapping surface parameters were analyzed following the above-mentioned method and classified as a "non-physiological" or "physiological" wrapping configuration.
6) Classification of "Non-Physiological" Wrapping Geometries: To prevent EMG-driven model calibration from converging to a non-physiological set of wrapping surface parameters, we needed to inform the optimizer in an iterative and efficient manner whether wrapping surface parameters were "physiological" or "non-physiological."In this study, support vector machine (SVM) models were adopted to recognize non-physiological combinations of wrapping surface parameters [32].For every wrapping surface, after all physiologically unrealistic combinations of wrapping parameters were identified and labeled, an SVM model was trained with 1000 "training samples" of wrapping surface parameters and corresponding labels, where all six wrapping parameters were taken directly as "features" and third-degree polynomial kernel functions were consistently applied across all wrapping surfaces.Next, to implement efficient classification within EMG-driven model calibration, an optimal sigmoid function was fitted for each wrapping surface to calculate the continuous probabilities defining whether a combination of wrapping parameters is likely to be classified as "non-physiological" from the SVM-predicted scores s j [33]: This transformation function was dependent on the slope M and the intercept B, which were fitted for each wrapping muscles individually.A threshold value was calculated when s j was set as 0 for the j th wrapping surface.If the computed probability was larger than the threshold value, the set of wrapping parameters was classified as a "non-physiological" configuration and was rejected as an unrealistic solution within EMG-driven model calibration.surfaces.Outer-level polynomial coefficients defined in Eqs.
(1) and (3) were calculated by solving linear least squares problems, where the known inputs were sampled joint angle, muscle-tendon length, and muscle moment arm time histories.This solution process fit sampled muscle-tendon length and moment arm data simultaneously.Next, inner-level polynomial coefficients defined in Eq. ( 4) were calculated by solving additional linear least squares problems, where the known inputs were all "physiological" combinations of sampled wrapping surface parameters.The accuracy of surrogate musculoskeletal geometric models was evaluated using 250 new combinations of wrapping surface parameters sampled in the same six-dimensional space (named "validation samples") and 500 new sets of model poses (named "validation poses").OpenSim muscle analysis was again performed for each validation sample of wrapping surface parameters at each validation combination of model poses.Non-physiological combinations of wrapping surface parameters in the validation set were also removed by following the same three criteria before comparing the model-predicted and OpenSimcalculated musculoskeletal geometries.

D. Model Calibrations With Geometric Adjustments
1) EMG-Driven Musculoskeletal Model: In this study, an EMG-driven model calibration process was employed to adjust musculoskeletal geometries while simultaneously calibrating EMG-driven model parameter values.The combined calibration process extends a previously published EMG-driven model calibration process [6], [34].Given processed EMG data, joint kinematics, and musculoskeletal geometries (muscle-tendon lengths, velocities, and moment arms) as inputs, the EMG-driven model predicted net joint moments in the lower extremities following three steps.First, a nonlinear muscle activation dynamics model with electromechanical delay was used to determine muscle activation profiles from muscle excitations.Muscle excitations were derived by multiplying processed EMG signals by a set of scale factors that account for unknown maximum excitation levels [6], [35], [36].Second, a Hill-type muscle-tendon model with a rigid tendon predicted muscle forces from muscle activations, muscle-tendon lengths, and muscle-tendon velocities [37], [38].Last, net joint moments were calculated by combining estimated muscle forces with muscle moment arms.The joint moment produced by a muscle spanning a particular joint were represented as: where M is the moment generated by the muscle about the joint, M A is the moment arm of the muscle about the same joint, F m o is the maximum isometric force of the muscle, a is the muscle activation, α is the pennation angle of the muscle, Lm is the normalized muscle fiber length, and ṽm is the normalized muscle fiber velocity.f l Lm and f v ( ṽm ) represent the normalized muscle active force-length relationship and the normalized muscle force-velocity relationship, while f p Lm defines the normalized passive force-length relationship.
2) Muscle Geometric Adjustment Approaches: This study explored three different approaches for adjusting musculoskeletal geometries within EMG-driven calibration.First, as the primary focus of this study, musculoskeletal geometries for muscles with cylindrical wrapping surfaces were predicted by the proposed two-level surrogate musculoskeletal geometric models defined by Eqs. ( 1) to (3), both of which are functions of joint kinematics and muscle wrapping surface parameters.We allowed cylindrical wrapping surface parameters to be altered during the optimization process such that the outer-level polynomial coefficients and muscle-tendon geometries were updated accordingly.This method was named physical geometric adjustment (henceforth called "PGA").Second, musculoskeletal geometries of all muscles in the EMG-driven model were calculated using outer-level surrogate models that were controlled by joint kinematics alone.Instead of being controlled by wrapping surface parameters, independent outer-level polynomial coefficients were directly tuned during optimization to achieve geometric adjustment.This method was named surrogate geometric adjustment (henceforth called "SGA") [6].Third, musculoskeletal geometries of only muscles with cylindrical wrapping surfaces were determined using outer-level surrogate models, where outer-level polynomial coefficients were again design variables in the optimization.However, unlike SGA, this approach did not adjust musculoskeletal geometries for muscles without a cylindrical surface.This method was named limited surrogate geometric adjustment (henceforth called "LSGA").In addition, to establish a control for comparison, we also performed optimizations where no musculoskeletal geometries were adjusted (henceforth called "NGA").
3) EMG-Driven Model Calibration: Within EMG-driven model calibration, a sequence of non-linear optimizations was performed such that the net joint moments produced by the model (M mod i ) matched experimental joint moments obtained from OpenSim Inverse Dynamics analyses (M exp i ) as closely as possible.The primary cost function term was defined as follows: where N is the total number of joints.The variables calibrated for each muscle-tendon actuator by the optimization process included 1) electromechanical delays, activation dynamics time constants, activation non-linearization shape factors, and EMG scale factors that described subject-specific EMG-toactivation relationships; 2) optimal muscle fiber length and tendon slack length that described subject-specific muscle force-generating ability; 3) variables involved in geometric adjustment: outer-level polynomial coefficients for LSGA or SGA, and wrapping surface parameters for PGA.The cost function not only aimed to minimize errors in joint moment curves but also included penalty terms as "soft constraints" to restrict deviations of model parameter values and curves away from the initial model.Furthermore, inequality constraints, referred to as "hard constraints," were incorporated into the problem formulation to ensure that normalized muscle lengths remained within realistic ranges, to prevent sign reversal of muscle moment arms for both SGA and LSGA approaches, to enforce customized bounds on the muscle wrapping surface parameters, and to ensure that wrapping surface parameter values were physiological for PGA.Further details can be found in the supplementary materials.
In this study, experimental data from ten gait cycles (five cycles per available speed) were adopted for model calibration (henceforth called "calibration trials"), while experimental data from an additional ten gait cycles were used for model validation (henceforth called "validation trials").The optimized model parameters, including activation dynamics model parameters, Hill-type muscle-tendon model parameters, and wrapping surface model parameters, obtained during the calibration process remained unchanged throughout the validation process.
For the PGA implementation, the upper and lower bounds for the parameters of each wrapping surface were selected to be consistent with the upper and lower limits for sampling.For a given wrapping surface, a nonlinear constraint was added to the optimization problem formulation to prevent the solution from converging to a "non-physiological" wrapping configuration, which required the probability of a set of parameters being classified as 'non-physiological' to be below the corresponding threshold value.

E. Hip Contact Force Calculation
This study also investigated how personalized muscle wrapping surfaces impacted estimated hip joint contact forces (HJRF) during walking.Hip joint contact forces were calculated using OpenSim for models with generic and personalized wrapping surfaces across all calibration walking trials.During analysis, muscle forces predicted using NGA-based EMGdriven models were used to actuate the model with generic wrapping surfaces, while muscle forces estimated using PGAbased EMG-driven models were employed to actuate the optimized model with personalized wrapping surfaces.

F. Evaluation Metrics
Several common metrics were employed to quantify the performance of surrogate musculoskeletal geometric models and the reliability of EMG-driven model predictions generated using geometric adjustment.First, root mean square error (RMSE) values between the estimated musculoskeletal geometries using the two-level surrogate geometric models and the generated musculoskeletal geometries using OpenSim muscle analysis were computed to evaluate the performance of the proposed two-level surrogate musculoskeletal geometric models.Second, mean absolute error (MAE) values between the personalized musculoskeletal geometries following EMG-driven calibration and the generic OpenSim-derived musculoskeletal geometries were calculated to quantify the changes of muscle-tendon lengths and moment arms with different geometric adjustment approaches.Third, MAE values between inverse dynamic and model predicted net joint moments were derived to evaluate the EMG-driven joint moment tracking performance using different geometric adjustment strategies.Last, to illustrate the difference between NGA-and PGA-based HJRFs, peak values in resultant HJRF near 20% of the gait cycle were calculated and compared.HJRFs could not be calculated for the SGA and LSGA methods since physical muscle lines of action are not produced by those methods.
Multiple statistical analyses were performed to assess whether the calculated metrics were statistically different across different geometric adjustment approaches.First, we performed a single-factor analysis of variance (ANOVA) with a Tukey-Kramer post-hoc analysis on MAE values to investigate whether the geometric adjustment approaches had a significant influence on the accuracy of joint moment prediction.Second, we performed paired t-tests to investigate whether joint moment estimation accuracy and magnitudes of muscle forces and muscle activations were significantly different between each geometric adjustment approach and the case where no geometric adjustments (NGA) were made within EMG-driven model calibration.Third, we ran paired ttests to compare changes in musculoskeletal geometries from original values for each of the three geometric adjustment methods.Last, we assessed whether peak values of HJRFs were statistically different between PGA and NGA using paired t-tests.All statistical analyses were performed in MAT-LAB, and significance levels were set at p < 0.05.

III. RESULTS
For both calibration and validation trials, PGA produced significantly more accurate joint moment estimates for all three hip DOFs than did NGA (Table I, p < 0.05).Knee moments estimated by PGA were also significantly more accurate than those from NGA for the fastest speed ( p < 0.05), while ankle moment estimates from PGA were not significantly different from those found by NGA.For joint moment estimates that were not statistically different between PGA and NGA, PGA always produced lower mean MAE errors.In contrast, SGA and LSGA produced significantly more accurate joint moment estimates for all six DOFs than did NGA (Table I and Fig. S5).The one exception was the ankle plantarflexion-dorsiflexion moment estimated by LSGA for the validation trials at selfselected speed, through the mean MAE error for LSGA was lower than for NGA.
Absolute changes of musculoskeletal geometries using different geometric adjustment methods were quantified by MAE values calculated between the optimized (SGA, LSGA, or PGA) and scaled generic (NGA) musculoskeletal geometries (Table S3).Paired t-tests indicated that for musculoskeletal geometries associated with the three hip DOFs, there were no significant differences in absolute changes among SGA, LSGA, and PGA ( p > 0.05).However, for musculoskeletal geometries associated with the knee DOF and two ankle DOFs, changes in moment arms were significantly larger for LSGA and SGA than for PGA ( p < 0.05).For changes in muscle-tendon lengths, average absolute changes were less than 0.17 cm for SGA, 0.41 cm for LSGA, and 0.88 cm for PGA.For changes in moment arms, average absolute changes were less than 0.85 cm for SGA, 1.22 cm for LSGA, and 0.84 cm for PGA.Taking the results of iliacus as an example (Fig. 2), there was an observable discrepancy Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.

TABLE I MEAN ± STANDARD DEVIATION JOINT MOMENT TRACKING ERRORS (QUANTIFIED USING MEAN ABSOLUTE ERROR) PRODUCED BY EMG-DRIVEN MODEL CALIBRATION (CAL) AND VALIDATION (VAL) PERFORMED USING DIFFERENT GEOMETRIC ADJUSTMENT APPROACHES
for hip arms among the three approaches, where LSGA and SGA produced larger changes than did PGA when comparing calibrated and original moment arms in terms of both shape and magnitude.Also, personalized locations, orientations, and radii of wrapping surfaces produced by PGA differed significantly from those of the scaled generic models ( p < 0.05, Fig. 2 and Table S4).
Hip joint reaction forces (HJRFs) were estimated using models with original and calibrated wrapping surfaces, respectively, when walking motions for analysis were matched (Fig. 3).The profiles of resultant overall hip contact forces calculated by the two methods were similar in shape.However, compared to the non-adjusted wrapping surfaces, PGA-calibrated wrapping surfaces yielded significantly lower average peak values near 20% of the gait cycle.

IV. DISCUSSION
This study incorporated personalization of cylindrical muscle wrapping surfaces into the EMG-driven musculoskeletal model calibration process, which produced subject-specific EMG-driven model parameter values and physical OpenSim models with personalized muscle wrapping surfaces simultaneously.On the one hand, our method could help avoid the need to collect medical imaging data and provide a new path for personalizing a large number of muscle wrapping surfaces simultaneously.On the other hand, scaled generic muscle-tendon lengths and moment arms have been reported to be one of the primary sources of joint moment matching errors for EMG-driven model calibration [6], [11], [39].Our method allowed wrapping surface parameter values to be adjusted during optimization, which significantly improved joint moment prediction due to the use of personalized musculoskeletal geometries.
We developed two-level surrogate musculoskeletal geometric models that allowed actual OpenSim musculoskeletal geometries to be updated when both joint kinematics and OpenSim muscle wrapping surface parameter values were changed concurrently (Fig. S2).Two-level surrogate models were employed for several reasons.First, two-level structures could reduce the number of possible terms that appeared in each polynomial function, which enabled us to fit musculoskeletal geometries without encountering the issue of overfitting.Second, outer-level polynomial functions were constructed using only joint kinematics as inputs, allowing physical constraints between muscle-tendon length and moment arms to be imposed in a straightforward manner [6].Third, we could flexibly adjust the power of variables in each category when implementing linear regression fitting at the corresponding level.Fourth, once muscle wrapping surface parameter values were personalized, outer-level polynomial coefficients became readily available, where outer-level surrogate models could be adopted for further applications with extremely low computational cost.
Our PGA-based adjustment of wrapping surface parameter values within EMG-driven model calibration was isolated from OpenSim to reduce computational cost, which required a continuous and accurate model capable of determining whether a wrapping surface configuration was "physiological" or "non-physiological."SVM models were advantageous for accurate classification of "non-physiological" combinations of wrapping surface parameters due to observations that in most prescribed six-dimensional spaces for wrapping surface parameters, "non-physiological" or "non-anatomical" wrapping surface configurations normally happened when multiple parameter values were close at the edges of the design space, where the boundary plane between two classes was feasible to describe.Moreover, posterior SVM score-to-probability models consistently output probability values between 0 and 1 across wrapping surfaces, which facilitated implementing nonlinear constraints during optimization.Fig. 2. Representative muscle-tendon geometries for the iliacus obtained from EMG-driven model calibration performed using three different geometric adjustment approaches for both subjects.The last column to the right illustrates the difference between original wrapping surfaces (semicylinders in green, NGA) and personalized wrapping surfaces (semi-cylinders in red, PGA) associated with the iliacus (named "IL_at_brim").Red lines in OpenSim models represent the muscle-tendon path of the iliacus.
Among the three geometric adjustment approaches investigated in this study, PGA yielded the smallest improvement in joint moment prediction accuracy while SGA produced the largest improvement (Table I).This observation can be explained by the number of design variable available for adjustment and number of muscles adjusted by each of the three methods.First, for any muscle, SGA and LSGA introduced substantially more design variables into the optimization problem formulation than did PGA.Taking iliacus as an example, its geometries were controlled by only six wrapping surface parameters with PGA, while 41 independent outer-level polynomial coefficients controlled them with SGA and LSGA for the left leg of S1.Second, geometric adjustments were performed on a larger number of muscles with SGA than with LSGA and PGA.SGA permitted muscle-tendon geometries to change for all 34 muscles for S1 and all 33 muscles for S2.However, LSGA and PGA modified geometries only for 14 muscles per leg muscles that contacted cylindrical wrapping surfaces during walking.
Muscle-tendon lengths and moment arms derived with PGA were close the original generic values in terms of both magnitude and shape and were more realistic physically than those produced by SGA or LSGA, especially at the hip (Fig. 2 and Table S3).Unlike SGA and LSGA, which were implemented with outer-level surrogate models alone, PGA employed two-level surrogate models, which maintained not only the physical relationship between muscle-tendon lengths and moment arms [30] but also physical constraints inherently imposed by realistic wrapping surface configurations.Also, PGA produced relatively large deviations of musculoskeletal geometries for only a small number of hip muscles (e.g., iliacus, psoas, and gluteus maximus middle), whereas SGA and LSGA spread out their adjustments over a larger number of muscles.
This study develop the PGA method to address two limitations found in the SGA and LSGA methods.First, the PGA method eliminates non-anatomical musculoskeletal geometry solutions which occur for some muscles with the SGA and LSGA methods.In the PGA method, the polynomial coefficients in the outer-level surrogate models are no longer treated as independent variables.Instead, these coefficients are controlled by a set of six wrapping surface parameters.This modification reduces the likelihood of the model converging to non-anatomical solutions for muscle-tendon lengths and moment arms, as evidenced in SGA and LSGA solutions.Second, the PGA method makes it possible to convert a musculoskeletal geometry solution into physical wrapping surfaces in an OpenSim musculoskeletal model, which is not possible with the SGA and LSGA methods.This capability could broaden the application of personalized musculoskeletal models to tasks such as custom implant design.
Overestimated hip joint reaction forces have been reported in several papers that estimated muscle forces using static optimization [12], [13], [40] or EMG-informed models [41] applied to musculoskeletal models with generic wrapping surfaces.Wesseling et al. [13] and De Pieri et al. [12] reported that the use of MRI data to increase the level of subject-specificity in hip muscle wrapping geometries resulted in decreased resultant hip contact forces during walking.We observed the same trend when our resultant hip joint contact forces calculated by PGA were compared with those produced by scaled generic models (Fig. 3), suggesting that some level of hip geometry personalization may be important for estimating hip joint contact forces reliably.There were two potential reasons behind these similar findings.First, our study and the two previous studies [12], [13] found that subject-specific muscle wrapping surfaces increased moment arms of hip flexor muscles (e.g., iliacus and psoas) in the sagittal and frontal planes (Fig. 2), where these muscles make major contributions to hip contact force [11].Larger moment arms with the same net hip moments mean that smaller forces are required from those muscles, thereby decreasing hip joint contact forces (Fig. S3) [13].Second, muscle wrapping surface configurations had a considerable effect on muscle lines of action and directions of muscle forces [12], [40].Therefore, the components and orientations of the resultant hip contact force were affected.
The improvement in joint moment estimation accuracy within our EMG-driven modeling method can be attributed to adjustment of musculoskeletal geometries rather than muscle activations and muscle force-generating properties.We formulated our optimization problem to penalize deviations of model parameter values away from designated reference values, which reduced the influence of musculoskeletal geometric adjustments on the solutions found for other model parameter values.The estimates of muscle activations (Fig. S4) and muscle-tendon model parameters (Fig. S6) demonstrated a high degree of similarity across various geometric adjustment approaches, suggesting that musculoskeletal geometric adjustements had little influence on the estimation of muscle activations and muscle force-generating properties.The primary cause of changes in muscle forces could therefore be attributed to alterations in muscle-tendon lengths, which lead to modifications in normalized muscle fiber lengths.
In addition to the future research directions suggested above, our study possesses several limitations that would be worth addressing in the future.First, we could not validate the personalized muscle wrapping surfaces with our current datasets, since the necessary medical imaging data were not available from our two subjects.Second, joint moment prediction errors were lower with SGA and LSGA than with PGA, which demonstrated that there is still room for improvement by personalizing other musculoskeletal geometric features, such as muscle origin and insertion points, or the number and paths of multiple heads used to represent a single muscle.Third, we modeled only a single task -walking at two speeds.Movement datasets collected during multiple daily activities functional activities other than walking have been used for EMG-driven model calibration [42], [43].If richer datasets were used for model calibration, such as movements with larger ranges of motion (e.g., squatting and stair climbing), then the variation of muscle activity and musculoskeletal geometry would be more pronounced, which in turn would improve the robustness and reliability of the calibrated model.Fourth, our developed approach concentrated on personalizing cylindrical wrapping surface geometry, which sped up musculoskeletal geometry calculations compared to ellipsoidal wrapping surfaces [28].It would be worthwhile to explore how well a similar approach would work for other types of wrapping surfaces (e.g., ellipsoids, spheres).Last, to increase computational efficiency, we used a rigid tendon model in our Hill-type muscle-tendon models.Previous researches have shown that rigid and compliant tendon models yield nearly identical muscle force estimates for slow movements like walking at a healthy speed but different muscle force estimates for faster movements like running [44], [45].Both of our stroke subjects walked at slow speeds, suggesting that use of a compliant tendon model was appropriate.Nonetheless, it would be valuable to extend our method to Hill-type muscletendon models with compliant tendons so that our EMG-driven model calibration process could be applied to tasks involving fast movements.

V. CONCLUSION
The EMG-driven model calibration method presented in this study made several valuable contributions.First, two-level surrogate musculoskeletal geometric models were developed such that actual OpenSim musculoskeletal geometries were accurately, reliably, and efficiently estimated when joint kinematics and OpenSim muscle wrapping surface parameters were changed.Second, SVM classification models were created automatically that could accurately determine whether a given combination of wrapping surface parameters would lead to "non-physiological" wrapping behaviors in OpenSim.Third, as the core of this study, an EMG-driven model calibration method was developed that adjusted muscle wrapping surface parameter values to improve the prediction of both joint moments and hip contact forces during walking.Fourth, physical OpenSim models were updated accordingly with EMG-driven model calibration results such that muscle wrapping surfaces represented increased subject specificity.Our method provided not only a new approach for personalizing a large number of muscle wrapping surfaces, which requires less effort compared with construction of geometric models from MR and/or CT data, but also improved joint moment prediction accuracy when performing EMG-driven modeling.

Fig. 3 .
Fig. 3. Hip joint contact force calculated for all calibration trials using the original OpenSim model with nominal wrapping surfaces (NGA) and the calibrated OpenSim model with personalized wrapping surfaces (PGA).
7) Training and Validation of Surrogate Models: Training of surrogate musculoskeletal geometric models began at the outer level for all muscles in the lower extremity musculoskeletal model, including those with and without wrapping Authorized licensed use limited to the terms of the applicable license agreement with IEEE.Restrictions apply.