Optimization of Material Coefficients in the Holzapfel-Gasser-Ogden Material Model for the Main Four Ligaments of the Knee Joint-A Finite Element Study ()
1. Introduction
Computational biomechanics, specifically FE analysis, has become an indispensable tool that assisted in-vivo and in-vitro experiments over the past few decades. Moreover, factors such as increased cost-effectiveness, the non-invasive nature of the studies, and ethical purposes of reducing the risk to human subjects have complemented its ability to decode the behavior of complex biological tissues [1] [2]. Nevertheless, a reliable FE model should accurately represent the joint geometry and material properties [3]. Determining the material properties, especially those for the soft tissues, has been a challenge in the biomechanical FE studies [4]. Knee joint computational modeling is not an exception; kinematics and the overall biomechanics of a knee joint are directly correlated to the properties assigned to the respective cruciate and collateral ligaments. Hence, it is essential to use an appropriate material model to capture the constitutive behavior of the ligaments.
A wide range of material properties has been used for ligament modeling in the literature [5]. Earlier ligament models employed one-dimensional elements with nonlinear elastic force-elongation equations for the ligaments’ behavior [6] - [11]. Other studies have used the neo-Hookean [12] [13] [14] [15] [16] or Mooney-Rivlin [17] [18] material models, which are mainly isotropic and might not accurately represent the direction dependency of the ligaments. Lately, anisotropic properties based on the strain energy density function were used in several studies to suitably represent the ligaments’ anisotropic behavior [19] - [25]. The equation consists of an incompressible neo-Hookean component to represent the ground substance of the connective tissues [26] and a component to model the fibrous behavior, such as in Holzapfel-Gasser-Ogden (HGO) [27] and other custom-developed constitutive models [21] [22] [23] [24]. However, none of the studies [19] [20] that used HGO-based constitutive models for the knee ligaments have accounted for the ligaments’ viscoelastic behavior [28] [29] [30]. This work aimed to overcome these limitations by obtaining the HGO coefficients through optimization while considering the ligaments’ viscoelastic properties and comparing knee behavior using the optimized properties and Neo- Hookean materials. This work differs from Kiapour et al. [19] and Beidikhti et al. [25] in that they both investigated the effect of ligament modeling techniques, i.e., using different ligament geometry representations and materials, on knee FE simulations, which makes it challenging to decide what portion of the changes in outcomes attributes to the geometry selection and what portion is related to the applied material model. The former [19] compared isotropic non-linear elastic properties while modeling the ligaments as a group of 2D truss elements and HGO material models applied to 3D knee ligaments. The latter [25] compared HGO in knee FE with 3D ligaments and simplified spring elements with non-linear stiffness. This study uses the same 3D geometries of the ligaments to better delineate the effect of different material models for both Neo-Hookean and HGO models. Viscoelastic effects were also included in both cases.
Mechanical tests have been used in the literature to obtain stress-strain curves for different hard and soft tissues, e.g., bones, ligaments, menisci, and cartilages [31]. Utilizing the accurate properties from mechanical tests for modeling ligament behavior in FE software platforms such as Abaqus requires three essential steps: 1) finding the appropriate test data, 2) choosing the proper material model to represent the anisotropic hyper-elastic behavior of the ligaments, and 3) obtaining the coefficients to use in the selected material model using curve fitting techniques. This study chose the HGO model for the main four tibiofemoral ligaments. HGO model coefficients for both pairs of cruciate and collateral ligaments were obtained using optimization schemes which utilized load-displacement data from three separate mechanical test studies for ACL [32], PCL [24], and MCL [33]. Limited representative test data were available for the LCL; thus, the data from the MCL mechanical tests were adopted to determine LCL coefficients.
The purpose of this study was to identify the coefficients in HGO model for the knee joint cruciate and collateral ligaments based on the mechanical test results. Viscoelastic behavior was also taken into account for a closer representation of the physiological biomechanics of the knee joint. The obtained material properties were then used in a dynamic FE model of the knee joint [34] under loading conditions from an in-vitro experiment [35]. Resultant joint kinematics and ligament strains were compared with FE model which uses neo-Hookean properties for the ligaments. We hypothesized that the FE model with optimized ligament properties would better match the experimental joint kinematics and ligament strains when used in a dynamic finite element model entailing other properties and characteristics from the literature. Having access to these optimized coefficients will help other researchers to be able to speed up their analyses by directly applying these values to the ligaments in their knee finite element simulations whenever they are using the HGO material model in combination with the ligaments’ visco-elastic properties, without having to perform the cumbersome optimization steps for all these four knee ligaments.
2. Materials and Methods
2.1. Material Property Optimization
MR images of the knee cruciate and collateral ligaments from a healthy female subject (23 yr, 1.71 m, 60.3 kg) were segmented in Mimics v15.0 (Materialise, Leuven, Belgium), smoothened in Geomagics Studio (3D Systems, Rock Hill, South Carolina), and meshed in Hypermesh v12.0 (Altair, Troy, MI, USA) using four-node tetrahedral elements. The 3D meshes were then imported to Abaqus/ explicit 6.14-5 (SIMULIA, Providence, RI, USA) for mesh convergence. Mesh was refined in localized regions until the stress and strain differences between two subsequent meshes were less than 5%. Mechanical test data from literature were used to optimize the hyper-elastic coefficients for the HGO material model [36]. For ACL, Woo et al. [32] experimental data was used, Wan et al. [24] data was used for PCL. In order to obtain MCL coefficients, Quapp and Weiss’s [33] experimental results were used, and for LCL same experimental data as MCL was used to derive the coefficients.
Ligaments were modeled as bone-ligament-bone structures in the FE simulations to replicate mechanical test conditions (Figure 1). Two block geometries were created in Abaqus to simulate the bone attachment sites for each end of the ligament. The distal bone block was fixed and the proximal block was placed under a controlled displacement via a reference point coupled to the nodes of the bone block. If literature test curves reported load-displacement data, they were converted to stress-strain curves based on the length and area of each ligament in the experiments to be applicable for optimizing the properties of ligaments with different lengths and geometries. Then they were converted back to load-displacement format based on the geometries of the ligaments used in the FE for optimizations. The displacement was applied in the same direction and orientation as per each experiment. Step time for applying the displacement in each case was calculated to reproduce the same loading rates as in the experiments. Reaction force outputs acting on the reference point were extracted to generate load-displacement curves for each ligament.
Figure 1. Bone-ligament-bone setup for ACL to simulate tension test experiments in Abaqus.
ACL and PCL were modeled with two fiber families, while one fiber family was considered for MCL and LCL. Simulia products (Abaqus and Isight) were used to perform the optimizations. Hooke-Jeeves optimization algorithm [37] was selected to fit the load-displacement curve to that of the literature by minimizing the root mean square error and achieving a correlation greater than 0.95. Hooke-Jeeves optimization algorithm was previously used to estimate material model coefficients of biological tissues, such as liver tissue [38] and residual limb bulk soft tissue [39]. The strain energy potential equation used in Abaqus for the anisotropic hyper-elastic HGO material model is based on Holzapfel et al. 2000 [27], and Gasser et al. 2006 [36] as follows.
with
,
where the five coefficients of C10, D, K1, K2, and κ are temperature-dependent material parameters that are the user-defined inputs to the software. C10 is the Neo-Hookean constant, D is the inverse of bulk modulus and controls the incompressibility, and κ determines the dispersion of the fibers.
Initially, the values of the C10 (Neo-Hookean constant) and D (the inverse of the bulk modulus) were kept fixed to the values available in the literature (Table 1), and optimizations were performed by varying the other three parameters. This approach worked well for both ACL and PCL. However, this did not allow for an appropriate fit for the load-displacement curves of the collateral ligaments so C10 and D were allowed to vary during the optimization process to obtain a satisfactory fit. The optimized coefficients were also tested in a simulated compression to monitor the ligaments’ behaviors under compressive loads. The volumes of the ligaments used in this study are also reported in Table 2 to give the reader an idea of how close they are to their models’.
2.2. FE Model Development
Details on the FE models’ development and validation were presented in our
Table 1. Neo-Hookean material coefficients for the knee ligaments (Pena et al., 2006).
Table 2. Ligaments volumes in mm3 for the model used in optimizations.
previous work [34] and is briefly described in this section. Following IRB approval, MR images of a healthy female subject (23 yr, 1.71 m, 60.3 kg) were used for creating the 3D geometry. All bony structures and soft tissues including femur, tibia, patella, fibula, medial and lateral menisci (MM & LM), articular cartilages (femoral; FC, medial tibial; MTC, lateral tibial; LTC, patellar; PC), and ligaments (anterior and posterior cruciate; ACL & PCL and medial and lateral collateral; MCL & LCL) were manually segmented in Mimics. The smoothing process of the 3D surfaces was performed in Geomagics Studio. Hypermesh was used for the subsequent mesh generation. Four-node tetrahedral elements were created and assessed for mesh quality [40] in Hypermesh, factoring into consideration the elements warpage, aspect ratio, Jacobean, and tet collapse. The final meshes had only less than 1% of elements for each part possessing warpage > 5, aspect ratio > 3, Jacobean < 0.7, and tet collapse < 0.3. Model assembly (Figure 2) and mesh convergence analysis were performed using Abaqus/explicit. The four main ligaments (ACL, PCL, MCL, and LCL) were coupled at their insertion sites to the bones. Femoral, tibial, and patellar cartilages were tied to the bones underneath them. Seven frictionless surface-to-surface interactions were defined at the contact between articular cartilages, among each other and with the menisci. Additional surface-to-surface interactions were assigned to contacts between ACL-PCL, ACL-femoral notch, and tibia-MCL. Meniscal horn attachments were modeled as kinematic couplings between the horns and the insertion site on the tibial plateau. Connector elements were used as peripheral attachments connecting the menisci to the tibia articular surface. Other capsule structures were modeled as connectors with nonlinear load-displacement data. Material properties were derived from the literature. In brief, bones were modeled as linearly elastic materials with different properties for cortical and cancellous bones [3]. Menisci were modeled as transversely isotropic materials [41]. The neo-Hookean material model was used for the articular cartilage with values
Figure 2. Knee finite element model including 3D structures for bones, cartilages, menisci, and four main ligaments, and 2D connectors for patellofemoral attachments and knee capsule.
calculated from young’s modulus and Poisson’s ratio [42]. The four main ligaments were modeled as hyper-elastic materials according to HGO [36], with coefficients obtained through the curve fitting process based on experimental data, as is described above in the material properties optimization section. Lastly, viscoelastic properties were applied to cruciate and collateral ligaments by assigning stress relaxation data using time-based prony series [43].
2.3. FE Simulations
Load cases examined included various combinations of knee abduction moment (KAM), internal tibial rotation moment (ITR), and anterior tibial shear force (ATS) followed by an axial compression force equal to the impact load of half a bodyweight dropping from a 30 cm height, simulating bipedal landing from a jump (Table 3). All analyses were done with the knee flexed to 25 degrees. Muscle forces of 441N were applied to quadriceps; for counterbalance, the same amount of force (441N) was considered for hamstring muscles, equally distributed among the lateral and medial hamstring groups. The biceps femoris long head was used for applying the lateral hamstring load, and the semitendinosus, semimembranosus, and gracilis were used for the medial side. Analyses were done using Abaqus explicit solver [44]. Readers may refer to Erbulut et al. 2021 [34] for more details on the loads and boundary conditions.
FE simulations were performed once using the optimized properties and once using the Neo-Hookean properties; viscoelastic behavior was considered in both cases. Outputs for knee kinematics (valgus/varus angle, internal/external tibial rotation angle, anterior/posterior tibial translation, and superior/inferior tibial translation) along with ACL and MCL strains were extracted and compared.
3. Results
3.1. Ligament Material Property Optimization
Comparison between the load-displacement curves obtained from literature with those of the optimized properties showed a high correlation (r > 0.95) for all the four knee ligaments (Table 4 & Figures 3-6). Moreover, these properties were
Table 3. Loading scenarios simulating sub-failure loadings of knee abduction moment (KAM), anterior tibial shear force (ATS), and internal tibial rotation moment (ITR), determined with regard to the in vivo population percentage (%) (Bates et al., 2017), followed by the axial impact of a drop load equal to half a bodyweight.
Table 4. Holzapfel-Gasser-Ogden material model coefficients for the knee main ligaments obtained via optimization, and the correlation coefficients.
Figure 3. ACL curve fitting from Simulia iSight. The horizontal axis shows displacement in mm, and the vertical axis displays force in N. The black curve is the experimental data from the literature, and the blue curve is the load-displacement curve from the optimization.
Figure 4. PCL curve fitting from Simulia iSight. The horizontal axis shows displacement in mm, and the vertical axis displays force in N. The black curve is the experimental data from the literature, and the blue curve is the load-displacement curve from the optimization.
Figure 5. MCL curve fitting from Simulia iSight. The horizontal axis shows displacement in mm, and the vertical axis displays force in N. The black curve is the experimental data from the literature, and the blue curve is the load-displacement curve from the optimization.
Figure 6. LCL curve fitting from Simulia iSight. The horizontal axis shows displacement in mm, and the vertical axis displays force in N. The black curve is the experimental data from the literature, and the blue curve is the load-displacement curve from the optimization.
used in five subject-specific FE knee models and produced the kinematics and ligament strains largely in range of 1.5 standard deviation (SD) of the mean in-vitro data [35], which supports the initial hypothesis. The outcomes looked at for validation included internal tibial rotation angles, knee abduction angles, anterior tibial translations, axial compressions and ACL and MCL strains. When used in a simulated compression, these properties showed no resistance under compressional loads, which ensures the incompressibility of the ligament behavior.
3.2. FE Simulations of Bipedal Landings
Kinematic outputs from the experimental data were reported with respect to baseline kinematics at 25˚ knee flexion. Results from the FE simulations follow this same convention (Table 5 & Table 6, and Figures 7-10). Results are presented at 33 ms, 66 ms, and 100 ms after axial impact loading, the determined time course for ACL rupture following landing [45] [46].
Table 5. Knee joint kinematics for load cases 1 - 4 at 33, 66, and 100 ms after initial ground contact. Green marks the values in the range of one standard deviation from the mean in in-vitro experiments, blue marks the values within the 1.5 standard deviation range, and orange marks the values outside this range.
Table 6. Knee joint ligament strains for load cases 1 - 4 at 33, 66, and 100 ms after initial ground contact. Green marks the values in the range of one standard deviation from the mean in in-vitro experiments, blue marks the values within the 1.5 standard deviation range, and orange marks the values outside this range.
Figure 7. Knee joint kinematics and ligaments strain values from FE simulations with optimized HGO materials for ligaments (blue), FE simulations with Neo-Hookean materials for ligaments (orange), and in-vitro experiments (red) with standard deviations error bars for load case #1 at 33, 66, and 100 ms after initial contact. (A) Internal tibial rotation angles, (B) Valgus angles, (C) Anterior tibial translations, (D) Axial compressions, (E) ACL strains, and (F) MCL strains.
Figure 8. Knee joint kinematics and ligaments strain values from FE simulations with optimized HGO materials for ligaments (blue), FE simulations with Neo-Hookean materials for ligaments (orange), and in-vitro experiments (red) with standard deviations error bars for load case #2 at 33, 66, and 100 ms after initial contact. (A) Internal tibial rotation angles, (B) Valgus angles, (C) Anterior tibial translations, (D) Axial compressions, (E) ACL strains, and (F) MCL strains.
Figure 9. Knee joint kinematics and ligaments strain values from FE simulations with optimized HGO materials for ligaments (blue), FE simulations with Neo-Hookean materials for ligaments (orange), and in-vitro experiments (red) with standard deviations error bars for load case #3 at 33, 66, and 100 ms after initial contact. (A) Internal tibial rotation angles, (B) Valgus angles, (C) Anterior tibial translations, (D) Axial compressions, (E) ACL strains, and (F) MCL strains.
Figure 10. Knee joint kinematics and ligaments strain values from FE simulations with optimized HGO materials for ligaments (blue), FE simulations with Neo-Hookean materials for ligaments (orange), and in-vitro experiments (red) with standard deviations error bars for load case #4 at 33, 66, and 100 ms after initial contact. (A) Internal tibial rotation angles, (B) Valgus angles, (C) Anterior tibial translations, (D) Axial compressions, (E) ACL strains, and (F) MCL strains.
Valgus angles and anterior tibial translations were within 1SD of experimental mean data for both materials in all four examined load cases. In total, the simulations with HGO materials led to a closer match with the outputs from in-vitro experiments, with only four data points outside 1.5SD, compared to the simulations with Neo-Hookean materials (22 data points outside 1.5SD).
For ITR angles, most data were in the range of experimental outputs for the simulations with HGO materials, except for load case 1 at 33 ms and load case 2 at 100 ms, which lay outside the 1.5SD and 1SD from the experimental mean, respectively. For the models with Neo-Hookean materials, most ITR angles were in the range of experiments as well; except at 33 ms for load cases 1 & 3, which were outside the 1.5SD range while load cases 2 & 4 at 33 ms, load cases 1 & 2 at 100 ms, and load case 1 at 66 ms time point had results outside 1SD from the mean experiments.
In the case of axial compressions, all outputs at 33 ms time points for both materials in all load cases were in the range of in-vitro outputs. At 66 ms after initial contact, load cases 2 & 3 had outside 1SD data for both materials. At 100 ms, all data were outside the 1SD range for both HGO and Neo-Hook; however, models with Neo-Hookean materials outputs at this time point were also outside the 1.5SD margin for load cases 3 & 4.
None of the ACL and MCL outputs from the simulations with Neo-Hookean materials were in the range of 1SD from the mean experimental values; and for load case #4, the highest loads, all strain values from the models with Neo-Hookean materials were outside the 1.5SD margin (orange cells in Table 6, bottom row).
4. Discussion
Accuracy in modeling soft tissues is of great importance in computational biomechanics [47], and small variations in their properties will result in significant changes in the FE outputs. Knee ligaments are not an exception, and their properties will directly affect the knee joint kinematics and soft and hard tissue stresses and strains. The considerable differences in joint kinematics and ligament strains after changing the hyper-elastic material model of the main knee ligaments in this study proved this point and confirmed the sensitivity of joint biomechanics to these properties in FE simulations.
In this study, the HGO hyper-elastic material parameters for the main knee ligaments were identified via optimization for use in finite element studies. The ACL showed the best match among all the optimized properties, most likely due to a more robust data set available in the literature for ACL mechanical testing. For example, limited studies evaluated mechanical properties of the LCL, specifically no mechanical test data could be found for the human LCL for the young population in the literature. Most studies that reported LCL mechanical properties [48] [49] [50] [51] performed tensile test experiments on elderly cadavers, e.g. (81 ± 11 yrs. old) in [48] (74 ± 7 yrs. old) in [49], (77.1 ± 9.6 yrs. old) in [50], and 71 yrs old in [51]. Therefore, we did not find those appropriate for this study that simulates landing from a jump in young athletes. Also, only a few studies were found for the PCL and MCL. Therefore, getting the accurate properties would be more difficult for these three ligaments relative to ACL. Nevertheless, all the four ligaments’ optimized coefficients presented high correlations (r > 0.95).
Special care should be given to using these properties in other knee FE models if the modeling technique varies substantially, such as not applying viscoelastic effects. Also, it should be considered that coefficients obtained from optimization of material properties are not unique [52] and the purpose of this work was to present a set of data that were validated against in-vitro experiments. Although the proposed coefficients belong to one specific subject, they have been tested in five different subject-specific models in the validation process. Other modelers can apply these properties to help with their progress; they can also use the tools presented here to perform optimizations of their own.
Another important aspect while simulating soft tissue material properties is the conditions of experimental tension tests used. For instance, failure force, failure elongations, and other tensile properties of femur-ACL-tibia complex have been compared for fresh cadaveric specimens versus formalin preserved or deep- frozen preserved, and it has been demonstrated that the preservation method causes significant changes in the ligaments’ mechanical properties [53]. The experiments used in this study utilized fresh frozen samples for MCL [33], −20˚C frozen preserved for PCL [24], and ACL [32]. This study used more current test data (2015) for PCL [24] and from younger specimens for both PCL (30 yr old) [24] and ACL (mean = 29 yr old, 22 - 35 yr) [32] as opposed to the relatively older data [54] used by other studies for ACL and PCL (38 yr old).
The comparison between our optimized HGO coefficients vs. Neo-Hookean coefficients from the literature showed that HGO with optimized coefficients could produce outputs closer to the experimental data; this was expected because of the directionally dependent behavior of ligaments that the Neo-Hookean model does not permit [55], along with the more physiologically relevant coupling of fiber and ground matrix contributions provided in the HGO model. Neo- Hookean material led to larger strains in the FE simulations; this could be because there is less resistance to force in isotropic behavior due to the uniform direction of material fibers. Therefore, Neo-Hookean material showed less resistance leading to greater kinematics at the joint and higher strains in the ligaments.
One limitation of this study was that only one specific geometry for each ligament was used. While these materials properties have shown satisfactory results in five different models, future efforts can increase accuracy by optimizing the properties for each model separately and using subject-specific material coefficients, or providing a statistically relevant aggregate set of coefficients for each ligament. Another limitation concerns the use of uniaxial tensile data for determining the model coefficients where biaxial data may provide more robust coefficients by considering transverse properties of the ligaments in tension. However, as has been shown by other researchers, there is a lack of biaxial test data in the literature for knee ligaments [19]. Lastly, error sources are associated with 3D geometry creation from medical images, such as inaccuracies in segmenting and smoothing, which could affect the results.
5. Conclusion
The coefficients of the HGO material model for the knee cruciate and collateral ligaments were presented in this study using optimization techniques for the use in the knee finite element models while including the viscoelastic properties. When used in finite element simulations of bipedal landing, these coefficients provided closer agreement with the in-vitro experiments for the joint kinematics and ligaments strains relative to the Neo-Hookean materials. Future works may expand upon this study by performing subject-specific optimizations. Also, having access to the tensile test experiment data from a younger population will lead to more accurate optimization results for investigating young athletes’ biomechanics. Therefore, future studies may perform updated tensile test experiments and then redo the optimization to explore the effects on the optimized coefficients.
Acknowledgements
The authors would like to thank Timothy Hewett, Mayo Clinic, Rochester, MN, and his coworkers for providing the CT and MRI images of the subjects and the results of their cadaver study that were used to validate our FE models.
The authors acknowledge partial funding support from the National Institutes of Health (NIH) Grant R01 AR0056259-05A1.