Prediction of local fixed charge density loss in cartilage following ACL injury and reconstruction: A computational proof‐of‐concept study with MRI follow‐up

The purpose of this proof‐of‐concept study was to develop three‐dimensional patient‐specific mechanobiological knee joint models to simulate alterations in the fixed charged density (FCD) around cartilage lesions during the stance phase of the walking gait. Two patients with anterior cruciate ligament (ACL) reconstructed knees were imaged at 1 and 3 years after surgery. The magnetic resonance imaging (MRI) data were used for segmenting the knee geometries, including the cartilage lesions. Based on these geometries, finite element (FE) models were developed. The gait of the patients was obtained using a motion capture system. Musculoskeletal modeling was utilized to calculate knee joint contact and lower extremity muscle forces for the FE models. Finally, a cartilage adaptation algorithm was implemented in both FE models. In the algorithm, it was assumed that excessive maximum shear and deviatoric strains (calculated as the combination of principal strains), and fluid velocity, are responsible for the FCD loss. Changes in the longitudinal T1ρ and T2 relaxation times were postulated to be related to changes in the cartilage composition and were compared with the numerical predictions. In patient 1 model, both the excessive fluid velocity and strain caused the FCD loss primarily near the cartilage lesion. T1ρ and T2 relaxation times increased during the follow‐up in the same location. In contrast, in patient 2 model, only the excessive fluid velocity led to a slight FCD loss near the lesion, where MRI parameters did not show evidence of alterations. Significance: This novel proof‐of‐concept study suggests mechanisms through which a local FCD loss might occur near cartilage lesions. In order to obtain statistical evidence for these findings, the method should be investigated with a larger cohort of subjects.


| INTRODUCTION
Anterior cruciate ligament (ACL) rupture is one of the most common traumatic knee joint injuries, which often involves damage to other tissues in the joint such as meniscus and cartilage. A rupture of the ACL mostly affects the young and healthy population causing not only pain and instability but it can also predispose the subject to posttraumatic osteoarthritis (PTOA). 1 Furthermore, there is a lack of evidence that ACL reconstruction is able to prevent the progression of PTOA. 2 The signs of PTOA include a loss of fixed charge density (FCD) of proteoglycans (PGs) from cartilage and surface disruption with lesions penetrating into the tissue. Consequently, cartilage swelling decreases around the lesion, reducing the ability of the collagen network to support tensile forces. However, the local changes in the composition of cartilage associated with these physiological changes remain unclear and are challenging to predict. 3 If a method was available and it would be capable of predicting these tissue changes with time, planning future activities of patients and/or interventions would become both more straightforward and more effective.
It has been suggested that local cartilage lesions might contribute to the development of PTOA following an ACL injury and reconstruction. 4 Computational knee joint modeling has been used to predict changes in the properties of articular cartilage due to abnormal loading. 5 Likewise, in vitro mechanobiological models have been shown capable of simulating the FCD loss around cartilage lesions as a function of time. 6 Suggested mechanisms leading to PTOA have been related to higher localized tissue strains (shear, deviatoric) or fluid flow around the lesions, leading to FCD loss. 7 In fact, it has been suggested that (a) FCD loss appears earlier than collagen damage 8,9 over a short follow-up time 10 and that (b) variations in the collagen fibril network organization are minor around cartilage defects. 11 Moreover, animal model studies have indicated that the collagen content does not change extensively in early PTOA, but rather follows other physical and compositional changes. 9,12 These findings form the basis for the analysis and prediction of cartilage FCD loss in injured joints. However, there are no studies that would have merged patient-specific knee joint models (including cartilage lesions) with adaptation algorithms in the quantitative and timedependent prediction of cartilage FCD loss in ACL injured and reconstructed knee joints.
Recent magnetic resonance imaging (MRI) studies have revealed that T 1ρ and T 2 relaxation times increase after the ACL injury and reconstruction surgery, suggesting that these are reflecting changes in cartilage composition. 13,14 However, there are no studies that have compared MRI follow-up information of ACL reconstructed patients with computational model predictions of FCD loss.
In this proof-of-concept investigation, we have applied an algorithm with (a) maximum shear strain, (b) deviatoric strain, and (c) fluid velocity-controlled tissue adaptation mechanisms into threedimensional (3D) patient-specific fibril-reinforced poroviscoelastic knee joint finite element (FE) models with swelling properties. The models are used to simulate alterations in the FCD content around cartilage lesions during the stance phase of the gait. The numerical predictions emerging from the model are compared to the longitudinal changes in T 1ρ and T 2 relaxation times from MRI at 1 to 3 years after ACL reconstruction. We hypothesized that the localized change of the FCD content around cartilage lesions in the ACL reconstructed knees would correspond to the increases in the relaxation times. It was postulated that an estimation of these compositional changes with the mechanistic knee joint model could improve the identification of lesions at a high risk of progression to PTOA and could be applied for simulating surgical interventions and rehabilitation procedures. Thus, the model could ultimately suggest an optimal intervention to slow down or prevent the progression of PTOA.

