Contribution of joint tissue properties to load-induced osteoarthritis.

Objective: Clinical evidence suggests that abnormal mechanical forces play a major role in the initiation and progression of osteoarthritis (OA). However, few studies have examined the mechanical environment that leads to disease. Thus, using a mouse tibial loading model, we quantified the cartilage contact stresses and examined the effects of altering tissue material properties on joint stresses during loading. Design: Using a discrete element model (DEA) in conjunction with joint kinematics data from a murine knee joint compression model, the magnitude and distribution of contact stresses in the tibial cartilage during joint loading were quantified at levels ranging from 0 to 9 N in 1 N increments. In addition, a simplified finite element (FEA) contact model was developed to simulate the knee joint, and parametric analyses were conducted to investigate the effects of altering bone and cartilage material properties on joint stresses during compressive loading.
Results: As loading increased, the peak contact pressures were sufficient to induce fibrillations on the cartilage surfaces. The computed areas of peak contact pressures correlated with experimentally defined areas of highest cartilage damage. Only alterations in cartilage properties and geometry caused large changes in cartilage contact pressures. However, changes in both bone and cartilage material properties resulted in significant changes in stresses induced in the bone during compressive loading. Conclusions: The level of mechanical stress induced by compressive tibial loading directly correlated with areas of biological change observed in the mouse knee joint. These results, taken together with the parametric analyses, are the first to demonstrate both experimentally and computationally that the tibial loading model is a useful preclinical platform with which to predict and study the effects of modulating bone and/or cartilage properties on attenuating OA progression. Given the direct correlation between computational modeling and experimental results, the effects of tissue-modifying treatments may be predicted prior to in vivo experimentation, allowing for novel therapeutics to be developed.

Introduction

Abnormal or excessive mechanical forces have been implicated in the initiation and progression of osteoarthritis (OA) (Felson, 2000; Arokoski et al., 1993; Lapveteläinen et al., 1995). A degenerative joint disease, OA is characterized by cartilage degradation, subchondral bone changes, and osteophyte formation, leading to radiographic joint narrowing and loss of joint function (Felson, 2000; Felson, 2006; Wieland et al., 2005; Hunter and Felson, 2006). Abnormal forces from joint instability due to injury (Felson, 2000; Lohmander et al., 2007), excessive physical activity (Kaila-Kangas et al., 2011; Cameron et al., 2011) and obesity (Felson et al., 1988; Messier et al., 2004) have been associated with higher risks of OA. On the other hand, moderate exercise and reduced loading may be beneficial treatment options for the attenuation of the disease (Messier et al., 2004; Felson, 1992; Hochberg et al., 2012). In addition, tissue changes such as subchondral bone sclerosis suggest that tissue properties also can play a critical role in the response of the joint to abnormal loads during OA progression (Radin et al., 1972; Radin and Rose, 1986; Radin et al., 1984; Burr and Schaffler, 1997; Burr, 2004; Li et al., 2013; Kawcak et al., 2010; Hayami et al., 2006; Botter et al., 2011).

Although preclinical models have been used to investigate the biological progression of OA due to induced joint instability (Glasson et al., 2007; Culley et al., 2015; Kamekura et al., 2005; Ma et al., 2007; Frank et al., 2012), few explored the relationship between mechanical forces, whole joint tissue material properties, and OA joint pathology initiation and progression in vivo (Burr and Schaffler, 1997; Brown et al., 1984; Poulet et al., 2013; Adebayo et al., 2016). Instead, most studies focused on isolating and exploring the functional relationship between cartilage and applied mechanical forces ex vivo (Ateshian et al., 1997; Lai et al., 1985; Mow et al., 1980; Wong et al., 2008). Elucidating the joint stresses associated with tissue degeneration in OA progression in vivo would provide valuable knowledge to more precisely differentiate the loads that are detrimental to the joint from the potentially beneficial loading regimens with therapeutic value. Furthermore, examining how changes to the joint tissue material properties affect joint stresses during disease would enable better understanding of how to moderate mechanical forces without worsening OA-related tissue damage in vivo.

