Modeling Musculoskeletal Dynamics during Gait: Evaluating the Best Personalization Strategy through Model Anatomical Consistency

and personalization strategies in musculoskeletal (MSK) models and gait analysis. A low anatomical consistency raises concern about MSK model outputs credibility. Abstract: No consensus exists on how to model human articulations within MSK models for the analysis of gait dynamics. We propose a method to evaluate joint models and we apply it to three models with different levels of personalization. The method evaluates the joint model’s adherence to the MSK hypothesis of negligible joint work by quantifying ligament and cartilage deformations resulting from joint motion; to be anatomically consistent, these deformations should be minimum. The contrary would require considerable external work to move the joint, violating a strong working hypothesis and raising concerns about the credibility of the MSK outputs. Gait analysis and medical resonance imaging (MRI) from ten participants were combined to build lower limb subject-speciﬁc MSK models. MRI-reconstructed anatomy enabled three levels of personalization using different ankle joint models, in which motion corresponded to different ligament elongation and cartilage co-penetration. To estimate the impact of anatomical inconsistency in MSK outputs, joint internal forces resulting from tissue deformations were computed for each joint model and MSK simulations were performed ignoring or considering their contribution. The three models differed considerably for maximum ligament elongation and cartilage co-penetration (between 5.94 and 50.69% and between − 0.53 and − 5.36 mm, respectively). However, the model dynamic output from the gait simulations were similar. When accounting for the internal forces associated with tissue deformation, outputs changed considerably, the higher the personalization level the smaller the changes. Anatomical consistency provides a solid method to compare different joint models. Results suggest that consistency grows with personalization, which should be tailored according to the research question. A high level of anatomical consistency is recommended when individual speciﬁcity and the behavior of articular structures is under investigation.


Introduction
Musculoskeletal (MSK) models allow the evaluation of quantities otherwise not directly measurable in vivo, such as muscle and joint internal forces, during gait and other daily tasks. Additionally, they bring the advantages of predictions into the process of designing and personalizing interventions and medical devices. Trusted by the recognized need for personalized medicine [1], the number of studies involving MSK models has increased exponentially in the last twenty years [2].
(BTS, Smart DX, MA, USA 100 Hz) with two force plates (Kistler, Winterhur, Switzerland, 1 kHz), and an 8-camera system (Vicon, MX, Oxford, UK, 200 Hz) and two force plates (AMTI, OR6, Albion, Australia, 1 kHz), respectively. For each participant, three walking trials at self-selected speed were used for the analysis. The marker set was a combination of the markers from the Vicon Plug in gait protocol (Vicon, Oxford, UK) and the modified Oxford Foot Model (mOFM) protocol [18], for a total of thirty-eight markers. All gait markers were marked on the skin of the participants and replaced with MRI-visible markers retained during the subsequent imaging acquisition for data registration.
Three MRI sequences were acquired: a lower-limb e-THRIVE scan with 1 mm in-plane resolution and a 1 mm slice thickness; a regional foot and ankle multi-slice multi-echo 3D Gradient Echo (3D_mFFE_WATS) scan with water-only selection in the sagittal plane with 1 mm slice thickness, 20.5 mm interslice gap and 0.5 mm in-plane resolution; a regional foot and ankle 3D short T1 inversion time inversion recovery fast field echo (T1w_TSE) scan in the sagittal plane with 2 mm slice thickness, 21 mm interslice gap and 0.6 mm in-plane resolution.
Experimental data were collected from both limbs, however in one case data from one limb were discarded due to issues with the gait analysis acquisition, whereas in three other cases data from one limb were discarded due to the quality of the foot imaging, resulting in bilateral information being only available for 6 participants, hence leading to a total of 16 datasets.
The 3D representations of the volunteers talus and tibia were obtained from manual segmentation of the T1w_TSE (through MITK freeware software). The bones and the cartilages of the tibio-talar articulation were reconstructed and the attachment regions of the tibio-calcaneal (TiCa) and fibulo-calcaneal (FiCa) ligaments were also identified.

