Towards a Transferable Modeling Method of the Knee to Distinguish Between Future Healthy Joints from Osteoarthritic Joints: Data from the Osteoarthritis Initiative

Computational models can be used to predict the onset and progression of knee osteoarthritis. Ensuring the transferability of these approaches among computational frameworks is urgent for their reliability. In this work, we assessed the transferability of a template-based modeling strategy, based on the finite element (FE) method, by implementing it on two different FE softwares and comparing their results and conclusions. For that, we simulated the knee joint cartilage biomechanics of 154 knees using healthy baseline conditions and predicted the degeneration that occurred after 8 years of follow-up. For comparisons, we grouped the knees using their Kellgren–Lawrence grade at the 8-year follow-up time and the simulated volume of cartilage tissue that exceeded age-dependent thresholds of maximum principal stress. We considered the medial compartment of the knee in the FE models and used ABAQUS and FEBio FE softwares for simulations. The two FE softwares detected different volumes of overstressed tissue in corresponding knee samples (p < 0.01). However, both programs correctly distinguished between the joints that remained healthy and those that developed severe osteoarthritis after the follow-up (AUC = 0.73). These results indicate that different software implementations of a template-based modeling method similarly classify future knee osteoarthritis grades, motivating further evaluations using simpler cartilage constitutive models and additional studies on the reproducibility of these modeling strategies. Supplementary Information The online version contains supplementary material available at 10.1007/s10439-023-03252-8.