The non-invasive tibial loading model provides a controlled mechanical setting to examine the relationship between mechanical forces and tissue degradation during disease progression (Adebayo et al., 2016; Ko et al., 2013; Christiansen et al., 2015; Poulet et al., 2011). This preclinical model offers a controlled platform with which to modulate the mechanical forces engendered at the knee joint, thus potentially allowing for exploration of loading regimes (Ko et al., 2016), treatments (Ziemian et al., 2021; Adebayo et al., 2017) and/or genetic models (Poulet et al., 2013) that would alleviate stresses on the cartilage surface and modulate load-induced OA development or progression. Previous studies using this model evaluated the biological and structural progression of the disease, which develops in response to controlled mechanical loading (Ko et al., 2013; Poulet et al., 2011) with consistent joint kinematics (Adebayo et al., 2016). Cyclic tibial compressive loading induces OA progression after a single loading bout (Ko et al., 2016; Ziemian et al., 2021) and after daily loading for 1, 2, and 6 weeks (Ko et al., 2013; Poulet et al., 2011; Christiansen et al., 2015). However, only one study to date investigated the contact stresses in the joint associated with load-induced OA progression in this model (Poulet et al., 2013), with contact stresses reported only at the peak compressive load of 12 N. Thus, contact behavior throughout each load cycle remains unknown. Furthermore, the effect of modifying tissue material properties on joint contact stresses under loading has not been explored in this model. Understanding how changes in tissue material properties affect joint mechanics before and during disease progression would allow for the development of tissue-specific therapies aimed to potentially alleviate load-induced OA development and progression.

In this study, we quantified the average and peak cartilage contact stresses and their spatial distribution during compressive tibial loading using experimental joint kinematics data and discrete element analysis (DEA). In addition, we evaluated how changes in bone and cartilage tissue material properties affected the peak cartilage stresses, using a simplified finite element (FE) contact model. We hypothesized that the spatial distribution and localization of the highest contact stresses calculated computationally would correlate to areas of greatest tissue damage observed experimentally. We further hypothesized that differences in bone and cartilage material properties would lead to changes in stresses induced at the cartilage surface, thus providing a potential explanation of the relationship between mechanical forces, tissue properties, and OA joint pathology progression.

Methods

Kinematic analysis

The left knee joints of 12 adult C57Bl/6 (B6, 32-wk old) cadaver male mice were labeled with 100 μm-diameter radiopaque bone fiducial markers on the medial proximal region of the tibia and anterior distal region of the femur close to the knee joint, with minimal disruption to the soft tissue, as previously described and validated (Adebayo et al., 2016). The joints then were subjected to compressive tibial loading levels ranging from 0 to 9 N (Fig. 1A). Roentgen stereophotogrammetric analysis (RSA) was performed to evaluate quasi-static joint kinematics at every 1 N increment during loading. Using a custom-made calibration cage and a dental x-ray, radiographs of the labeled joint were taken in 2 planes at each increment of loading. Each joint was subjected to 3 loading trials. The three-dimensional location of each bone marker was then calculated. The Eulerian method for rigid body kinematics was used to quantify the absolute translation and rotation of the tibia and the relative motion of the femur to the tibia at loading levels from 0 to 9 N.

Fig. 1

A) Schematic of the mouse tibial loading device, and B) the loading protocol for each stepped loading trial applied to the joint. C) Knee joint kinematics were analyzed using roentgen stereophotogrammetric analysis (RSA) (Adebayo et al., 2016). Arrows denote bead locations on the tibia and femur. Reference frame beads evident along perimeter of image. Scale bar = 5.0 mm. D) One sample was scanned by microCT, manually contoured, and E) bead locations from RSA and microCT were aligned to produce F) point clouds of each bone geometry with cartilage (red) inserted between the two surfaces to calculate the contact forces at each joint position by discrete element analysis.

Discrete model development

