Mechanical Metrics of the Proximal Tibia are Precise and Differentiate Osteoarthritic and Normal Knees: A Finite Element Study

Our objective was to identify precise mechanical metrics of the proximal tibia which differentiated OA and normal knees. We developed subject-specific FE models for 14 participants (7 OA, 7 normal) who were imaged three times each for assessing precision (repeatability). We assessed various mechanical metrics (minimum principal and von Mises stress and strain as well as structural stiffness) across the proximal tibia for each subject. In vivo precision of these mechanical metrics was assessed using CV%RMS. We performed parametric and non-parametric statistical analyses and determined Cohen’s d effect sizes to explore differences between OA and normal knees. For all FE-based mechanical metrics, average CV%RMS was less than 6%. Minimum principal stress was, on average, 75% higher in OA versus normal knees while minimum principal strain values did not differ. No difference was observed in structural stiffness. FE modeling could precisely quantify and differentiate mechanical metrics variations in normal and OA knees, in vivo. This study suggests that bone stress patterns may be important for understanding OA pathogenesis at the knee.

properties, cyst formation, and alignment changed stress/strain distributions and structural stiffness of bone. Although the previous research has suggested a potential structural role of bone in OA initiation and development, most FE research-to-date has been simulation-based with a focus on individual alterations (e.g., studying the effect of cyst size on the stress levels by McErlain et al. 22 and studying the effects of subchondral bone and cartilage stiffness on the stress levels by Dar and Aspden 24 ). To examine the combined effects of these individual bone alterations (i.e., bone mineral density (BMD) distribution, alignment, presence of cysts, etc.) and represent accurate bone geometry and material properties, we need to develop subject-specific FE models that can incorporate various bone alterations. Quantitative computed tomography images (QCT) are commonly used to create subject-specific FE models 22,[25][26][27] . In this method, the geometry of the bony tissue is obtained using segmented QCT images and image-based densities are converted to elastic moduli using published density-modulus relationships. With subject-specific FE modeling, we can assess various mechanical metrics of bone (which cannot be measured otherwise), and study them in OA and normal knees. However, information regarding precision (i.e., repeatability) errors for FE-based mechanical metrics is important when applying these methods to analyze and compare OA and normal knees.
Our objective was to identify precise mechanical metrics of the proximal tibia that differentiated OA and normal knees. For this, we used QCT scans of the knees of fourteen participants. Seven of the knees were classified as normal while the remaining seven were classified as having OA (details of participant and image acquisition is provided in the Methods section). Each knee was imaged three times for the purpose of calculating precision errors. Subject-specific FE models of the knees were created, and single-leg stance loading was simulated (see details in Method section). FE-based minimum principal and von Mises stress and strain were obtained, as well as medial and lateral proximal tibial stiffness. Regional stress and strain were assessed across the proximal tibia ( Fig. 1). Short-term precision errors were assessed by calculating root mean square coefficients of variation (CV% RMS ) for all outcomes and in different regions of the proximal tibia. As well, pairwise comparisons of FE-based outcomes as well as BMD, bone mineral content (BMC) and bone volume were performed between OA and normal knees.

Results
Precision Study. CV% RMS for both OA and normal knees ranged from 3.7% to 10.5% (average: 6.1%) for minimum principal stress (Table 1) and 3.2% to 7.6% (average: 5.5%) for minimum principal strain ( Table 2). Similar precision errors were noted for von-Mises stress (Supplement Table 1) and strain (Supplement Table 2). Precision error for structural stiffness for both OA and normal participants was 3.6% in the medial compartment and 5.0% in the lateral compartment (Table 3).