| Gait analysis and musculoskeletal modeling
The gait information of the patients was obtained at the 1-year follow-up time point using a published protocol. 19 The motion capture system consisted of 10 cameras (Vicon, Oxford Metrics) and 41 retroreflective markers to obtain segment position data.
Simultaneously, two in-ground force plates (AMTI, Watertown) were utilized to obtain ground reaction forces. Then, a generic muscu- Knee joint and quadriceps forces were simulated for every trial.
Then, the musculoskeletal model results were used to drive the knee joint FE models, similarly as has been done previously 21,22 (Figure 2). using eight-node hexahedral poroelastic (C3D8P) elements and modeled using fibril-reinforced poroviscoelastic and fibril-reinforced poroelastic materials, respectively, including swelling. 23 Specifically, the constitutive model assumes that the tissue is composed of solid and fluid matrices. The solid matrix is separated into a porous nonfibrillar part, representing the proteoglycan matrix, and a fibrillar network (viscoelastic in cartilage; elastic in meniscus), describing the collagen fibrils, and the influence of swelling caused by FCD. The total stress is given by

| Finite element knee joint model
where tot σ is the total stress tensor, s σ and fl σ represent the stress tensors of the solid and fluid matrices, respectively, p and ∆π are the hydrostatic and swelling pressures, respectively, I is the unit tensor, where E nf and nf ν are the Young´s modulus and the Poisson's ratio of the nonfibrillar matrix. The stresses in the viscoelastic collagen fibrils are defined with the damping coefficient η, the initial fibril network modulus E 0 , and the strain-dependent fibril network while the stresses in the elastic collagen fibrils are given by where f σ and f ε represent the fibril stress and strain, respectively. Therefore, collagen fibrils support primary tension. The fibril network stress emerges from the sum of primary and secondary collagen fibril OROZCO ET AL. stresses, which are computed individually for each integration point in each element. 23 Hence, stresses for these fibrils in tension are where i f, p σ and i f, s σ are the stresses for primary and secondary fibrils, respectively, z ρ is the relative collagen density, and C is the density ratio between primary and secondary fibrils. Then, the total stress tensor of the fibrillar network is defined as the sum of the stresses in where totf is the total number of fibrils, e i where cFCD 0 is the initial fixed charge density and nfl 0 is the initial fluid volume fraction. In addition, the chemical expansion stress T c is determined as where a 0 and κ are material constants 23 and − c is the mobile anion concentration. Moreover, the fluid flow in the nonfibrillar matrix is simulated according to Darcy's law where q is the flux in the nonfibrillar matrix, ∇p is the hydrostatic pressure gradient vector across the region, and k is the hydraulic permeability. The hydraulic permeability (k) is defined to be strain- where k 0 is the initial permeability, M is a positive constant, and e and e 0 are the current and initial void ratios, respectively. 24 The void ratio e is expressed by the ratio of the fluid to solid where n s is the solid volume fraction and n fl is the fluid volume fraction. The fluid velocity v fl through the porous medium is established as The depth-dependent primary collagen fibril orientations with split-line patterns, fluid fraction, and FCD distribution were implemented in the cartilage tissues. 25 In the menisci, the primary collagen fibrils were oriented circumferentially, and the fluid fraction and FCD content were homogeneously distributed. A complete list of the material parameters used is given in Table 1.
The meniscal roots were fixed to the bone using linear spring elements with the total stiffness of 350 N/mm at each root. 31 The quadriceps (QT) and patellar (PT) tendons, as well as collateral ligaments (MCL and LCL) and cruciate ligaments (ACL and PCL), were modeled using spring elements with a bilinear behavior. 32 These have recently been shown to provide reasonable FE modeling results, as compared to the models with solid ligaments. 33 The ligaments were assumed to be pre-elongated (MCL and LCL = 4%, ACL and PCL = 5%; of the initial length) at the segmented distance by using the bilinear spring option in Abaqus. The stiffness of the ligaments (MCL 100 N/mm, LCL 100 N/mm, ACL 380 N/mm, and PCL 200 N/mm,) and tendons (QT 475 N/mm and PT 545 N/mm) of the knee joint were selected from previous studies. 34,35 Similarly to previous studies, 22

| Cartilage adaptation algorithm
A previously presented iterative adaptation algorithm was utilized to predict changes in the FCD content of cartilage. 7 where max ε and min ε are the maximum and minimum principal strains, respectively.
The criteria regarding the strain thresholds were initially based on earlier studies. 6,36,37 A sensitivity analysis was also conducted (a strain threshold range was set from 0.15 to 0.6) to determine the correspondence between the simulation results and observed changes in the follow-up MRI maps. The initial fluid velocity threshold (0.03 mm/s) was chosen based on a previous study. 7 The sensitivity of the fluid velocity threshold to the modeling results was also tested (a threshold range was set from 0.01 to 0.2 mm/s). See the sensitivity analysis findings in Section 3. We also tested diverse linear and nonlinear approaches to study the rate of FCD loss. 7 where ε β is the strain-variable (either deviatoric (β = dev) or maximum shear strain (β = shr)), ,thres ε β is the strain threshold at which the cell death and the non-fibrillar matrix damage are assumed to initiate,