One labeled joint was scanned by microcomputed tomography (microCT, GE eXplore CT-120) at a 25 μm resolution, which was sufficient to image the surface geometry of the tibia, femur, and fiducial markers. This joint served as the model geometry for discrete element analysis. The scan was converted to a solid model (Materialise Mimics Research software, Plymouth, MI). The femur, tibia, and bone markers were manually contoured, resulting in three separate surface geometries (Fig. 1B). These surface geometries were meshed with a triangular surface and exported as point clouds for DEA. DEA has been used extensively to study articular contact mechanics in which the articulating bones are considered rigid, and the cartilage tissues are regarded as a set of individual compression springs (Guess et al., 2011; Volokh et al., 2007; Anderson et al., 2010; Bei and Fregly, 2004). This numerical technique has been validated by experimental measurements and finite element models and shown to be within 4 % and 5 % of physical contact stress measurements and finite element analyses, respectively (Anderson et al., 2010; Li et al., 1997; Anderson et al., 2007). DEA provides a simple yet accurate framework to calculate compressive cartilage contact stresses, especially in situations for which experimental measurement techniques are difficult, as is the case here due to the small animal size.

Using custom Matlab code (2016b, The Mathworks, Natick, MA), the tibia and femur were registered and aligned to RSA joint locations at 0 N compressive load, as determined by the common location of each bone marker (Fig. 1C, D). With kinematics data from RSA for all samples, DEA was conducted at every 1 N increment of compressive loading, and point clouds of the tibia and femur were transformed and oriented to correct locations at loading levels from 0 to 9 N, accordingly. Because the cartilage was not visualized in the microCT scan, the tissue was assumed to be a uniform thickness of 100 μm on both the tibial and femoral surfaces. Using previously established DEA protocols (Volokh et al., 2007; Anderson et al., 2010; Gardner-Morse et al., 2013), an elastic spring contact model was used to determine contact stresses in the tibial cartilage. Spring stiffness (k) and contact stress (σ) were calculated as follows:where E and ν are the Young's modulus and Poisson's ratio of the cartilage, respectively; h is the assumed total cartilage thickness (200 μm); and δ is the deformation in the tibial cartilage, as determined by the distance between the tibial and femoral surfaces and resolved along the normal vectors on the tibial surface. Based on previous literature (Wong and Carter, 1990; Beaupré et al., 2000), the Young's modulus and Poisson's ratio of cartilage were assumed to be 6 MPa and 0.47, respectively. The average forces acting on the tibial surface also were calculated by multiplying the total contact area by the average compressive stress calculated for each tibial plateau.

Sample measurements for finite element model

A simplified FE model was developed to validate the results from the discrete element model and to conduct a parametric analysis to evaluate the effects of changing joint tissue properties on cartilage stresses during loading. Right knee joints from adult C57Bl/6 male mice (n = 6) were scanned using microCT at 10 μm resolution (μCT35, Scanco: 55 kVp, 145 μA, 600 ms integration time). Each scan was converted to a solid geometry (Mimics, Materialise). We measured the approximate radius of curvature of the medial tibial plateau and femoral condyle by fitting circles to each surface. We also measured the radius of concavity on the tibial medial plateau and the depth of the concavity from the intercondylar eminence (Fig. 2A). We specifically selected the medial section of the joint as we have previously demonstrated that cartilage degradation and osteophyte formation localize to this area during loading (Ko et al., 2013). Cartilage and subchondral bone thickness measurements were obtained from Safranin O-stained histological slides and approximated for the purposes of the FE model.

Fig. 2

A) Tibial (yellow) and femoral (red) radii of curvature were measured on the medial aspect of the joint by fitting of spheres. The radius of curvature and depth of the tibial concavity (yellow arrow) was also measured. C) Simplified geometric contact model with the noted boundary conditions for finite element analysis based on the geometry measurements made of the contacting surfaces. B) Mesh convergence analysis concluded that approximately 15,000 cartilage elements (red arrow) were required for accurate contact pressure results.

Simple contact finite element model

