Biomechanics of the medial meniscus in the osteoarthritic knee joint

Geometry

A model of the left knee joint was prepared on the basis of magnetic resonance imaging (MRI) of a healthy volunteer (female, 43-year-old, 1.75 m, 68 kg) in a full extension position. A 1.5 T MRI scanner (Magnetom Aera, Siemens, resolution of 0.28125 mm) was used to obtain the images of the knee in the axial plane (slice spacing = 0.7 mm) and the sagittal and coronal planes (slice spacing = 3.6 mm). Written informed consent was obtained from patient to use data from her medical records in research. The study was approved by Independent Bioethics Committee for Research at Medical University of Gdańsk, Poland. Then, the MRI scans were imported into Mimics (Materialise NV, Leuven, Belgium), where semi-automatic segmentation was performed by an experienced user and verified by a radiologist. The obtained three dimensional (3D) model of the knee joint geometry consisted of bones (femur, tibia, fibula), menisci and articular cartilage (femoral and tibial).

Finite element models

An intact finite element (FE) model of the knee joint was developed in Abaqus 6.14-2 (Dassault Systemes Simulia Corp., Providence, RI, USA) based on the 3D geometry imported from Mimics. Ligaments were modelled using two node, 3D axial connector elements CONN3D2. The following ligaments were included in the FE model: anterior cruciate ligament (ACL), anterior intermeniscal ligament (AIML), anterolateral ligament (ALL), lateral collateral ligament (LCL), medial collateral ligament (MCL), posterior cruciate ligament (PCL) and posterior oblique ligament (POL). The attachment points of the specific ligament bundles were identified in the MRI images. The medial collateral ligament was divided into deep bundle (dMCL) and three superficial bundles: anterior, middle and posterior. Two portions of the dMCL were distinguished: meniscotibial ligament (MTL) connecting the medial meniscus to the tibia and meniscofemoral ligament (MFL) connecting the medial meniscus to the femur. The MTL, MFL and AIML ligaments were modelled by the six axial connector elements to minimize stress concentrations at the attachment points to the medial meniscus.

Three additional models were prepared to examine the influence of medial knee OA on the biomechanics of the medial meniscus. Model G was introduced to assess the effect of only altered geometry of the articular cartilage corresponding to KL score 3 and model M to evaluate the influence of only decrease in material parameters characteristic of medial knee OA. Degenerative changes in the OA model were simulated by both non-uniform reduction in cartilage thickness and reduction in material parameters of cartilage and menisci.

The crucial factor in medial compartment knee osteoarthritis is the degradation of articulate cartilage (Wirth et al., 2010; Arno et al., 2012). In the present study alterations of medial knee OA geometry were introduced into model G and the OA model based on difference thickness maps from an MRI study (Favre et al., 2017) comparing the cartilage thickness of the KL 3 subgroup to the non-OA subgroup. The cartilage thickness was reduced in specific areas of the medial femoral condyle and both compartments of the tibial plateau, where significant differences were observed in the study (Favre et al., 2017). The cross-section of the menisci and articular cartilage in the coronal plane for models with intact geometry and models with altered geometry is presented in Figs. 1C and 1D, respectively. The continuous non-uniform offset of contact surfaces of the articular cartilage resulted in a decrease of 9.2 and 7% in the volume of tibial cartilage in the medial compartment and the lateral compartment, respectively.

Materials

The Yeoh model was used to model an isotropic, hyperelastic, nearly incompressible behavior of the articular cartilage. The strain energy function for the Yeoh model has the following form (1)${\displaystyle \Phi =C_{10}\left({\overline{I}}_1-3\right)+C_{20}{\left({\overline{I}}_1-3\right)}^2+C_{30}{\left({\overline{I}}_1-3\right)}^3+\frac{1}{D_1}{\left(J-1\right)}^2+\frac{1}{D_2}{\left(J-1\right)}^4+\frac{1}{D_3}{\left(J-1\right)}^6,}$