| Comparison between FCD loss predictions and MRI maps
The

| Sensitivity analysis for MRI relaxation time thresholds
The sensitivity analysis showed that when the relaxation time threshold was reduced from 60 to 50 ms, the volumes of cartilage with altered cartilage composition (as assumed from the altered MRI signal) reached the areas of healthy cartilage (based on WORMS evaluation included in Supporting Information Material and Figures S2 and S3).

| Sensitivity analysis for adaptation algorithm thresholds
The sensitivity analysis for the numerical threshold values revealed that with smaller strain and fluid flow velocity threshold values, the FCD loss was faster and occurred also in other locations than those near to the cartilage lesions (eg, deep cartilage and central tibial cartilage surface). In contrast, for larger values, the FCD loss was either prolonged or negligible in the entire cartilage volume ( Figures   S4 and S5).

| Patient-specific FE models
For patient 1, the FE model revealed those areas in the tibial cartilage where the maximum shear strain, deviatoric strain, and fluid velocity exceeded the thresholds of FCD loss and were at a maximum during the second peak of the stance phase (∼0.75 stance fraction) (Figure 4).
For patient 2, primarily, the maximum shear strain exceeded the FCD loss threshold in small areas located in both the lateral and medial tibial cartilage during the first peak of the stance phase (∼0.35 stance fraction) ( Figure 5). However, around the defect in the lateral femoral cartilage ( Figure 5C), only the fluid velocity exceeded the threshold in a small area during the loading response of the stance.

| Qualitative and quantitative comparison between T 1ρ and T 2 maps and mechanobiological model predictions
For patient 1, volumes of T 1ρ and T 2 relaxations times over the 60 ms threshold increased by ∼18% and ∼17% in the lateral tibial cartilage from 1-to 3-year follow-up time points, respectively ( Figures 6A and 7). For patient 2, the volumes of T 1ρ and T 2 relaxation times over the 60 ms threshold increased by ∼7% and ∼2% in the medial and lateral tibial cartilages from the 1-to 3-year follow-up time points, respectively ( Figures 8A and 9). Numerical predictions were consistent with the MRI results and showed that the volume of FCD loss was 3.5% and 2% in the lateral and medial tibial cartilages, respectively, when the model was driven by the maximum shear strain, while the deviatoric strain mechanism predicted a 0.5% FCD loss volume in the lateral tibial cartilage.
The FE model driven by the fluid flow velocity revealed a 0.2% FCD loss around the lesion in the femoral cartilage ( Figures 8B-D and 9).

| DISCUSSION
In the present proof-of-concept study, we applied a mechan-  Sixth, the FCD content was not allowed to decrease to zero during the iterative process. The minimum value allowed was set to 10% of the smallest initial FCD content to avoid computational instabilities. The FCD distribution reached an equilibrium state after the iterative process. This might not be fully realistic since a possible further FCD loss might occur due to other factors that were not considered in the model, such as different physical activities, thresholds for the initiation of FCD loss, and biochemically driven degradation due to diffusion of inflammatory cytokines. However, our numerical predictions led to a good match with the experimental MRI follow-up data.
Due to the resolution of MRI, probably small uncertainties were present during the segmentation process. The voxel size was 0.55 × 0.55 × 4 mm which might mean that synovial fluid has contributed to the slightly increased relaxation times due to the partial volume effect, especially for patient 2. However, since the changes were rather small, this uncertainty should not affect our conclusions.
For both patients, full-thickness cartilage lesions (WORMS score 2.5) were diagnosed at the 1-year follow-up time point by three experienced musculoskeletal radiologists. However, the cartilage lesions in both patients were challenging to identify and segment due to the size of the lesions and inherent limitations of MRI (ie, the contrast between fluid and cartilage, spatial resolution). In particular, the lesion segmentation of patient 2 was challenging, and the final lesion geometry appears to be ideally symmetric. However, a musculoskeletal radiologist assisted us in identifying and segmenting the lesion geometries. It is also worth mentioning that our model predictions and T 1ρ and T 2 relaxation times did not reveal any evidence of alterations around the lesion in patient 2. Thus, even if the segmentation could be performed more accurately, the main conclusions would not change due to the small lesion size.
In conclusion, we suggest that the FCD decrease following ACL injury and reconstruction, including cartilage injury and subsequent tissue loading, might be caused by a large tissue deformation around the defect and extensive leakage of PGs through the damaged surface by high fluid outflow. 7 In the future, we will expand this proofof-concept study to investigate also other mechanisms (biomechanical and biochemical) leading to changes in cartilage composition and structure after a traumatic joint injury, and compare the model results with the data obtained from a larger cohort of patients. After careful validation, the model could be applied in the planning of joint F I G U R E 7 Relative volumes of FCD loss, as predicted by each mechanism in the FE models, and altered T 1ρ and T 2 relaxation times during the follow-up of patient 1. The relative volume of FCD loss in the FE models was estimated as the volume of the elements associated with FCD loss in the last iteration, divided by the total volume of each compartment. The relative volume of cartilage with altered relaxation times was estimated by subtracting the calculated volume over the threshold of 60 ms at