Using the average geometric measurements from the samples, semi-spherical shell geometries were created with the appropriate radii for the tibial and femoral surfaces (Fig. 2C). The tibia also included a concavity with the appropriate depth as measured from the microCT scans (Fig. 2A). The outermost layer for each surface was cartilage, which was modeled as a linear elastic tissue with a uniform thickness of 100 μm, Young's modulus of 6 MPa, and Poisson's ratio of 0.47, as for the DEA (Table 1). Adjacent to the cartilage, the subchondral cortical bone layer was assigned a uniform thickness of 130 μm, Young's modulus of 18GPa, and Poisson's ratio of 0.3, based on values from the literature (Rho et al., 1993). The interior semi-spherical shell geometry of the epiphysis was modeled as cancellous bone with a Young's modulus of 568 MPa and Poisson's ratio of 0.3, calculated using the average epiphyseal bone volume fraction in each tibia (Vijayakumar and Quenneville, 2016). Tied constraints were modeled between both the epiphyseal and subchondral cortical layers and the cortical and cartilage layers, and frictionless contact was modeled between the tibial and femoral cartilage surfaces. To decrease computational expense and take advantage of the symmetry of the geometry, only a quarter model was created, and symmetric boundary conditions were applied to appropriate surfaces (Fig. 2C).

Table 1

Geometric and material property values for the parametric analysis conducted on the simple contact finite element model.

	Cartilage			Subchondral cortical plate			Epiphyseal cancellous bone
	Thickness (μm)	Young's modulus (E, MPa)	Poisson's ratio (ν)	Thickness (μm)	Young's modulus (E, MPa)	Poisson's ratio (ν)	Young's modulus (E, MPa)	Poisson's ratio (ν)
Normal geometry	100	6	0.47	130	18,000	0.3	568	0.3
0.5 cartilage thickness	50	6	0.47	130	18,000	0.3	568	0.3
2× cartilage thickness	200	6	0.47	130	18,000	0.3	568	0.3
0.5 cartilage E	100	3	0.47	130	18,000	0.3	568	0.3
2× cartilage E	100	12	0.47	130	18,000	0.3	568	0.3
0.5 subchondral plate thickness	100	6	0.47	65	18,000	0.3	568	0.3
2× subchondral plate thickness	100	6	0.47	260	18,000	0.3	568	0.3
0.5 subchondral plate E	100	6	0.47	130	9000	0.3	568	0.3
2× subchondral plate E	100	6	0.47	130	36,000	0.3	568	0.3
0.5 epiphyseal cancellous E	100	6	0.47	130	18,000	0.3	284	0.3
2× epiphyseal cancellous E	100	6	0.47	130	18,000	0.3	1136	0.3

All tissues were meshed with linear 8-node brick elements. The element size for each cartilage surface was refined during a mesh convergence analysis. A 0.5 N total force was applied to the proximal surface of the cancellous and cortical bone, and the distal surface of the tibia remained fixed in all directions. Mesh convergence was conducted to determine the minimum element density on each cartilage surface that would result in the most accurate measures of contact pressures. An element density was chosen such that any further increase in element density would result in less than 1 % change in resulting contact pressure and maximum principal stress. Thus, the cartilage on the tibial and femoral surfaces was meshed with approximately 8100 and 6825 elements, respectively (Fig. 2B).

Validation of discrete element analysis by finite element model

To validate the compressive stresses determined using DEA, the average contact forces calculated on the medial tibial plateau at loading levels from 0 to 9 N were applied to the simplified contact geometry. At each load, a quarter of the average contact force on the medial plateau was applied to account for the quarter geometry in the FE model. Peak contact stresses then were compared in the two models at loading levels from 0 to 9 N.

Parametric analysis

A parametric analysis was conducted on the simple contact FE model to examine the effects of changes in tissue properties and mass on the stresses induced at the cartilage surface. Variations in cartilage modulus and thickness, subchondral bone modulus and thickness, and epiphyseal cancellous modulus were examined (Fig. 3, Table 1). Specifically, in addition to the normal geometry, altered geometries included a 50 % decrease and 100 % increase in cartilage and bone thicknesses and moduli. Changes were made on both the tibial and femoral sides of the model. Changes in contact pressure and maximum principal stress and strain at the cartilage surface and in the cortical and cancellous bone were compared.

Fig. 3

Geometric property values for the parametric analysis conducted on the simple contact finite element model. Cartilage layer in purple; subchondral cortical plate in green; epiphyseal geometry in yellow.

Results

Compressive loading resulted in two distinct contact behaviors