Scaled-Generic Hinge Joint (SGJ)
The joint model was obtained from the generic OpenSim model gait2392 [19]. The original model was modified adding a set of virtual markers matching the experimental marker set previously described, and then scaled relying on the Scaling Tool in Open-Sim [20] (Figure 1). Markers and adjusted weighting factors were used to scale the model using one static gait trial, in compliance with OpenSim best practice recommendations [3].

Morphologically Fitted Hinge Joint (MFJ)
The morphological fitting approach (described in [15,16]) was adopted to fit a cylinder to the superior surface of the 3D talus geometry segmented from the MRI. The cylinder axis was therefore assumed as the ankle axis of rotation under the assumption of hinge-like behavior of the articulation (Figure 1).

Maximum-Congruence Floating Axis Joint (MCJ)
Following the approach proposed in [21], the individual 3-dimensional tibio-talar motion was computed as the trajectory that optimizes the transmission of contact forces, i.e., that maximizes joint congruence [22]. The result is a 1 degree-of-freedom floating-axis joint (Figure 1), in which all components of the talar position and orientation with respect to the tibia are coupled with the ankle flexion angle.

Evaluation of Anatomical Consistency and Computation of Joint Internal Forces
The three ankle models were used to quantify ligament elongation and maximum cartilage co-penetration, and to estimate ligament and contact forces for each model and each foot.
Both TiCa and FiCa ligaments were represented by their single most isometric fiber [13], which was shown to guide the ankle natural motion [23]. For the same reason, only these two ligaments were considered in this analysis. The fibers were defined based on the MCJ motion, as this model proved to predict the individual natural motion [21]. For the Appl. Sci. 2021, 11, 8348 4 of 14 computation of ligament elongation and cartilage co-penetration, talus and calcaneus were moved rigidly as a single bone with respect to the fixed tibia and fibula, considering each ankle model; co-penetration was evaluated only between the talar and tibial cartilages.
i. 2021, 11, x FOR PEER REVIEW 4 of 15 Figure 1. Definition of the personalized MSK model: anatomical information from a regional MRI taken at the foot and a whole lower-limb MRI were merged by registering the segmented tibia bones. Muscle attaching points and marker location with respect to the bones were derived from the lower-limb MRI, while MFJ and MCJ were defined based on detail representation of the ankle articulating surfaces.

Morphologically Fitted Hinge Joint (MFJ)
The morphological fitting approach (described in [15,16]) was adopted to fit a cylinder to the superior surface of the 3D talus geometry segmented from the MRI. The cylinder axis was therefore assumed as the ankle axis of rotation under the assumption of hinge-like behavior of the articulation (Figure 1).

Maximum-Congruence Floating Axis Joint (MCJ)
Following the approach proposed in [21], the individual 3-dimensional tibio-talar motion was computed as the trajectory that optimizes the transmission of contact forces, i.e., that maximizes joint congruence [22]. The result is a 1 degree-of-freedom floating-axis joint (Figure 1), in which all components of the talar position and orientation with respect to the tibia are coupled with the ankle flexion angle. Definition of the personalized MSK model: anatomical information from a regional MRI taken at the foot and a whole lower-limb MRI were merged by registering the segmented tibia bones. Muscle attaching points and marker location with respect to the bones were derived from the lower-limb MRI, while MFJ and MCJ were defined based on detail representation of the ankle articulating surfaces.
Simplified ligament and cartilage representations were adopted in the attempt to simplify the comparison among joint models. Ligaments were modeled using a bilinear elastic characteristic: where L, L 0 and k lig are the ligament length, resting length and stiffness, respectively. In particular, k lig was taken equal to 164.3 N/mm [24] for both ligaments, since no linear stiffness values was found in the literature for the TiCa. This was done under the assumption that the two ligaments have the same function in guiding the ankle motion [25]. The resting length L 0 was assigned as the ligament length measured at the MRI scan. Cartilages were modeled using the elastic foundation approach, according to which where V is the co-penetration volume and k cont is the cartilage spring stiffness [22]. k cont was determined according to [26] as: here, the Young modulus E was assumed equal to 12 MPa [27], the Poisson ratio ν equal to 0.42 [28] and the combined cartilage thickness h = 2 * 1.3 = 2.6 mm [29]. The resultant force and the resultant moment (with respect to the joint fixed or floating axis) of the ankle internal forces associated with the tissue deformation were computed according to the same ligament and contact models above described.