where ${\overline{I}}_1$ is the first invariant of the modified left Cauchy–Green tensor, J is the elastic volume ratio and C₁₀, C₂₀, C₃₀, D₁, D₂, D₃ are independent Yeoh material constants. The values of the parameters C₁₀, C₂₀ for normal and degenerated knee cartilage were taken from an experimental study (Robinson et al., 2016), see Table 1. In the paper, the constant C₃₀ was assumed to be zero due to a lack of experimental data. The value of the parameter D₁ was calculated based on the value of C₁₀ and the value v = 0.45 of Poisson’s ratio (Table 1). Moreover, the last two terms in the Eq. (1) were neglected because a nonzero value of D₁ was sufficient to model a nearly incompressible behavior of the articular cartilage.

A linear, elastic and transversely isotropic material model was used for the menisci. The Young’s moduli in the circumferential direction E_θ and in the cross sectional plane E_p, Poisson’s ratios ν_θp, ν_p and the shear modulus G for healthy menisci were taken from papers (Lechner, Hull & Howell, 2000; Haut Donahue et al., 2003; Yang, Canavan & Nayeb-Hashemi, 2010), see Table 2. The material parameters of osteoarthritic menisci in models M and OA were assumed based on experimental data for the menisci with a grade 3 of degeneration (Fischenich et al., 2015). A significant decrease in the instantaneous compressive modulus was observed in both osteoarthritic menisci, therefore four times smaller Young’s modulus E_p was assumed in the axial and radial directions in which menisci are compressed (Fischenich et al., 2015; Fischenich et al., 2017), see Table 2. The unchanged Young’s modulus was used in the circumferential direction, because degenerated menisci retain a tensile modulus similar to that obtained for healthy menisci (Lechner, Hull & Howell, 2000; Fischenich et al., 2015). The same values of the Poisson’s ratios and the shear modulus were assumed in all models due to a lack of experimental study about effect of medial knee OA on these parameters.

A constitutive relation in the ligaments modelled by the axial connector elements was defined by relationship between the ligament tensile force F and its elongation u. The force in the ligaments is zero when the following condition is fulfilled: (2)${\displaystyle F=0,u<-{\varepsilon }_r{L}_0,}$

where ɛ_r is the reference (initial) strain and L ₀ is the zero-length of a ligament bundle. The length L ₀ was computed based on the ligament length L_r and the strain ɛ_r in the reference configuration (3)${\displaystyle L_0=L_r/\left(1+{\varepsilon }_r\right).}$

For low values of elongation the relationship F(u) is nonlinear and given by (4)${\displaystyle F=\frac{k{\left(u+{\varepsilon }_r{L}_0\right)}^2}{4L_0^2{\varepsilon }_l},-{\varepsilon }_r{L}_0\le u\le \left(2{\varepsilon }_l-{\varepsilon }_r\right)L_0,}$ where k is the ligament stiffness and ɛ_l is the transition strain, assumed as 0.03 (Butler, Kay & Stouffer, 1986). The values of k and ɛ_r are given for the whole ligaments and the specific ligament bundles in Table 3. For greater strains the ligament force is proportional to its elongation (5)${\displaystyle F=k\left(\frac{u}{L_0}+{\varepsilon }_r-{\varepsilon }_l\right),u>\left(2{\varepsilon }_l-{\varepsilon }_r\right)L_0.}$

The meniscal horn attachments were modelled by nonlinear spring elements carrying only tensile forces. Insertion nodes for each meniscal attachment were identified on the basis of the MRI scans. The stiffness of each meniscal horn attachment was assumed based on the experimental linear stiffness, considering difference between its length in the FE models and an experimental study (Hauch, Villegas & Haut Donahue, 2010). The total tensile stiffness k_h and the spring stiffness k_s are reported for the meniscal horn attachments in Table 4.

The use of reduced material parameters for degenerated cartilage and menisci is the authors’ contribution to the description of material parameters for FEM models. The density for all soft tissues was assumed as 1,500 kg/m³ (Erdemir & Sibole, 2010) and for the bones as 2,000 kg/m³ (Maharaj, Maheswaran & Vasanthanathan, 2013).

Finite element mesh