The results of the kinematics analysis of intact vs. transected joints under loading were published previously (Adebayo et al., 2016). Briefly, for intact joints, the tibia translated primarily anteriorly and proximally with little rotation around all axes as loading increased. In terms of relative femur to tibia motion, knee flexion increased as load magnitudes increased from 0 to 9 N. Using the tibial and femoral locations for loads ranging from 0 to 9 N, corresponding contact stresses were calculated using DEA. Due to the assumption of 100 μm uniform cartilage thicknesses on the tibial and femoral surfaces, 0 N load engendered an initial average compressive stress of 1.00 MPa with a peak compressive stress of 4.60 MPa on the tibial cartilage surface.

As tibial loading increased from 0 to 9 N, two distinct behaviors were observed. Regardless of individual samples, 19 of 36 trials exhibited primarily compression of the cartilage between the tibial and femoral surfaces as knee flexion increased with load magnitude. In contrast, 17 trials exhibited “rolling” of the femoral surface away from the tibia as loading increased, and cartilage compression did not increase with load (Fig. 4). To ensure that the two loading behaviors were not a function of creep errors due to quasi-static kinematic measurements at every 1 N increment of load, we also evaluated 5 samples that were loaded directly from 0 to 9 N. Four of five samples exhibited compressive behavior with an average peak stress of 7.70 ± 1.83 MPa at 9 N, while one exhibited rolling behavior with a peak stress of 4.19 MPa at 9 N, thus confirming the two contact behaviors during tibial compressive loading. Henceforth, in this manuscript, we will focus on the trials that exhibited primarily compressive behavior. Results of the trials with rolling behavior can be found in the supplemental data (Fig. S-1).

Fig. 4

Side and top views of a trial exhibiting compressive behavior and a trial exhibiting rolling behavior as seen by comparing the 0 N and 9 N positions. The compressive behavior maintained contact between the tibia and femur.

For trials that exhibited primarily compressive behavior, the average compressive stress almost tripled from 1.00 MPa at 0 N to an average of 2.82 MPa at 9 N, with peak compressive stress doubling from 4.60 MPa to 8.98 MPa as the load increased (Fig. 5A, B). Strain values averaged 2.78 % (12.75 % peak strain) at 0 N and increased by 181 % to 7.82 % (24.91 % peak strain) at 9 N (Fig. 5C, D). The average forces acting on each tibial plateau increased similarly to 2.19 N and 2.71 N at 9 N load for the medial and lateral plateaus, respectively (Fig. 5E). Peak contact stresses/strains occurred in the middle of each tibial plateau at 0 N, and moved posteriorly, increasing in magnitude as loading increased (Fig. 5F).

Fig. 5

In trials exhibiting compressive behavior, A) mean and B) peak contact stresses increased with load magnitude (mean ± SD). C) Mean and D) peak contact strains also increased with load magnitude. E) Contact forces increased in magnitude with load, and F) peak contact stresses translated posteriorly on the tibial plateau as loading increased.

Simple finite element contact model validated DEA peak compressive stresses

To validate the contact results from DEA, peak loads calculated specifically from trials exhibiting compressive behavior were applied to the simple FE model. Peak compressive stresses were similar between the FE and DEA models (Fig. 6A). The greatest difference in calculated peak compressive stresses occurred at 0 N with a 21 % difference between DEA and FEA results. These differences reduced as the load increased, with a minimum difference in compressive stress of 3.6 % at 8 N between the two models. The small differences in calculated stresses between the two methods validate the contact results from the discrete element analysis.

Fig. 6

A) Peak contact pressures calculated using finite element analysis validated contact stress values determined by discrete element analysis (mean ± SD). B) A 0.5 N compressive load resulted in a peak contact stress of 6.37 MPa in the middle of the tibial cartilage surface.

Only cartilage changes affected cartilage contact pressure in the simplified finite element model