Lower-Limb Musculoskeletal Model
The full lower-limb MSK models including 4 segments (pelvis, femur, tibia, foot) were built in NMS Builder (Bologna, Italy) [30] according to the step-by-step pipeline proposed by Modenese et al. [15] (Figure 1). 3D bone geometries and soft tissue were segmented using a statistical shape model approach developed after Steger et al. [31] from the 3D_mFFE_WATS (foot and distal tibia) and from the e-THRIVE (proximal tibia, femur and pelvis), and regional and lower-limb MRIs were registered by aligning the tibial bones segmented in the two scans via iterative-closest-point algorithm in MeshLab. Volume of bone and soft tissue segmentations were used to assign inertial properties to the model [15,32].
Within each KC, the hip and knee were modeled as ideal ball-and-socket and hinge joints, respectively. The ankle joint was modeled using the three methods described in Section 2.2, which led to three lower limb models for each dataset, differing only for the definition of the ankle axis ( Figure 1). Muscle paths were mapped from the generic gait2392 model [19] using a supervised atlas registration technique [15] and muscle parameters were scaled from gait2392 model default values according to body mass (maximal isometric force) or musculotendon unit length (optimal fiber length, tendon slack length) [15]. Reserve actuators were included for each dof of the model in order to support muscle forces by contributing to the joint torque generation when equilibrium could not be reached. Reserve actuators were assigned a maximum force equal to 1 N to favor the recruitment of muscles.
In this way, a first family of three MSK models were built for each individual (M1), differing by the ankle joint representation (M1 SG , M1 MF , M1 MC ). To evaluate the impact of anatomical inconsistency, a second family of three additional models were built for each individual (M2), in which the ankle internal resultant and resultant moment computed at II.3 for each ankle model (M2 SG , M2 MF , M2 MC ) were also added to the dynamic simulations ( Figure 2).

Dynamic Simulations of Gait
A total of 96 (16 datasets × 3 ankle joint models × 2 family of MSK models) monolateral lower-limb MSK models underwent a standard inverse pipeline in OpenSim employing the MATLAB (R2020b, The Mathworks, Portola Valley, CA, USA) API. Joint angles were estimated using the Inverse Kinematics Tool. Inverse Dynamics analysis was performed considering only the ground reaction force (GRF) and inertial contribution for models belonging to M1, while ankle internal forces associated with tissue deformation were also added to the simulations for models belonging to M2. Finally, the absolute differences between the output of M1 and M2 were quantified. assigned a maximum force equal to 1 N to favor the recruitment of muscles.
In this way, a first family of three MSK models were built for each individual (M1), differing by the ankle joint representation (M1SG, M1MF, M1MC). To evaluate the impact of anatomical inconsistency, a second family of three additional models were built for each individual (M2), in which the ankle internal resultant and resultant moment computed at II.3 for each ankle model (M2SG, M2MF, M2MC) were also added to the dynamic simulations ( Figure  2).

Figure 2.
Visualization of difference loading acting on M1 (on the left) and M2 (on the right). Green arrow represents ground reaction force while blues arrows represent internal forces, computed as explained in Section 2.3.

Dynamic Simulations of Gait
A total of 96 (16 datasets × 3 ankle joint models × 2 family of MSK models) monolateral lower-limb MSK models underwent a standard inverse pipeline in OpenSim employing the MATLAB (R2020b, The Mathworks, Portola Valley, CA, USA) API. Joint angles were estimated using the Inverse Kinematics Tool. Inverse Dynamics analysis was performed considering only the ground reaction force (GRF) and inertial contribution for models belonging to M1, while ankle internal forces associated with tissue deformation were also added to the simulations for models belonging to M2. Finally, the absolute differences between the output of M1 and M2 were quantified.