The menisci and the articular cartilage were discretized by the 4-node tetrahedral elements (C3D4) with linear shape functions. The rigid shell triangular (R3D3) and quadrilateral (R3D4) elements were used to model the bones. The average size of 0.8 mm, 1.1 mm and 1.2 mm was assumed for the menisci, the tibial cartilage and the femoral cartilage, respectively. The intact knee model consisted of the total number of 69,808 nodes and 267,479 finite elements. The mesh convergence test was performed to verify the correctness of FE mesh. A refined discretization (319,274 nodes) was generated by reduction of finite element size by a factor of two in the menisci and articular cartilage. The differences in contact pressures and hoop stresses between the refined discretization and the one used in the paper were less than 4%. Thus, the assumed FE mesh was considered adequate for numerical simulations (Zienkiewicz & Taylor, 2002).

Finite element analysis

Only stance phase of the gait cycle was simulated in the FE code Abaqus, because it has been well documented in-vivo (Gilbert et al., 2014; Bergmann et al., 2014) and through modelling (Shriram et al., 2019) that tibiofemoral contact forces and stresses are small during the swing phase compared to the stance phase. In the initial static step the equilibrium state between contacting surfaces was obtained for the imposed reference strains in the ligaments (Table 3). Then the initial loads from the gait cycle were applied to the model within one second in a dynamic step. The dynamic implicit solver was used in the nonlinear FE analysis of the stance phase that accounts for 60% of the duration of the gait cycle of 1.0674 s (Bergmann et al., 2014).

Boundary conditions and loads

The femur, tibia and fibula were modelled as rigid bodies due to significantly greater stiffness of bones than soft tissues. Kinematics of the bones is defined by the six degrees of freedom of their so-called reference points (RPs). The tibia RP was located on the intersection of the tibial mechanical axis and the surface of the tibial plateau. The femur RP was assumed at the midpoint between the medial and lateral epicondyles. The orthogonal joint coordinate system was defined based on the tibial mechanical axis and the femur flexion axis. The positive directions for translational and rotational degrees of freedom in the knee joint are presented in Fig. 1B.

Forces and moments acting in the knee during the gait cycle were taken from experimental studies (Bergmann, 2008; Bergmann et al., 2014), because gait data was not available for the analyzed knee. First, the standardized average loads (file Walking_AllLoadLevels from database OrthoLoad (Bergmann, 2008)) were multiplied by 68/75 to re-calculate them from an average body weight of 75 kg to the volunteer weight of 68 kg. Then the forces and moments were transformed from the implant-based system to the knee joint coordinate system. The obtained loads (Fig. 2) were applied to the tibia RP and they followed rotations of the tibia RP during the gait simulation.

The greater reduction in thickness of articular cartilage in the medial compartment resulted in approximately 1° greater varus rotation in models G and OA than in the intact model in gait simulations. It is well established that the greater varus alignment causes an increase in knee varus moment in patients with knee OA (Moyer et al., 2010; Sharma, 2007). However, our model of the knee joint cannot simulate an increase in varus moment as the effect of medial movement of the mechanical axis of the leg with respect to the knee. Therefore, corrected varus moment (denoted by dashed line in Fig. 2) computed for 1° greater varus alignment was used in models G and OA to directly take into account this effect. The loading cycles from 7 subjects (files K1L, K2L, K3R, K5R, K7L, K8L, K9L from the database OrthoLoad (Bergmann, 2008)) were averaged intra-individually using basic average method (Bergmann et al., 2014) to determine corrected standardized average varus moment for the volunteer weight of 68 kg. Load data for K6L subject with valgus alignment (from the database OrthoLoad (Bergmann, 2008)) was not used in these computations to obtain the increased varus moment for 1° greater varus alignment than in the intact model. The transformation of the corrected varus moment to the knee joint coordinate system resulted in a slight change in the internal moment, see Fig. 2.

All degrees of freedom (DOFs) were fixed at the femur RP. During the gait simulation five degrees of freedom were released at the tibia RP, except flexion rotation which was prescribed (see Fig. 2D). The hard contact approach was used to model interaction between external surfaces of the articular cartilage and menisci. A friction with the very low friction coefficient of 0.02 was assumed between the contacting surfaces based on a previous study (Warnecke et al., 2019). Moreover, an automatic stabilization factor of 0.2 was used to improve the convergence of computations in the initial step. A “tie” constraint was defined between the bone RPs and cartilage surfaces to enforce the rigid connection at the cartilage-bone interfaces. The fibula RP was constrained with the tibia RP using a “kinematic coupling” option.

Model validation