Parametric analyses of the effects of cartilage and bone material and geometric properties on contact stresses revealed that changes to the cartilage had significant effects on contact mechanics at the joint (Table 2). Normal joint geometry and properties produced a peak contact pressure of 6.37 MPa in the middle of the tibial cartilage surface (Fig. 6B). A 50 % decrease in cartilage thickness increased contact pressure by 29 % (8.20 MPa), whereas doubled thickness decreased contact pressure by 23 % (4.92 MPa). Furthermore, a 50 % decrease in the Young's modulus of the cartilage accounted for a 12 % decrease in contact pressure, and doubling of Young's modulus resulted in a 14 % increase in contact pressure. Changes in the Young's modulus and the cortical thickness of the subchondral cortical and cancellous bone did not significantly alter contact pressures at the cartilage surface. Differences in contact pressures associated with bone changes ranged from a 1.05 % decrease to a 0.42 % increase depending on bone material or thickness change.

Table 2

Peak contact pressure and max principal stresses and strains induced on the cartilage, subchondral plate and cancellous bone due to changes in tissue geometry or Young's modulus. Negative values indicate compression.

	Cartilage				Subchondral cortical plate		Epiphyseal cancellous bone
	Peak contact pressure (MPa)	Change in contact pressure (%)	Max principal stress (MPa)	Max principal strain (%)	Max principal stress (MPa)	Max principal strain (%)	Max principal stress (MPa)	Max principal strain (%)
Normal geometry	6.37	–	−6.36	−0.39	−20.24	−8.59E-04	−2.23	−3.61E-03
0.5 cartilage thickness	8.20	+29	−8.20	−0.34	−25.14	−9.59E-04	−2.73	−4.48E-03
2× cartilage thickness	4.92	−23	−4.56	−0.49	−9.09	−4.93E-04	−1.54	−2.51E-03
0.5 cartilage E	5.61	−12	−5.59	−0.68	−18.84	−8.04E-04	−2.06	−3.32E-03
2× cartilage E	7.24	+14	−7.22	−0.23	−21.12	−8.73E-04	−2.35	−3.84E-03
0.5 subchondral plate thickness	6.31	−0.9	−6.29	−0.39	−20.46	−1.04E-03	−4.30	−6.59E-03
2× subchondral plate thickness	6.40	+0.5	−6.38	−0.39	−9.28	−3.35E-04	−0.64	−1.09E-03
0.5 subchondral plate E	6.35	−0.3	−6.33	−0.39	−12.76	−1.17E-03	−2.94	−4.76E-03
2× subchondral plate E	6.39	+0.3	−6.37	−0.39	−28.10	−5.31E-04	−1.49	−2.41E-03
0.5 epiphyseal cancellous E	6.36	−0.2	−6.34	−0.39	−27.41	−1.04E-03	−1.48	−4.79E-03
2× epiphyseal cancellous E	6.39	−0.3	−6.37	−0.39	−12.88	−5.95E-04	−2.96	−2.40E-03

When bone material properties were varied, the effects on cortical and cancellous bone stresses differed (Table 2). Decreased cancellous Young's modulus led to increased peak stresses in cortical bone and decreased stresses in cancellous bone. Conversely, decreased cortical Young's modulus decreased stresses in the cortical bone and increased stresses in the cancellous bone.

Because bone is responsive to dynamic strain values (Radin et al., 1984), strain values are important to consider. Trends in bone strains generally followed the same patterns as the bone stresses, except for changes in bone material properties. Decreased cancellous and cortical Young's moduli increased cancellous and cortical peak strains, respectively (Table 2).

Discussion

This study is the first to examine contact stresses in the joint under compressive loading using both discrete element and finite element analyses. We demonstrated that compressive loading of the mouse tibia resulted in two contact behaviors – a primarily compressive behavior and a femoral “rolling” behavior. In trials that exhibited compressive contact behavior, maximum contact pressure increased from 4.60 MPa at 0 N to approximately 9 MPa at 9 N. The failure or flexural strength of cartilage previously calculated in cyclic compression studies ranged between 15 and 50 MPa, and depended on loading frequencies (Sadeghi et al., 2017; Kerin et al., 1998). Although these stress values are higher than the 9 MPa calculated in our study, cartilage surface fibrillations can occur with as little as 8 MPa under static compression (Fick and Espino, 2011; Fick and Espino, 2012), or from 3 to 7 MPa over 10,000 compressive cycles applied with an indenter (Sadeghi et al., 2015). Thus, the application of 9 MPa over multiple bouts of 1200 compressive cycles at 4 Hz is likely sufficient to induce cartilage surface damage in vivo. These values also may explain the absence of macroscopic cartilage damage after a single tibial loading bout. A single bout of 1200 cycles may not induce mechanical damage at the cartilage surface but instead may induce cellular or signaling changes that lead to downstream degradation of cartilage and OA development. The applied contact pressure at 9 N over multiple loading bouts likely leads to cartilage damage at the surface that may progress through the thickness with the application of each additional loading bout.