Introduction
Articular cartilage is a fibril-reinforced biphasic tissue that offers frictionless contact and impact absorption functions in the knee joint for many years of an individual's life. Unfortunately, a variety of risk factors such as aging, joint injury, and being overweight [3] may impair cartilage optimal functioning, leading to a degenerative condition called osteoarthritis (OA). Therefore, understanding OA driving mechanisms is key to developing treatment strategies for this painful condition [11].
Currently, there is no clinical method for personalized risk estimations for the onset and progression of knee OA. In this regard, computer-based approaches have been recently utilized to simulate the biomechanics of the knee in healthy and diseased conditions [27]. The results of those simulations, e.g., stress and strain distributions, can be further utilized to optimize surgical procedures or preventive actions [12]. For instance, by identifying loadings that cause abnormal cartilage hip mechanics [13] or minimizing bone stress concentrations after ankle fracture reductions [18]. However, regarding conservative measures, such as weight loss or gait retraining, validated computational approaches are missing to evaluate the personalized effects of those measures with simulated disease progression.
Modeling the biomechanics of living tissues in the human body is challenging due to the complex geometries and loads they bear. For the knee joint, several strategies based on finite element analysis (FEA) have been developed, accounting for the complex tissue structure, subject-specific geometries, and loading conditions of the joint [2,8,24,30]. In addition, when the studies involve cohort data, generating and running a large number (N > 100) of personalized models are time-consuming processes [32]. In this regard, Mononen et al. (2019) proposed a template-based approach to rapidly simulate the biomechanics of the medial compartment of the knee using FEA. Yet, they modeled only 21 subjects, which could represent a limitation for the reliability of the approach. Moreover, they used a user-defined cartilage formulation in ABAQUS, which is a commercial software not readily available to many users, possibly limiting the use of the method by researchers worldwide. Approaches like this should be replicated with a large patient population and by using more accessible tools, including freely available FE software like FEBio [25], which is receiving increasing attention to model multiphasic problems in biomechanics due to its flexibility given by a large library of non-linear constitutive equations that are easily implementable. More accessible methods would enable faster development of novel strategies that could help in clinical decision-making on subject-specific conservative measures. For instance, by visualizing abnormally loaded regions in the joint and simulating preventive interventions to alleviate those harmful conditions, as well as by offering means to quantify the risk for OA.
In this study, we extended the template-based method [26] to over 100 joints to evaluate its capability for predicting the development of knee osteoarthritis. We generated the models and reproduced the results by ABAQUS and FEBio finite element softwares by using baseline and follow-up information from subjects from the Osteoarthritis Initiative (https:// nda. nih. gov/ oai/). We hypothesize that the template-based FEA pipeline, implemented in two different finite element softwares, distinguishes subjects at high risk for knee OA from those at low risk, even though all subjects have similar demographic characteristics and healthy radiographic knee conditions at baseline. Figure 1 overviews the workflow of this study. First, medial knee compartment FE models were generated using the template-based approach to simulate the stance phase of the gait [26]. For that, we used anatomical measurements from the distal femur and tibiofemoral joint space obtained from clinical MRI, and demographic information (age and weight) of the subjects at a common healthy baseline. Then, we used age-dependent thresholds of maximum principal stresses, based on reported experimental observations of monotonic tensile tests of human cartilage samples [17,26], to compute the volume of degenerated cartilage, assuming that the degeneration initiates in the collagen network [14,21]. Finally, we compared degenerated volumes between healthy and OA knees, classified based on radiographic evidence at the 8th year of follow-up, using two different finite element (FE) softwares: FEBio [25]

Subjects Sample
We obtained baseline and 8-year follow-up information from 78 subjects from the Osteoarthritis Initiative database. This number of subjects corresponds to 154 knees since the information from both knees was not available for all the subjects. Figure 2 shows the inclusion criteria for the subject selection. The data used in the study included MRI images and Kellgren-Lawrence (KL) scores for the left and right knees, as well as the subjects' weight and age.
Anatomical dimensions of the distal femur and tibiofemoral joint space were measured for all the knees at the healthy baseline (KL grade 0 or 1) and used to scale a template finite element model (Fig. 1a) [26]. The dimensions correspond to the medial-lateral maximum condylar distance (M-L), the maximum anterior-posterior distance (A-P) measured in the sagittal plane of the medial condyle, and the joint space (JS) measured in the same plane as the A-P parameter. The knees were grouped based on their KL status at the 8-year follow-up. We will refer to individuals with KL = 0 (KL0) and KL = 1 (KL1) as healthy, KL = 2 (KL2) as mildly affected, and KL = 3 and KL = 4 (KL34) as severely affected by OA. Table 1 summarizes the subject characteristics and measurements (a visual representation of the information can be found in Supplementary Material, Section 1).

Finite Element Models
We made assumptions and simplifications aiming for the rapid simulation of subject-specific knee models under a gait loading condition. Regarding the geometries, FE models comprised the cartilages of the medial compartment of the tibiofemoral joint [26] since knee OA is more frequent in the medial than in the lateral compartment [3,19]. We did not include ligaments, tendons, and muscles to avoid increasing computational costs [7]. However, their constraining effects were considered through appropriate but simplified boundary conditions, by fixing anterior-posterior and medial-lateral translations, and internal-external rotations. The meniscus was excluded from the simulations, but its contribution was accounted for indirectly by subtracting the meniscus support force from the force acting through the tibiofemoral joint [26]. This assumption 1 3 conveys similar stress and fluid pressure distributions in cartilage-cartilage contact regions during walking as if we would physically consider the meniscus in the model [26]. Regarding loading, we prescribed generic effective joint axial force and knee flexion over an analytical point located at the mid-distance between the epicondyles of the distal femur. We assumed the medial compartment of the knee bears half of the joint reaction force during gait. This is because we did not account for the individual motion data of the subjects, and this assumption fits into the experimental observations for the medial/lateral load sharing [22,43].

Template Model Scaling
We followed the method proposed by Mononen et al. (2019). We scaled the FE mesh of one existing template based on the fractional differences in M-L, A-P, and JS parameters between the subject of interest and the template (Fig. 1). A cartesian scaling was used for anatomical dimensions, except for femoral cartilage thickness, where a cylindrical scaling was utilized. Regarding loading, we scaled the axial joint contact force [23] based on the body weight of each subject and kept the knee flexion angle trajectory [6] the same for all models.

Framework Comparisons (ABAQUS vs FEBio)
Finite element software helps model complex problems involving multiphasic physics in biomechanics. ABAQUS (Dassault Systèmes Simulia Corp., Providence, RI) is a wellknown commercial software in computational biomechanics with useful characteristics for researchers, like user-defined functions for material formulations and boundary conditions. On the other hand, FEBio is a freely available option for solving these non-linear biphasic problems [25]. In this regard, we aimed at exploring the transferability of the template-based approach using both of these highly known programs.
In a previous study, we showed how to impose equivalent boundary conditions and similar constitutive equations in a knee joint model using both ABAQUS and FEBio softwares [33]. Regarding cartilage material, in ABAQUS standard v.2018, we used an experimentally validated fibril-reinforced poroviscoelastic formulation [15,41,42], while in FEBio 3.3.2 we used a fibril-reinforced biphasic description [5,33], with material properties obtained to reproduce the instantaneous response of the material used in ABAQUS. In both frameworks, we implemented the same depth-dependent orientation for collagen fibrils. Table 2 summarizes the cartilage model parameters.
We compared the results of both finite element softwares in terms of time and space. These comparisons include average and peak values over the contact area, stress maps at different fractions of the stance, and volumes of tissue overstressed for different KL grade groups. The "Statistical Analysis" section describes these comparisons.

Sensitivity Analysis
We performed a factorial design to assess the sensitivity of the peak tensile stress occurring at the first load peak of the stance phase of gait with respect to the parameters used to generate the compartment models. This is because our degeneration model relies on age-dependent thresholds of cartilage tensile stress to define the volume of tissue at risk (see "Damage Model" section). The four factors were the weight of the subject and measurements of knee M-L distance, medial knee A-P distance, and medial JS width. The low and high levels of these factors were the minimum and maximum values observed in the sample, respectively. Table 3 summarizes the 16 combinations of the four factors and the two levels, e.g., model 1 was constructed using the high level for all parameters. Pain frequency > 2 indicates daily to constant knee pain, and pain severity > 5 refers to scores in the upper half of a self-reported pain level from 0 ("no pain") to 10 ("pain as bad as you can imagine"). The final number of subjects was 78 (154 knee models). Namely, 88 knees with KL grade 0, 15 with KL1, 36 with KL2, and 15 accounting for knees with KL3 and KL4 To assess the sensitivity of the peak stress with respect to each factor, we first calculated the overall mean peak stress from the 16 models ( Peak ). Then, we calculated the mean peak stresses from the eight models sharing either the highest ( + ) or the lowest (−) level i of each factor j , e.g., Peak +,j or Peak −,j for the high and low levels, respectively. Finally, we computed the relative difference of each of these peaks with respect to the overall mean peak stress,

Damage Model
We used tensile stress thresholds to define tissue degeneration. In this way, we assume the damage is initiated in the fibril collagen network since this constituent is mainly in charge of the tensile stresses [14,21]. We considered the threshold for damage ( T f ) according to experimental observations of tensile tests of human articular cartilage over a wide range of ages [17,26]. Then, T f was defined as, (2a)  Weight Subsequently, we hypothesized that larger volumes of overstressed tissue positively correlate with a higher risk of developing knee OA. Thus, in our analyses, the volume of tissue at risk was equal to the volume occupied by the elements that exceeded the thresholds during the simulated loading cycle.

Classification Capabilities
We tested our method's ability to distinguish future osteoarthritic knees from healthy knees using simulated degenerated tissue volumes. To this end, we modeled the knees at the healthy baseline (KL0,1) and computed the volumes at risk. We then grouped the volumes using the KL grades determined in the 8th year of follow-up and reported the area under the curve (AUC) of receiver operating characteristic (ROC) curves between pairs of KL grades.

Statistical Analysis
We compared the results of two finite element softwares as well as the volumes of overstressed cartilage from the different KL groups of simulated knees.
To compare the stress results of both FE softwares, we used one-dimensional (1D) and two-dimensional (2D) statistical parametric mapping (SPM{t}). [31] In 1D statistical analyses, we compared the average and peak values of the maximum principal stress over the tibial contact area during the load cycle. We used the two-tailed 1D nonparametric SPM{t} (α = 0.05) since these variables did not follow normal distributions. In 2D statistical analyses, we compared the spatial stress distributions in the superficial cartilage layer for different time points of the stance fraction. For that, we staked the histories of stress maps from FEBio and ABAQUS models and used the two-tailed 2D nonparametric SPM{t} (α = 0.05). This was possible because we used the same meshes and matched the time points between the softwares.
Regarding the comparison of the volumes of overstressed tissue from each software for the same KL groups, we used the Wilcoxon signed-rank test. To compare the (2b)

Results
The temporal and spatial distributions of the maximum principal stress showed some variation that was dependent on the software that was utilized. While the principal stress averaged over the contact area was essentially the same in both programs (Fig. 3), the stress peaks during the load cycle using ABAQUS were on average 27% higher than those from FEBio (Fig. 3). Comparisons for maximum principal strain and fluid pressure can be found in Section 2a of the Supplementary Material. Regarding the spatial distribution, Fig. 4 highlights the regions where the maximum principal stress statistically differed between softwares at different stance fractions. The red zones indicate that ABAQUS offered higher values compared to FEBio, and the blue zones indicate the opposite trend. The gray color depicts regions with no differences. In Section 2b of the Supplementary Material, the reader can appreciate the influence of the scaling geometry parameters (M-L, A-P, JS) on the differences between programs.
Regarding the sensitivity analysis, Fig. 5 shows higher peak tensile stresses in ABAQUS models compared to FEBio models. Both FE software implementations of the template-based approach yielded similar relative sensitivity of stresses to the parameters used to generate the models. The weight and joint space (JS) were the factors that had the greatest influence on the peak stress, with the steepest slopes, as shown in Fig. 5. In contrast, the anterior-posterior distance (A-P) caused a minimal variation in the peak stress between the low and high levels.
Regarding the simulated degenerated volumes, for the corresponding KL groups, the ABAQUS volumes were about 43% greater than those from FEBio (Fig. 6). Nevertheless, both softwares yielded significant differences between healthy (KL0,1) and severely affected (KL3,4) knees, when the knees were separated into three or four groups (Fig. 6).
Finally, looking at the layers of the overstressed tissue (Fig. 7), the superficial zone is responsible for ~ 80% of these volumes, while the middle zone contributes ~ 15% and the bottom zone contributes ~ 5%. These proportions remained constant between KL groups, independently of the software utilized.

Discussion
A template-based finite element approach was implemented in commercial (ABAQUS) and freely available (FEBio) FE software to generate knee joint models of 154 healthy knees. We used the simulated levels of maximum principal stresses in the knees to predict personalized cartilage degeneration. The predicted volumes of degenerated tissue were compared against KL grade changes after an 8-year follow-up. Though there were notable differences in the simulated absolute values for the maximum principal stresses, similar OA classification outcomes were predicted for the healthy and osteoarthritic knee joints with both softwares based on the ROC analysis. This suggests that simpler constitutive models for cartilage may be  Volume comparisons and ROC curves between KL grades using ABAQUS and FEBio frameworks. a Classification using 3 KL grades and b using 4 KL grades. Statistical significance was evalu-ated using Kruskal-Wallis and Dunn's tests between KL grades, and the Wilcoxon rank test for the same KL grades. *p < 0.05, **p < 0.01 used to simulate OA degeneration. In addition, we investigated the contributions of subject-specific factors (joint shape and weight) on the simulated peak stresses in cartilage, as this is utilized to generate a prediction for cartilage degeneration as a function of age. From this analysis, we found that peak stresses in both software implementations were similarly sensitive to the parameters used to generate the models.
Results showed notable differences in the tensile stress over the tibiofemoral contact area during the gait loading between softwares. The stress distributions showed regions where either ABAQUS or FEBio simulated higher values. These differences may arise from the different material formulations between ABAQUS and FEBio. In ABAQUS, we implemented an FRPVE cartilage model that we mimicked with a fibril-reinforced poroelastic cartilage model in FEBio [33]. Therefore, the variability in loading conditions with respect to those used to calibrate the FEBio material properties exposes this difference in formulations.
The sensitivity analysis indicates that the peak of maximum principal stress was similarly sensitive to the anatomical and loading factors used in the template-based FE modeling in both softwares. This finding supports the idea that the method could work as a means for studying these biomechanical risk factors and translating them into quantifiable modeled variables such as stress or volumes of tissue at risk for the onset of OA development. Weight positively influenced peak stress, indicating that being overweight increases the chances of the tissue failing, which agrees with known biomechanical knee OA risks [3]. On the other hand, the JS negatively influences the peak stress, suggesting that thinning of tissue also increases the stress when it bears the same weight compared to a control sample, which is in agreement with previous studies correlating biomechanics and tissue adaptation or damage [4,29]. The medial-lateral (M-L) and anterior-posterior (A-P) distances showed a lower influence compared to weight and joint space. However, we considered that M-L and A-P dimensions may affect the congruency of the joint and may have a greater influence if the loading involves larger relative rotations of bone ends, like in deep squatting [22].
Regarding the predictions of cartilage degeneration, we found that the simulated volumes of overstressed tissue differed between software implementations. Although, both cases acceptably classified healthy (KL01) from severely affected joints (KL34), suggested by an AUC ~ 0.73. In addition, despite the promising difference observed between KL1 and KL3,4, using 4 KL groups, the proposed approach did not show differences between KL0 and KL1 groups. It should be noted that KL1 does not indicate degenerative changes in cartilage, but rather doubtful joint space narrowing (and possible osteophytic lipping) [20], which can be interpreted as normal joint space in many cases. Therefore, it might be an unrealistic goal to develop a method that could predict future joint conditions between KL0 and KL1.
We based our predictions on the short-term biomechanical response of the collagen network, simulated by the maximum principal stresses (tensile stresses). However, other failure criteria and loading conditions should be explored for predicting the onset and progression of cartilage damage with cohort data. For instance, the fatigue behavior of cartilage under tension [39,40], shear, compression [16], or sliding [9], as well as tissue remodeling induced by the biochemical signaling of cells [28,36]. Combining these damage mechanisms may yield more accurate physics-based computational approaches to determine not only the severity but the location and pathway of tissue damage. Although, validating such models is challenging since it is hard to obtain suitable experimental measurement data to compare with, especially when using knee-level models.  Figure 7 revealed that most of the volume of overstressed tissue is located in the superficial layer of cartilage. This distribution is explained by the fibril-reinforced formulation and the depth-dependent orientation of the collagen network. In addition, this finding suggests that, for predictive purposes, focusing our attention only on the cartilage surface may represent a saving in computational time. In this regard, using contact stress from discrete element analysis (DEA) [1,38] seems to be a possibility. DEA approaches have been evaluated in predicting the future stages of knee OA [35] and ankle post-traumatic OA [18] with promising results. Although we acknowledge that this method is not appropriate to study structure-function alterations in the tissue.
The time required for simulating the 154 knee models was 34 h in FEBio and 123 h in ABAQUS. This difference in time is explained by the differences in material model formulations between ABAQUS and FEBio, and the limited parallelization allowed by the licensing of ABAQUS. Although the original cartilage degeneration algorithm was developed with ABAQUS using the FRPVE material formulation for cartilage [26], the use of FEBio is also recommended based on the similar classification accuracy (KL01 vs KL34 knees) we observed. In FEBio, users can run as many models as computing power allows, while in ABAQUS, users are limited by the number of licenses (tokens) purchased. The benefit of using FEBio also offers a more cost-effective route to prediction results with large datasets.
The limitations of the present work include the absence of the lateral tibiofemoral compartment, restricting the identification of regions at risk to the medial compartment. In addition, we assumed a fixed load distribution between medial and lateral compartments, since the Osteoarthritis Initiative database does not have motion data. We assumed a scaled generic knee loading and a fixed flexion trajectory, limiting the study of the effects of subject-specific gait patterns that may overload or protect different contact regions in the knee. We did not use subject-specific material properties in the joint, which may show the sitespecific onset of OA [24]. Another characteristic that we consider a limitation was not including local morphologies in the cartilage or subchondral bone that can induce the onset and progression of damage [37], making the scaling of the template geometry sound simplistic. However, in the original paper by Mononen et at. 2019, they compared 21 knee models using the actual subject-specific geometries, obtained from manual segmentations, and scaling templates. They found that the model with the scaled geometries produced better results. Therefore, we did not consider other solutions to implement template scaling. We additionally consider that the low number of OA knees can be a limitation if we aim to calibrate predictive models using, for instance, machine learning tools. Furthermore, we compared only two softwares to give insights into the transferability of the method. In this case, it would be worth exploring the reproducibility not only between softwares but also among researchers. [10,34] In the future, we will implement our template-based approach to consider the lateral compartment of the tibiofemoral joint. This improvement should include the geometry of the lateral compartment and the differentiated medial/lateral loading during the stance phase of gait. Moreover, we will explore the capability of other biomechanical responses, such as the maximum shear strain, to identify the future knee OA stage, as they have been associated with damage in other tissue constituents than the collagen network [14]. Finally, considering local defects in the models can be a valuable feature to expand the use of our approach in cohort studies of post-traumatic knee OA. This could be done by rapidly defining affected zones by imposing adequate boundary conditions and depleting the mechanical properties of the tissue (e.g., low elastic properties, high permeability), as well as by performing proper local mesh refinement.
In conclusion, the template-based method to simulate the biomechanical behavior of the medial compartment of the knee was extrapolated to two widely used FE softwares obtaining similar acceptable classifications of the future knee KL grades. Through this work, we encourage the development of reproducible methods to study and predict the degeneration of human body structures, aiming for tools to help clinicians elucidate conservative treatments for multifactorial musculoskeletal disorders.

Supplementary Information
The online version contains supplementary material available at https:// doi. org/ 10. 1007/ s10439-023-03252-8. no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http:// creat iveco mmons. org/ licen ses/ by/4. 0/.