The FE models were validated against the expected results from literature. First, the passive knee flexion test up to 70 degrees was performed for the intact model and the OA model to verify correctness of the tibiofemoral kinematics. The tibial kinematics with respect to the femur predicted by the FE models were compared with the experimental tibiofemoral kinematics measured in 23 cadaveric lower limbs (Germain et al., 2016). Maximal tibial contact pressures obtained in the intact model for the first and second peaks of total force were verified against the experimental values measured with pressure sensors in a study (Bedi et al., 2010). Moreover, the change in location of the peak contact pressure between the first peak and the second peak of total force observed in the intact model was compared with the results of experimental studies (Bedi et al., 2010; Gilbert et al., 2014). The absolute difference in width of the medial meniscus midbody between the OA model and the intact model obtained for the first peak was verified against the corresponding difference measured in a MRI study (Wenger et al., 2013) between the subjects with knee OA and reference subjects. Moreover, the medial translation of the medial meniscus at point N computed in the intact model was compared with values of meniscal extrusion determined in studies (Zhang et al., 2019; Li et al., 2020).

Outcome variables

Outcome variables were presented in the results section for the first and second peaks of total force corresponding to 25% and 81%, respectively, of the stance phase of the gait cycle (Fig. 2A). The total force was defined as the sum of the resultant forces in the medial and lateral compartments of the tibiofemoral joint. The von Mises stresses and hoop stresses were analyzed for the medial meniscus to investigate the location of potential meniscal tear initiation. The von Mises stress is the equivalent stress that allows the comparison of three dimensional stress state with the uniaxial yield stress. The hoop stresses are tensile stresses along the circumferential collagen fibers of menisci that are generated by axial forces in the knee joint. The hoop stresses may be compared with maximum tensile stress determined for the circumferential specimens. Strains were not reported in the results section, because they are proportional to stresses in the menisci. The meniscal deformation and tibiofemoral kinematics were presented to assess the influence of medial knee OA on the biomechanics of the medial meniscus. Curves for medial/lateral translation were not given because it was low in all models.

The hoop stresses and translations of the medial meniscus were measured at specific points and edges shown in Fig. 1E. The average and maximal hoop stresses were determined on the inner edge A of the anterior horn and on the inner edge P of the posterior horn. The meniscal translations were computed relative to the tibia. The posterior meniscal translation v_K was measured at point K and the translations u_K and u_N in the medial direction at points K and N, respectively. Point K was chosen on the inner edge of the medial meniscus body as the closest point to the location of the maximal hoop stress in the OA model. Point N was selected on the outer edge of the medial meniscus body in the location of maximal medial translation. Moreover, medial translation at point N is close to a meniscal extrusion of the medial meniscus body. The changes in medial meniscal width Δw_B and Δw_PH were calculated in the meniscal midbody and the posterior horn, respectively.

The relative and absolute changes in values of outcome variables for models G, M and OA were computed relative to the values obtained for the intact knee model. The maximal contact pressure between articular cartilage, the maximal von Mises stress and hoop stress in the medial meniscus, medial translation and change of the medial meniscus width, the tibial translations and rotations were compared between the intact model and the OA model to assess the influence of osteoarthritis on the biomechanics of the medial meniscus. Moreover, the differences in contribution of the medial meniscus in carrying of compressive load, the maximal hoop stress, posterior and medial translation of the medial meniscus were compared between models G and M to evaluate the influence of changes in geometry and material parameters on the specific results, respectively.

Stress state

Similar distribution of contact pressure on the surfaces of the tibial cartilage and menisci was observed in the intact knee model and the OA model for the first peak in Fig. 3. A 23% increase and a 29% decrease in the maximal value of the contact pressure between the articular cartilage was obtained in the medial and lateral compartment, respectively. The contour plots for the second peak (Fig. 3) show a 10% decrease in maximal contact pressure and anterior shift of its location on the medial tibial cartilage in the OA model. These changes are the effect of the greatest reduction in cartilage thickness at the central part of the medial tibial cartilage (Favre et al., 2017).