Statistical Analysis
All output variables were tested for normality and non-parametric statistics were conducted in case of non-normal distribution.
Ligament elongation and cartilage co-penetration estimated with the three ankle models over the ankle range of motion (ROM) were compared by means of Kruskal-Wallis test (α = 0.05) and Wilcoxon post hoc test with Bonferroni correction (α = 0.0167) in MATLAB.
For the three different models in the M1 family, joint angles, moments (normalized by body mass and height [33]), muscle forces (normalized by body weight, BW) of seven muscle crossing the ankle joint, and magnitude of the resultant joint contact forces (normalized by BW) were compared by means of a non-parametric one-way repeated measure ANOVA (α = Figure 2. Visualization of difference loading acting on M1 (on the left) and M2 (on the right). Green arrow represents ground reaction force while blues arrows represent internal forces, computed as explained in Section 2.3.

Statistical Analysis
All output variables were tested for normality and non-parametric statistics were conducted in case of non-normal distribution.
Ligament elongation and cartilage co-penetration estimated with the three ankle models over the ankle range of motion (ROM) were compared by means of Kruskal-Wallis test (α = 0.05) and Wilcoxon post hoc test with Bonferroni correction (α = 0.0167) in MATLAB.
For the three different models in the M1 family, joint angles, moments (normalized by body mass and height [33]), muscle forces (normalized by body weight, BW) of seven muscle crossing the ankle joint, and magnitude of the resultant joint contact forces (normalized by BW) were compared by means of a non-parametric one-way repeated measure ANOVA (α = 0.05). This was implemented using the MATLAB package SPM1D [34] for time-dependent statistics based on Statistical Parametric Mapping. The same analysis was repeated for the output of M2 (except for the joint kinematics and hip and knee moments which are not affected by the ankle internal forces). Post hoc analyses were based on non-parametric two-tailed paired t-test with Bonferroni correction (α = 0.0167). The same analysis was repeated for the output of M2 and for the differences between M1 and M2 (except for the joint kinematics and hip and knee moments which are not affected by the ankle internal forces).

Model Verification
Inverse Kinematics output was assessed by verifying that peak tracking error (i.e., the error between virtual and experimental markers) was below 2 cm for all the markers in the lower limb [3]. Static optimization was considered successful if the contribution of reserve actuators was below 10% of the peak joint moment for the specific dof [35].

Anatomical Consistency
When considering the average ligament elongation across the cohort, no differences were observed between models for the TiCa (Figure 3), with values being overall within the 5% limit defined by the isometric behavior [23]. For the elongation of the FiCa, MFJ and MCJ performed similarly with significant differences only observed during dorsiflexion ( Figure 3). Conversely, the SGJ showed significant differences from the other models for most of the ROM, resulting in a non-physiological behavior with an average maximum elongation of 20% at the extremities of the flexion/extension ROM (Table 1).
in the lower limb [3]. Static optimization was considered successful if the contribution of reserve actuators was below 10% of the peak joint moment for the specific dof [35].

Anatomical Consistency
When considering the average ligament elongation across the cohort, no differences were observed between models for the TiCa (Figure 3), with values being overall within the 5% limit defined by the isometric behavior [23]. For the elongation of the FiCa, MFJ and MCJ performed similarly with significant differences only observed during dorsiflexion ( Figure 3). Conversely, the SGJ showed significant differences from the other models for most of the ROM, resulting in a non-physiological behavior with an average maximum elongation of 20% at the extremities of the flexion/extension ROM (Table 1). When considering maximum values for single individuals, SGJ exceeded 30% of elongation and MFJ showed elongations above the threshold of 5% for both FiCa and TiCa (Table 1). Conversely, the MCJ respected the isometric behavior for both ligaments and all participants, remaining within physiological working conditions.  When considering maximum values for single individuals, SGJ exceeded 30% of elongation and MFJ showed elongations above the threshold of 5% for both FiCa and TiCa (Table 1). Conversely, the MCJ respected the isometric behavior for both ligaments and all participants, remaining within physiological working conditions.
In terms of average co-penetration across the cohort, the three models differed from each other during dorsiflexion (Figure 3). MCJ showed co-penetration below −0.2 mm and preserved contact during the whole joint motion, as dictated by the model definition. SGJ and MFJ showed excessive co-penetration at the extremities of the ROM (up to −1.5 and −4.1 mm, respectively).
When looking at individual cases, MCJ still performed in compliance with the assumption of cartilage non-penetration, whereas MFJ and SGJ showed highly non-physiological values up to −2.4 and −5.4 mm, respectively (Table 1).

M1 Outputs
The results of the ANOVA confirmed small but significant within-subject differences in joint angles and moments ( Figure 4A,B and Figure S3) associated with the use of different ankle models. These differences were not limited to the ankle joint but propagated proximally through the KC. Post hoc analysis highlighted that the three models differed for the hip kinematics, mostly during mid-and late stance and swing; also, M1 MC kinematics differed from the others during mid-and late stance and mid-swing at the knee, and from M1 MF during mid-and late stance at the ankle. Joint moments showed biggest differences at the ankle joint, where the M1 SG differed from the others throughout the whole stance; M1 MF knee moments differed from those of the other models during mid-and late stance; hip moments presented differences across the three models limited to small portions of the gait cycle.

M2 Outputs
Anatomical inconsistency resulted in considerable ankle internal forces ( Figure S1 for the M1SG, causing the static optimization not being able to converge to a solution in 6 cases. For the remaining cases, ankle torque, muscle forces and ankle contact force changed significantly reaching non-physiological patterns and values. Consequently ANOVA and post-hoc test highlighted significant difference between M2SG and the othe two models throughout most gait cycle (Figures 5 and S3). Internal forces also induced significant differences between M2MF and M2MC in the Tibialis Posterior and Peroneu For what concerns muscle forces ( Figure 4C), post hoc analysis highlighted significant differences mostly attributable to the M1 SG , i.e., Medial Gastrocnemius, Soleus, Lateral Gastrocnemius muscles, during stance and late swing. Peronei muscles presented differences between each of the three models during mid-and late stance. M1 SG showed significant differences with respect to the other two models also in the joint contact forces at both the knee (mid-and late stance) and ankle (stance and late swing), also differing from the M1 MF during hip mid-and late stance ( Figure 4D).

M2 Outputs
Anatomical inconsistency resulted in considerable ankle internal forces ( Figure S1) for the M1 SG , causing the static optimization not being able to converge to a solution in 6 cases. For the remaining cases, ankle torque, muscle forces and ankle contact forces changed significantly reaching non-physiological patterns and values. Consequently, ANOVA and post-hoc test highlighted significant difference between M2 SG and the other two models throughout most gait cycle ( Figure 5 and Figure S3). Internal forces also induced significant differences between M2 MF and M2 MC in the Tibialis Posterior and Peroneus Longus and Brevis, as well as at the ankle contact force during swing phase ( Figure 5  Absolute differences between M1 and M2 outputs were over one order of magnitude larger in SGJ than in MCJ and MFJ, as reported in Table 2 for both average cohort and Individuals (Peak values for significant output for the two family of MSK models on average cohort and for individuals are reported in Tables S2 and S3, respectively). For example, differences between M1 and M2 were 0.0, 0.1, and 3.7 N/kg for average ankle joint moment, 0.1, 2.2, and 43.5 BW for average ankle joint contact force, and 0.0, 0.2 and 3.6 BW for average Soleus muscle force for the MCJ, MFJ and SGJ, respectively. An even higher impact was found when looking at differences for individual cases, with a maximum difference for the ankle joint moment of 0.0, 0.1, and 13.3 N/kg, for the ankle joint contact force of 0.8, 6.8, and 120.8 BW, and for the Soleus muscle force of 0.0, 0.2, and 2.5 BW for the MCJ, MFJ and SGJ, respectively.

Discussion
This study compared three ankle joint models with an increasing level of personalization to evaluate their anatomical consistency, defined as their capability of fulfilling the hypothesis of workless motion, and hence assess their credibility.
The reported results proved that anatomical consistency increases with the level of personalization of the joint model. On the cohort, SGJ showed elongation of the FiCa ligament up to 20%, associated with a highly non-physiological behavior. Both the fixed rotational axis of MFJ and SGJ resulted in considerable cartilage co-penetration during dorsiflexion, with the SGJ reaching up to −4.1 mm. The effects of joint personalization were more evident when looking at single individuals: MCJ showed good anatomical consistency for each of the modeled dataset; MFJ exceeded the experimental threshold of 5% of ligament isometric elongation, still remaining within physiological range, however, resulting in a −2.4 mm maximum of cartilage co-penetration; SGJ simply resulted in ligament elongation and cartilage co-penetration not compatible with the ankle physiology.
When included within a subject-specific MSK model of the lower-limb, the three ankle joints caused only small, yet significant, differences in the overall limb kinematics ( Figure 4A). Similarly, small but significant differences could be found for joint torques ( Figure 4B) Differences were more evident when analyzing the muscles: the M1 SG predictions were overall significantly higher than those from the other two models, with a maximum difference above 1.4 BW at the Soleus. M1 MF and M1 MC performed similarly on the cohort, with minor differences for the peroneal muscles ( Figure 4C).
When the effects of anatomical inconsistency are made explicit by introducing the joint internal forces associated with ligament and cartilage deformations ( Figure S1), the MSK model output changed considerably. M2 SG did not converged in six cases and outputs were altered to non-physiological values and patterns. Some muscle forces (Tibialis Posterior and Peroneal), ankle joint contact forces and ankle moments associated to M2 MF and M2 MC changed significantly as a result of accounting for the tissue deformations ( Figure 5C,D). The introduction of internal forces also altered the results for the individual participants at the joint moments and the muscle and joint contact forces (Table 2). For SGJ, maximum individual differences between M1 and M2 was 13.3 N/Kg for the ankle moment, 2.5 BW for the Soleus force and 120.8 BW for the ankle contact force. MFJ was more stable, reaching however a maximum individual difference for the ankle contact forces of 6.8 BW. MCJ was less affected, with an maximum individual difference for the ankle contact forces of 0.8 BW. These differences provide a quantification of the credibility of each joint model: differences between M1 and M2 should in fact be negligible for a behavior consistent with joint anatomy.
The results of the analysis presented in this paper show that similar joint models may lead to similar MSK output, and yet differ considerably in terms of estimate of the deformation of the articular tissues. The evaluation of anatomical consistency can thus be used to discriminate among different joint models and ultimately choosing the most suitable level of personalization of the kinematic chain within MSK models. Our results suggest that scaling of a generic KC provides acceptable results for the general characterization of a motor task, the pattern and the amplitude of all relevant quantities being indeed comparable with the other two models here considered. If a specific cohort has to be investigated and the analysis is extended to joint reactions, morphological fitting appears as a more reliable method, offering also a good compromise between accuracy and easiness in model definition. However, if a characterization of the single individual and its articular structures is needed, a level of personalization similar to the one obtained by congruence maximization is strongly recommended.
Anatomical consistency may also help to extend the knowledge of joint models' sensitivity. SGJ and MFJ are both simple hinges, differing only for the location of the revolute axis with respect to the patient articulation. As previously reported in the literature, MSK models are quite insensitive to small variation in the definition of axis location of hinge joints [11]. Indeed, the two models led to similar joint kinematics ( Figure S2), resulting in similar muscle lever arm and therefore in similar muscle forces. However, we showed that small differences in relative bone displacements may give rise to considerable ligament and cartilage contact forces: if these forces are considered all model outputs change considerably. This implies that, when investigating the articular response to load, a high kinematic accuracy needs to be achieved through an adequate model personalization. This may have important consequences for the design of prostheses and other medical devices, often relying on numerical estimation of joint ligament and contact forces. A highly anatomically consistent joint model would guarantee more credible and accurate MSK outputs, possibly opening the way to a full customization of the implant to the patient needs. In this perspective, MRI based personalization provided the best results in terms of anatomical consistency, showing the crucial role of an accurate representation of the patient articular surfaces.
The role of rigid joints within MSK multibody approach deserves to be discussed. Real articulations are non-rigid and an accurate investigation of internal forces requires detailed models of ligaments and cartilage deformation under loads. These models are however computationally demanding. In addition, the number of parameters that must be determined and tuned grows with the model complexity. On the other hand, rigid joint models may provide accurate estimate when properly defined. Marra and co-workers [36] showed indeed that a simple hinge predicted knee joint contact forces with the same fidelity of a 6-dof model with deformable ligaments for a prosthetic knee. Moreover, rigid body models can be used as preliminary step to inform deformable models, simplifying their definition and application. For example, the maximum congruence approach here used for the MCJ models may provide individual joint kinematics, as well as insertion and reference length for the isometric fibers, which are needed for investigation of the ligaments [37]. The same data can be used to design more complex joints such as rigid mechanism [38,39], whose constraints can eventually be released to cope with deformable joints [40,41].
This work has its limitations. The experimental design of the study, with data collected in two different laboratories, possibly introduced additional variability in the experimental data. Collected data were however post-processed consistently. From a methodological perspective, several assumptions were made. First, simplified ligament and cartilage representations were adopted in the attempt to simplify the comparison among joint models. Indeed, only the two major ankle ligaments, TiCa and FiCa, were included, they were represented by their most isometric fiber, assuming a bilinear elastic characteristic with the same stiffness for both ligaments; articular contact was modeled through the elastic foundation approach. When evaluating anatomical consistency, the effects of joint distraction or slack ligaments were ignored. Both would result in joint instability, hard to consider within a rigid multibody model. Second, the lower-limb MSK models had muscle geometry and musculotendon parameters scaled from average literature values obtained from older cadavers [19,42]. Their application to study a juvenile cohort is therefore questionable. However, muscle properties were kept identical across the three modeling methods, therefore not affecting the comparison carried out in the study. Furthermore, marker-based scaling of a generic MSK models assumes that skeletal anatomy scales linearly between individuals. While this is accepted for adults, its suitability is still unclear for juvenile individuals [43]. Therefore, the adoption of a scaled ankle axis could explain the significant differences observed for the SGJ. Third, a comprehensive validation of the MSK models could not be carried out based on the existing data. In fact, no electromyography data nor dynamic fluoroscopy data were available as part of the retrospective study. However, MCJ was previously validated [21] and the respect of ligament isometry found in this study can be considered as an indirect validation of the model. Last, while our analysis focused on the effect of anatomical inconsistency on the ankle joint, a comprehensive approach should look at the whole KC including joints with more complex kinematics, such as the knee or the articulations within the foot.

Conclusions
Quantification of anatomical consistency in terms of ligament and cartilage deformation associated with joint motion provides a solid criterion for the comparison of different joint models and personalization strategies. The smaller this deformation, the more consistent the model is with the patient anatomy and MSK model hypotheses. A low anatomical consistency raises concern about MSK outputs credibility.
Among the investigated personalization methods, our analysis suggests that markerbased scaling is suitable for characterizing generic muscle behavior associated to a given motor task. If interested in joint analysis on a specific population, morphological fitting appears a more reliable approach. Finally, if interested in individual specificity and the behavior of articular structures, a level of anatomical consistency similar to the one achievable trough congruence maximization is recommended.

Supplementary Materials:
The following material is available online at https://www.mdpi.com/ article/10.3390/app11188348/s1, Table S1: Anthropometric parameters of the participants, Table S2: Peak values for significant output over the gait cycle on average cohort for M1 and M2 models and maximum absolute variation between them, Table S3: Peak values for significant outputs over the gait cycle for individuals for M1 and M2 models and maximum absolute variation between them, Figure S1: Resultant force (a) and resultant torque (b) of the ankle internal actions due to ligament and cartilage deformation for the average cohort, plotted versus ankle flexion angle for the three ankle models. Bottom gray bars represent the flexion angles where post-hoc tests were significant, Figure S2: Motion componentes of the tibio-talar motion (abduction/adduction-AA; internal/external rotation-IEt; antero/posterior translation-X; prossimo/distal translation-Y; medio/lateral translation-Z) for the average cohort, plotted versus ankle flexion angle for the three ankle models. Bottom bars represent the flexion angles where post-hoc tests were significant, Figure S3: Differences between M1 and M2 outputs for the average cohort: forces exerted for significant muscles at the ankle (A); joint moments (B); joint contact force (C). Bottom gray bars represent the region of the gait cycle where post-hoc tests were significant.