Previous clinical (Fukubayashi and Kurosawa, 1980; Segal et al., 2009; Segal et al., 2012) and preclinical (Gardner-Morse et al., 2013; Brown et al., 1991) studies measured cartilage contact pressures in the healthy knee during ambulation from 0.5 to 4 MPa (Brand, 2005; Morrison, 1970). Thus, the pressures predicted in this study, particularly with compressive behavior, exceeded cartilage pressures measured in normal healthy tissue, also possibly explaining the cartilage damage observed in vivo (Segal et al., 2009; Segal et al., 2012). These pressure values also are comparable to those reported in a FE model with similar compressive loading (Poulet et al., 2013). Using an average adult male mouse body mass of 30 g, the average contact forces calculated with DEA under 9 N compression corresponded to approximately 17 times body weight. Clinically, these forces are more than those engendered during full knee extension landing (Makinejad et al., 2013), vertical jumping (Cleather et al., 2013), kneeling and squatting (Nagura et al., 2006), activities that have been suggested to lead to increased cartilage contact pressures, injury, and potentially subsequent cartilage degradation (Felson, 2000; Lohmander et al., 2007; Makinejad et al., 2013; Felson et al., 1991; Coggon et al., 2000).

Although the reason for the two distinct contact behaviors is not understood, we hypothesize that joint and mouse positioning in the compressive loading device may be a potential explanation. Relative femoral-to-tibial proximal and posterior translations were significantly different between the two contact behaviors at 9 N (Fig. S-2). Because knee flexion is not restricted in this model, the relative location of the femur to the tibia at 0 N could cause two distinct contact behaviors. Further in vivo studies are needed to examine which behavior occurs during multiple loading cycles. The damage evident in the knee articular cartilage post-loading suggests that compression behavior is dominant.

Despite differences in contact behaviors, the peak compressive stress location in all trials translated to the posterior aspect of the joint as loading increased from 0 to 9 N. The location of peak compressive stress at 9 N correlated with the area in which the most severe cartilage damage occurs in this model. Using histological scoring on Safranin O-stained slides, several studies demonstrated that loading induces the most severe cartilage damage on the posterior aspect of the tibial plateau (Ko et al., 2013; Ko et al., 2016). As loading increased from 0 to 9 N, peak compressive stresses translated to the posterior aspect of the joint, thus confirming the correlation between highest contact stresses and the most severe cartilage damage. Furthermore, contact pressures were not evenly distributed on medial and lateral surfaces, indicating varus or valgus rotation as loading increased. These rotations can produce tension in ligaments at the joint margins and could potentially explain bone formation at the entheses of these ligaments (Benjamin et al., 2006; Rogers et al., 1997). Further studies are needed to examine the loads on the ligaments, confirm excessive tension on the joint margins and differentiate between the formation of osteophytes and enthesiophytes in this model.

To validate the results of the DEA model, the simplified FE contact model confirmed the contact pressures calculated during loading. The normal geometry demonstrated comparable contact pressures and stresses on the cartilage surfaces. Parametric evaluation of the effects of joint tissue changes on stresses in cartilage surfaces revealed that only cartilage changes significantly affected contact pressures. Similar to previous results (Poulet et al., 2011; Li et al., 2001; Huber-Betzer et al., 1990), increased cartilage thickness and reduced cartilage modulus decreased cartilage surface contact pressures. Similar to results from other computation studies (Burr and Schaffler, 1997; Brown et al., 1984), bone changes did not significantly affect cartilage contact pressures; however, they affected the peak stresses and strains engendered in the cortical and cancellous bone.