The contribution of the medial meniscus in carrying of compressive load increased from 48% to 53.7 in model G and decreased to 44.5% in model M for the first peak (see Table 5). Consequently, a minor reduction in forces between the tibial cartilage and the menisci was obtained in the OA model. The maximum von Mises stress of 22.5 MPa for the first peak was obtained in the anterior horn of the medial meniscus body in the intact model (Fig. 4). While, the maximum von Mises stress of 26.8 MPa was observed in the posterior part of the medial meniscus body in the OA model.

The maximal and average hoop stresses in the medial meniscus are compared in Fig. 5 and Table 6. The maximal hoop stress was 83% greater for the first peak and 58% greater for the second peak on edge P (see Fig. 1E) in the OA model compared to the intact model. Similar values of hoop stress were obtained for all models in the anterior horn of the medial meniscus for the first peak. Table 6 and Fig. 5 show that reduction in articular cartilage thickness (model G) caused an increase of approximately 70% in the maximal value and an increase of 100% in the average value of hoop stress on edge P for the first peak. Decrease in compressive stiffness (model M) resulted in lower maximal and average hoop stresses on both edges than in model G (Table 6). However, a considerably greater increase in hoop stress from 2.6 MPa to 12.1 MPa was observed near point K in model M than in model G (6.64 MPa) for the first peak. Approximately seven times higher hoop stress was obtained near point K in the OA model (Table 6). A substantial change in the distribution of hoop stress was observed in the sagittal cross section of the posterior horn of the medial meniscus in the OA model, see Fig. 6.

Meniscal deformation

An increase of 74% and 25% in medial(-) translation of the medial meniscus was obtained at point K and point N (compare Fig. 1E), respectively, in the OA model for the first peak (see Table 6). The greatest increase of 57% in the medial translation at point N was determined in model G. Considerably greater deformation in the medial direction was observed in the inner part of the medial meniscus body (Fig. 7) in models M and OA due to the reduction of the meniscal body width w_B by 1.13 mm and 1.88 mm, respectively. A decrease in width w_B was approximately 4 times greater in the OA model compared to the intact model (Table 6). The greatest reduction of the medial meniscal posterior horn width by 1.34 mm was observed in model M. The posterior translation at point K v_K decreased by 31% in model G and increased by 30% in model M (Table 6). The similar values of posterior translation v_K were observed in the intact model and the OA model. The differences in meniscal deformation between the models are generally lower for the second peak than for the first peak (Table 6).

Knee kinematics

The tibial translations and rotations during the stance phase of the gait cycle are presented in Fig. 8. The posterior translation of the tibia is smaller for almost the entire stance phase in the OA model than in the intact model (Fig. 8A). The reduction in articular cartilage thickness caused an increase in superior translation, internal rotation and varus rotation in models G and OA (Figs. 8B–8D). The greater increase in the varus rotation in the OA model was observed for the first peak (1.34°) than for the second peak (0.66°) (see Table 6).

Validation

The tibial translations and rotations obtained for the intact model and the OA model in the passive knee flexion simulation (see Fig. 9) were within range measured in the physical experiment (Germain et al., 2016). Only the medial translation of the tibia relative to the femur obtained in our study was beyond the physiological corridor from the previous study (Germain et al., 2016), but it was in the experimental range from an in vitro study (Belvedere et al., 2011). Moreover, tibial translation in the medial direction was reported in the majority of studies (see i.e., Li et al., 2007; Seisler & Sheehan, 2007; Belvedere et al., 2011).

Maximal tibial contact pressures of 5.9 MPa and 8 MPa obtained in the intact model for the first and second peaks, respectively, are consistent with values of 6 ± 0.5 MPa and 7.4 ± 0.6 MPa from the experimental study (Bedi et al., 2010). The anterior shift of the maximal contact pressure location on the medial tibial cartilage between the first peak and the second peak of total force (Fig. 3) was also observed in the previous studies (Bedi et al., 2010; Gilbert et al., 2014).

The medial meniscal width was 1.43 mm narrower in the OA model than in the intact model for the first peak, see Table 6. This value is in good agreement with absolute difference of 1.39 mm in meniscal width reported in the study (Wenger et al., 2013). The values of 1.92 mm and 2.72 mm of translation u_N for the first peak and second peaks, respectively, are consistent with values of 1.8 mm and 3.16 mm of the medial meniscus extrusion from the studies (Zhang et al., 2019; Li et al., 2020).