Preliminary Comparisons of OA and Normal FE Outcomes.
Overall, stress (after scaling for bodyweight) was higher in OA versus normal bone. The percent difference in minimum principal stress was 115% in the medial peripheral cortical bone (~20x CV% RMS ) and 107% in the medial epiphyseal cortical region (~16x CV% RMS ) (Table 1). Similarly, von-Mises stress values were higher in OA knees compared to the normal knees (Supplement Table 1). Mapped representations of minimum principal stress results onto CT images of the proximal tibia are shown in Fig. 2. Qualitative analyses of the stress patterns indicated that the load was primarily transferred through the medial compartment for individuals with OA (Fig. 2).
Regarding the strain values, no difference was observed in minimum principal or von-Mises strain values between OA and normal knees (p > 0.05) ( Table 2 and Supplement Table 2). The highest strain was observed in the epiphyseal central region for both normal and OA proximal tibia. Similar to the strain results, no differences were observed in the structural stiffness of medial and lateral condyles of proximal tibia between OA and normal  (Table 3). When we compared BMD, BMC and bone volume values between OA and normal knees (Supplement Tables 3-5), no differences were noted, except BMD for the medial peripheral cortical region (p = 0.039) (Supplement Table 3).

Discussion
This is the first study to report in vivo precision errors for subject-specific FE-based mechanical metrics for various regions of the proximal tibia. This is also one of few studies to differentiate OA and normal knees using FE-derived estimates of mechanical behavior. The developed FE method provided precise estimates of von-Mises stress, minimum principal stress, von-Mises strain, and minimum principal strain in the different regions of the proximal tibia as well as structural stiffness of the medial and lateral compartments (average CV% RMS < 6%). Our precision errors were in the same range as previously reported CV% RMS precision errors for FE modeling using high-resolution peripheral QCT imaging (HR-pQCT) to measure bone stiffness and stress [28][29][30] .
On average, minimum principal stress was 75% larger in the medial side of the OA bone compared to normal bone. The difference in stress was as high as 20 times the precision error for minimum principal stress, indicating that the difference is far greater than repeatability error. These findings, combined with large effect sizes (Cohen's d) suggest that FE-based mechanical metrics can differentiate OA and normal knees. On a related note, we also calculated the in vivo precision of regional BMD, BMC and bone volume measures of the proximal tibia and when we compared these values between OA and normal knees, no differences were noted (except for BMD for the medial peripheral cortical region). This further emphasizes the importance of exploring the combined effects of individual OA-related bone alterations.
The obtained mechanical metrics are in agreement with previously published FE analyses 22,26,31 . The reported von-Mises stress values were within the same range (around 0.9 MPa with one body weight for the medial subchondral cortical region) as previous subject-specific FE modeling of the knee joint by McErlain et al. 22 and an FE study of the proximal tibia by Tuncer et al. 26 with similar loading conditions. Also, the obtained strain values in this research (around 700 microstrain with one body weight for the medial subchondral trabecular region) compared favorably with the strain range reported for FE modeling of the proximal tibia 26 31 . Although speculative, it is worthwhile hypothesizing on why stress was higher in OA knees as this may provide a better understanding of the mechanical behavior of bone in response to OA. We believe stress was higher due to the combined effects of various factors, including: (1) higher weight; (2) malalignment; and (3) BMD difference between OA and normal knees.
First, OA participants in this study had significantly higher weight compared to the healthy participants (+40%, p-value < 0.05), which could contribute to the observed higher stress in OA knees. Though, it is important to note that the difference in weight of OA and normal groups is much lower than the observed differences in stress values (up to +115%). Second, with regards to alignment, OA knees were slightly varus which appeared to contribute to higher stress values in the medial side of the OA proximal tibia. Third, although the differences in BMD values were not significant, there was a trend for higher BMD with OA proximal tibiae (Supplement Table 3). Because Goulet's equation for material mapping relates BMD to elastic modulus by a power of 2.1 32 , a slightly higher BMD may lead to significant differences in elastic modulus. Importantly, higher elastic moduli will change how load is distributed through the tibia, with highest elastic moduli carrying the most load. Accordingly, higher elastic moduli will lead to higher local stress. Related to this point, simply adjusting all density measures  Table 2. Minimum principal strain comparison between OA and normal proximal tibia. Mean and SD of repeated scans for both OA and normal, CV% RMS , mean and SD for OA knees, mean and SD for normal knees, the difference between OA and normal knees (absolute and percent relative to normal), 95% confidence of interval, p-value, and effect size (Cohen's d) of minimum principal strain in different regions of proximal tibia. *Shows regions which were not normally distributed whereby median value used in Mann-Whitney U-tests for statistical comparison, and confidence intervals were calculated using Hodges-Lehmann estimator.  Table 3. Medial and lateral stiffness comparison between OA and normal proximal tibia. Mean and SD of repeated scans for both OA and normal, CV% RMS, mean and SD for OA knees, mean and SD for normal knees, the difference between OA and normal knees (absolute and percent relative to normal), 95% confidence of interval, p-value, and effect size (Cohen's d) of structural stiffness in medial and lateral compartments.  Table 6). It is important to be cognizant that, with subject-specific FE modeling, individual factors such as higher weight, altered alignment and altered BMD are inherently integrated, thereby providing insight regarding their combined effect. Statistical approaches could be used to account for these factors to derive surrogate estimates of stress; however, they could be erroneous due to the employed assumptions (e.g., independency of observations, linear vs. non-linear associations, distribution of residuals). Also, statistical approaches typically require large datasets whereas subject-specific FE modeling can be employed at the individual-level. Interestingly, while stress was different between two groups, no difference was observed in the strain values. Similar strain levels could be an indication of bone adaptation in response to altered loading in the OA joint 33 . Further research with individuals at early stages of OA is needed to assess whether this is a progression-specific phenomenon.
This study had certain strengths worthy of discussion. First, this study met the conservative number of the patients and repeated scans per patient, as proposed by Gluer et al., which allowed us to establish reliable precision errors with an upper 90% confidence limit having less than 30% error (e.g., if the precision error was 2%, then we were 90% confident that the true precision error is less than 2.6%) 34 . Second, we assessed in vivo precision errors in order to incorporate uncertainty associated with subject motion. This approach offers a more realistic estimate of repeatability error than an ex vivo precision study of motionless cadavers or repeated analyses of a single scan. Third, we used CT images of study participants to create FE models as opposed to idealized geometries that resemble the knee joint. Fourth, the heterogeneity of bone mechanical properties was considered using a density-modulus equation which has been shown to affect the accuracy of QCT-based FE models 25,35 . Fifth, we aligned all knees in similar 3D orientations relative to landmark boundary points and best-fit planes. This minimized differences due to dissimilar orientations during image acquisition, thereby permitting dependable comparisons between OA and normal FE outcomes. Sixth, we provided precision errors for each mechanical metric as well as differences between OA and normal knees. As these differences were, on average, 14 times greater than precision errors, these differences are trustworthy.
Limitations of this study pertain to small sample size for statistical comparison, the use of a single density-modulus equation for both cortical and trabecular bone, and simplified modeling of cartilage and menisci. First, although the number of study participants was sufficient for calculating precision errors, we had low statistical power for detecting differences between OA and normal knees. As this was a preliminary study to detect potential differences between knees, we believe the use of a small sample size for statistical comparisons was justified. Further analyses with a larger sample size are needed to corroborate study findings. Second, a single density-modulus equation was used to derive elastic moduli for both cortical and trabecular bone in the FE models. While applying cortical-specific and trabecular-specific density-modulus relationships offers improved predictions of stiffness (R 2 = 0.77 versus R 2 = 0.70 with a single equation), these improvements are moderate. As such, we believe a single equation is justified to provide reasonable estimates of proximal tibial mechanical behavior 25,36 . Related to this point, we did run a pilot study with cortical-specific and trabecular-specific equations for 3 samples and FE outcomes were similar to those reported with a single equation. For this reason, we applied the simpler single-equation approach for converting BMD to E. Third, as only CT images were available, realistic subject-specific modeling of cartilage and menisci was not feasible. To overcome this limitation, a cylinder of tissue with an isotropic homogeneous material was used to model soft tissue of the knee joint, as per McErlain et al. 22 , with the same material properties for both OA and normal knees. Our approach to model soft tissues has both negative and positive aspects. As participants were not bearing weight during scanning, our model of soft tissue may not mimic the relative thickness of medial versus lateral cartilage, which could lead to erroneous estimates of how load and stress is distributed through the proximal tibia. On the other hand, without weight bearing our model likely overestimates cartilage thickness. This is important because prior research indicates that thinner cartilage leads to higher bone stress 24 . The approach employed for modeling cartilage also ignores cartilage defects seen with OA 37 , which would lead to a smaller contact area and higher bone stress in the subchondral cortical region. Accordingly, we believe this approach makes our analysis conservative. In other words, if subject-specific cartilage geometry and material properties were used, we would reasonably expect to see higher differences in stress between OA and normal bone due to thinner cartilage and the smaller contact area between articulating bones. Also, incorporation of estimated cartilage and meniscal geometry would have added new variability to the model. As this analysis was focused on bone, we aimed to create a reasonable, functional in vivo model of the knee which could be used to investigate the structural role of bone in OA. We believe we have achieved that aim; though, further research and development are needed to incorporate cartilage and meniscal structures with subject-specific geometry and mechanical properties (most likely via registration with magnetic resonance (MR) images) and lessen model development time (each model takes approximately 3 hours to develop and process).
Results of this study indicate that subject-specific FE modeling has potential to precisely quantify and differentiate mechanical metrics variations in normal and OA knees, in vivo. Preliminary findings suggest that OA bone exhibits higher stress levels compared to normal bone. In contrast, strain and structural stiffness did not differentiate between OA and normal bone. Subject-specific FE modeling may reveal insight into the structural role of bone in OA pathogenesis and warrant application in future studies.

Methods
Study Participants and OA status. In this study, we used CT images of 14 participants who were recruited for a previous study (3 men, 11 women, mean age 49.9, SD 11.9 years) 38 . Previously, OA status of the knees was assessed using a modified Kellgren-Lawrence (KL) OA severity scoring system 39  QCT Acquisition. Knee selection for CT imaging and imaging parameters were described previously 38 . If the participants had knee pain, the more painful knee was selected for imaging; otherwise, a random selection of left or right knee was scanned. Participants were scanned three times over two consecutive days. The knee of interest for each participant was imaged via single-energy QCT using a clinical CT scanner (Lightspeed 4-slice, General Electric, Milwaukee, WI, USA). Scanning was performed in the supine position of participants while the knee of interest was centered within the CT gantry. A solid QCT reference spine phantom (Model 3 T; Mindways Software Inc, Austin, TX, USA) was included in the images in order to convert grayscale CT Hounsfield units (HU) to equivalent apparent volumetric BMD (mg/cm 3 K 2 HPO 4 ). Scanned image volumes contained the distal femur, patella, proximal tibia, and fibula; though, image volumes were cropped to exclude the patella for this analysis. CT scanning parameters included: 120 kVp tube voltage, 150 mAs current-time product, axial scanning plane, 0.625 mm isotropic voxel size (0.625 mm slice thickness, 0.625 × 0.625 mm in-plane pixel size), ~240 slices, ~90 s scan time. Edge enhancement and post-processing were done using a standard bone reconstruction kernel (BONE). Effective radiation dose was estimated at ~0.073 mSv per scan using shareware software (CT-DOSE, National Board of Health, Herley, Denmark). This value is comparable to the average effective radiation dose during a transatlantic flight from Europe to North America (~0.05 mSv).
CT Image Analysis. We separated the proximal tibia, distal femur, and fibula from the surrounding soft tissue in the QCT images using semi-automatic segmentation and manual corrections. A subject-specific bone threshold, obtained using the half maximum height (HMH) technique, was used for segmenting each image 40,41 . In the HMH method, the density of a voxel with 50% cortical bone and 50% joint space is used as the minimum threshold for segmenting subchondral bone. We performed the segmentation using a commercial image processing software (ANALYZE10, Mayo Foundation, Rochester, MN, USA) (Fig. 3a,b), a stylus, and an interactive touch-screen tablet (Cintiq 21uX, Wacom, Krefeld, Germany). FE Modeling. Geometry. As CT imaging was performed in the supine position, the knees were reoriented in a neutral standing alignment (MATLAB, MathWorks, Natick, MA, USA) to model single-leg stance.
Two vectors were defined for re-alignment; one vector was a best-fit line passing through the centroid of tibial cross-sections at different levels of the tibia (from distal to proximal: ankle; midshaft (50% length of tibia); 66% shaft site; proximal tibia). The other vector was a best-fit line passing through 32 femoral cross-sections starting just proximal to the intercondylar fossa and extending 2 cm proximally towards the center of the femoral head. The images were rotated such that the average of the two vectors was vertical. With this approach, both the femoral axis and the tibial axis make the same angle with the vertical axis (e.g., the varus angle was 175 °, the angle between the femoral axis and the vertical axis would be 2.5 °, and the same for the angle between the tibial axis and the vertical axis). Visually, the new alignment of the CT images was similar to generic standing MR images (Fig. 4).
We converted the segmented, re-aligned dataset to a 3D polygonal surface mesh using a marching cube algorithm (ANALYZE10, Mayo Foundation, Rochester, MN, USA) and imported the obtained 3D object into reverse engineering software (GEOMAGIC STUDIO 12, Systems, Rock Hill, SC, USA). We smoothed the surface of bones to ensure the geometry was topologically valid and did not contain holes or rough edges (Fig. 3c). To maintain geometric complexity, the maximum deviation of the smoothed surface from the original meshed surface  was less than one voxel size (0.625 mm). We imported the resultant smoothed 3D volume of the knee into FE software (ABAQUS, Providence, RI, USA). Because soft tissues (e.g., cartilage, menisci) were indistinguishable in CT images, these tissues were modeled by an incompressible cylindrical medium, a method previously employed by McErlain et al. 22 . We then meshed all structures in the FE model using 10-noded, quadratic, tetrahedral elements (generic element size of 2 mm for bone and a generic element size of 20 mm for the soft tissue cylinder) (Fig. 3d). To identify an appropriate element size, we performed a mesh convergence study on seven different knee models. With this, FE results of stress and strain did not change by more than 1% when changing the element size from 2 mm to 1.8 mm. In total, FE models consisted of ~500000 to ~700000 elements depending on the size of the knees. Bonded contact was used between bone (the femur, tibia and fibula) and the soft tissue cylinder interface.
Material Properties. We converted imaged BMD to elastic modulus using a density-modulus equation proposed by Goulet et al. 32 (Supplement Table 7). This equation has been shown to explain 70% of variance in proximal tibial subchondral bone stiffness 25 , and two studies recommended its use for modeling the proximal tibia 25,36 . We used a custom algorithm (MATLAB), reported by Nazemi et al. 25 , to map the material properties to tetrahedral elements 25 (Fig. 3e). Elastic moduli ranged from 1 MPa to ~25 GPa for elements of the proximal tibia. Unfortunately, images of the distal femur were not full-length. This lead to localized, unrealistic load transfer from the most proximal cross-section of femur (where the load was applied) to the tibial subchondral bone (where the FE-results were analyzed) when using a flexible distal femur. To address this limitation, we modeled the distal femur as a rigid body with E set to 500 GPa. Since this is a common method used with subject-specific FE modeling 17,42,43 , and our analysis is focused on the proximal tibia, we believe this decision is justified. All of the elements of bone tissue were modeled with isotropic linear material properties and a Poisson's ratio of 0.3 22 . For surrounding soft tissues, homogeneous, incompressible, and isotropic material properties were applied (E = 10 MPa, Poisson's ratio = 0.495) 22 .
Loading and Boundary Conditions. We fixed the proximal femur in all directions except the longitudinal axis of the femur, where we applied a uniform displacement of 1 mm. The most distal sections of the tibia and fibula were constrained for all degrees of freedom. To normalize the results based on the weight of each participant, the vertical reaction force at the top surface of femur was obtained for each FE model in the FE software. As linear elastic models were used, the stress and strain results were adjusted based on the ratio of the derived reaction force to the weight of each participant. This is similar to applying one body weight to the most distal section of the proximal femur to simulate single-leg stance position. FE Outcomes. FE-based stiffness as well as minimum principal (the most compressive) and von-Mises stress and strain distributions were acquired for the proximal tibia. Although some argue that von-Mises stress and strain should not be used for analyzing bone 44 , we included these measures for comparison with previous research 22,26 . Different regions including peripheral cortical, subchondral cortical, subchondral trabecular, epiphyseal cortical, epiphyseal trabecular, metaphyseal cortical, and metaphyseal trabecular were defined for regional analysis of FE-based mechanical metrics (MATLAB) (Fig. 1). The regions were defined using a custom MATLAB code. The depth of each region was defined as per previous research 20,45 .
1. cortical, 0-2.5 mm from the outer surface of the tibia, 2. subchondral cortical, 0-2.5 mm from the surface of tibia in the plateau and spine, 3. subchondral trabecular, 2.5-5 mm from the surface of tibia in the plateau and spine, 4. peripheral, 0-5 mm depth from the tibial surface along the outer cortical region, 5. epiphyseal, 5-15 mm from the surface of tibia in the plateau and spine, and 6. metaphyseal, 15-35 mm from the surface of tibia in the plateau and spine.
Stiffness was calculated for each medial and lateral compartment by developing two additional models for each knee. The distal femur was divided into two compartments in these models, and a low elastic modulus (E = 10 MPa) was assigned to the opposite compartment of interest (i.e., when calculating medial proximal tibial stiffness, the lateral distal femur compartment was assigned a low elastic modulus) (Fig. 3f). This was done to ensure that all load was transferred to the compartment of interest in the proximal tibia (i.e., it avoided load sharing between the two compartments). Although this approach means that both halves of the tibia contribute (to some extent) to stiffness, we believe this stiffness measure is a reflective representation for individuals suffering from varus or valgus alignment, where all load could be transferred through either the medial or lateral compartment, respectively. Stiffness was calculated as the applied vertical load divided by the average vertical displacement of the subchondral bone surface nodes in the respective compartments.
Statistical Analysis. We assessed short-term precision by calculating the CV% RMS for each FE outcome as well as BMD, BMC and bone volume from 3 repeated scans of 14 individuals 34 . We compared minimum principal stress, von-Mises stress, minimum principal strain, von-Mises strain, and stiffness values in different regions of the proximal tibia between OA and normal knees using unpaired t-tests for normally distributed data and non-parametric Mann-Whitney U-tests for data that were not normally distributed. A variable with skewness or kurtosis Z-score outside of ±1.96 limits was considered to have a non-parametric distribution. We reported p-values and confidence intervals (CI) for each statistical test and considered an alpha level <5% to be statistically significant (IBM SPSS Statistics, Version 24.0. Armonk, NY, USA). 95% CI were reported for normal distributions while Hodges-Lehmann estimator was used to calculate confidence intervals for non-parametric distributions 46 . We also calculated Cohen's d effect sizes for between-group mean differences in FE outcomes in relation to SD 47 . Cohen's d > 0.8 was considered to be a large effect size with clinical significance 48 .