Bone adapts to dynamic changes in applied mechanical strains (Meakin et al., 2014; Rubin and Lanyon, 1985; Lynch et al., 2011). Based on our parametric analysis, the differences in peak bone stresses and strains would result in bone adaptation in vivo. Therefore, although changes in bone thickness and properties may not directly affect cartilage surface stresses during loading, changes in bone adaptation and remodeling could result and, thus promote changes in bone/cartilage crosstalk (Burr, 1998; Findlay and Kuliwaba, 2016; Yuan et al., 2014). Bone remodeling has been implicated in exacerbating bone/cartilage crosstalk during OA disease progression, in which biomolecular signals can be easily transferred from one tissue to the other and bone tissue and cells advance through vascular invasion into the overlying cartilage (Findlay and Kuliwaba, 2016; Yuan et al., 2014). Thus, bone tissue changes may drive bone adaptation, resulting in subchondral bone sclerosis and subsequent tissue crosstalk between the cartilage and subchondral bone during abnormal mechanical loading.

While the DEA and FE contact models allowed us to examine the contact mechanics of the joint during loading, the use of these models to understand complex biological states has some limitations. Both models only considered compressive stresses. Knee joints are subjected to shear and tensile forces and exhibit more complex stress states in vivo. Thus, future studies should examine the shear and tensile stress states under tibial loading. Furthermore, our discrete element analysis assumed that the bone tissue was rigid. While our FE model eliminated this assumption, the model is linearly elastic, which may not accurately reflect tissue behavior, particularly for cartilage. Cartilage has been modeled extensively as a linear elastic, nearly incompressible material, for situations capturing relatively short or instantaneous time-frames (Volokh et al., 2007; Ateshian et al., 2015). However, biologically, cartilage tissue is biphasic (Ateshian et al., 1997; Lai et al., 1985; Mow et al., 1980). Further studies are needed to account for other joint motions and the cartilage fluid phase during longer loading durations.

In these models the menisci and other soft tissues in the joint were not considered, yet play major roles in the mechanical stability of the joint. For instance, because the DEA did not account for the thickness of the menisci, we may potentially overestimate the contact stresses in areas of articular cartilage covered by menisci. In addition, we did not accommodate for a calcified hyaline cartilage zone, which would effectively reduce the thickness of the functional cartilage zones and/or increase the overall stiffness of the cartilage. Based on our parametric analysis, such changes in cartilage properties would significantly increase peak contact pressures on the cartilage surface. The use of a single model geometry and uniform cartilage thickness is another limitation of this study. However, the FEA parametric analysis provides insight into how variations in the articulating radii of curvatures and/or articular cartilage thickness would affect cartilage contact pressures. Specifically, changes in the bone geometry had minimal effect on cartilage contact pressures, whereas changes in the cartilage thickness would have a significant effect. Nonetheless, considering the size of the mouse joint and the limited imaging capabilities, this study provides an initial quantification of contact stress magnitudes and locations during in vivo loading of the mouse knee. More complex, customized geometric FE models are needed to quantify the contributions other soft tissues to the joint mechanical environment.

In conclusion, peak contact stress magnitude and location correlated with areas of severe cartilage damage during tibial compressive loading. Unlike changes in cartilage thickness and properties, the changes in bone properties and mass did not directly affect contact stresses on the cartilage surface. However, changes in all tissues did affect stresses and strains induced in the cortical and cancellous bone. Thus, while material-related changes to the bone may not directly influence stresses on the cartilage surface, these changes may lead to differences in bone adaptation and potentially promote subsequent tissue crosstalk between bone and cartilage during loading. Further in vivo studies are needed to investigate the molecular and cellular crosstalk associated with increased strains engendered on the bone. These studies would help to determine the compressive loads and contact pressures that could potentially promote cartilage health and attenuate OA progression.

CRediT authorship contribution statement

Conception and design: OOA, SRG, MBG, TMW, MCHM.

Acquisition, analysis, and interpretation of the data: OOA, SRG, MBG, TMW, MCHM.

Drafting and critical revision of the article for important intellectual content: OOA, SRG, MBG, TMW, MCHM.

Final approval of the article: OOA, SRG, MBG, TMW, MCHM.

Declaration of competing interest

The authors have no conflict of interest related to this work.