Front. Bioeng. Biotechnol. Frontiers in Bioengineering and Biotechnology Front. Bioeng. Biotechnol. 2296-4185 Frontiers Media S.A. 1214693 10.3389/fbioe.2023.1214693 Bioengineering and Biotechnology Original Research Signatures of disease progression in knee osteoarthritis: insights from an integrated multi-scale modeling approach, a proof of concept Mohout et al. 10.3389/fbioe.2023.1214693 Mohout Ikram 1 * Elahi Seyed Ali 1 2 Esrafilian Amir 3 Killen Bryce A. 1 Korhonen Rami K. 3 Verschueren Sabine 4 Jonkers Ilse 1 1 Department of Movement Science, Human Movement Biomechanics Research Group, Leuven, Belgium 2 Mechanical Engineering Department, Soft Tissue Biomechanics Group, Leuven, Belgium 3 Department of Technical Physics, Biophysics of Bone and Cartilage Research Group, University of Eastern Finland, Kuopio, Finland 4 Department of Rehabilitation Science, Research Group for Musculoskeletal Rehabilitation, Leuven, Belgium

Edited by: Andrew T. M. Phillips, Imperial College London, United Kingdom

Reviewed by: Andreas Martin Seitz, Ulm University Medical Center, Germany

Pinaki Bhattacharya, The University of Sheffield, United Kingdom

 *Correspondence: Ikram Mohout, ikram.mohout@kuleuven.be 27 07 2023 2023 11 1214693 30 04 2023 13 07 2023 Copyright © 2023 Mohout, Elahi, Esrafilian, Killen, Korhonen, Verschueren and Jonkers. 2023 Mohout, Elahi, Esrafilian, Killen, Korhonen, Verschueren and Jonkers

This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) and the copyright owner(s) are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.

Introduction: Knee osteoarthritis (KOA) is characterized by articular cartilage degeneration. It has been widely accepted that the mechanical joint environment plays a significant role in the onset and progression of this disease. In silico models have been used to study the interplay between mechanical loading and cartilage degeneration, hereby relying mainly on two key mechanoregulatory factors indicative of collagen degradation and proteoglycans depletion. These factors are the strain in collagen fibril direction (SFD) and maximum shear strain (MSS) respectively.

Methods: In this study, a multi-scale in silico modeling approach was used based on a synergy between musculoskeletal and finite element modeling to evaluate the SFD and MSS. These strains were evaluated during gait based on subject-specific gait analysis data collected at baseline (before a 2-year follow-up) for a healthy and progressive early-stage KOA subject with similar demographics.

Results: The results show that both SFD and MSS factors allowed distinguishing between a healthy subject and a KOA subject, showing progression at 2 years follow-up, at the instance of peak contact force as well as during the stance phase of the gait cycle. At the peak of the stance phase, the SFD were found to be more elevated in the KOA patient with the median being 0.82% higher in the lateral and 0.4% higher in the medial compartment of the tibial cartilage compared to the healthy subject. Similarly, for the MSS, the median strains were found to be 3.6% higher in the lateral and 0.7% higher in the medial tibial compartment of the KOA patient compared to the healthy subject. Based on these intersubject SFD and MSS differences, we were additionally able to identify that the tibial compartment of the KOA subject at risk of progression.

Conclusion/discussion: We confirmed the mechanoregulatory factors as potential biomarkers to discriminate patients at risk of disease progression. Future studies should evaluate the sensitivity of the mechanoregulatory factors calculated based on this multi-scale modeling workflow in larger patient and control cohorts.

 osteoarthritis articular cartilage finite element modeling musculoskeletal modeling adaptive modeling collagen fibril degeneration proteoglycans depletion multiscale modeling section-at-acceptance Biomechanics 1 Introduction

Osteoarthritis (OA) is the most common joint disease, with the knee being the most frequently affected joint. This complex multi-factorial disease is characterized by degenerative changes in articular cartilage, subchondral bone and by synovial inflammation (Baliunas et al., 2002; Martel-Pelletier et al., 2016; Syed and Makris, 2020; Katz et al., 2021). Mechanical degeneration of articular cartilage is a hallmark of OA and a major contributor to joint dysfunction (Sinusas, 2012; Katz et al., 2021). Therefore, understanding biomechanical factors involved in cartilage degeneration is critical for developing effective therapeutic strategies to prevent or treat this debilitating disease.

Articular cartilage, an avascular connective tissue, is composed of specialized chondrocyte cells embedded in an extracellular matrix (ECM), which consists mainly of water, proteoglycans (PGs) and a collagen fiber network. The interactions between these constituents dictate the mechanical behavior of cartilage and enable it to distribute mechanical load across the joint surface (Sophia Fox et al., 2009; Heinegård and Saxne, 2011; Mukherjee et al., 2020). However, changes in mechanical loading during locomotion will alter stresses and strains within the tissue, which—upon exceeding a critical threshold—is associated with irreversible changes to the collagen network and depletion of the PG content (Speirs et al., 2014; Tanska et al., 2018; Elahi et al., 2021; Elahi et al., 2023). Indeed, the mechanical joint environment is considered to be a significant contributor to the disease process, with pathomechanical loading potentially leading to maladaptive responses of chondrocytes and ECM, as well as increased inflammation, further aggravating cartilage tissue degeneration (Hurwitz et al., 2000; Baliunas et al., 2002; Sun, 2010).

To non-invasively study the mechanical joint environment, in silico models have been developed. In particular, musculoskeletal (MSK) models have proven useful in estimating knee joint loading during locomotion, which has been found to be altered in patients with OA (Lloyd and Besier, 2003; Lenhart et al., 2015; Meireles et al., 2017). By combining MSK models with finite element (FE) models of the articular cartilage, it is possible to simulate the effects of mechanical joint loading on cartilage strain distributions (Adouni and Shirazi-Adl, 2014; Halonen et al., 2014; Halonen et al., 2017; Erdemir et al., 2019; Esrafilian et al., 2020; Esrafilian et al., 2021). In addition, articular cartilage FE models can be further enhanced by assigning complex material models that mimic the main constituents of cartilage ECM, allowing an even more detailed evaluation of cartilage tissue responses to mechanical loading (Wilson et al., 2005). Such in silico models provide a valuable platform for investigating the complex interplay between mechanical loading, cartilage structure and function, and constituent damage and depletion.

Cartilage degeneration algorithms, also referred to as adaptive algorithms have emerged that mimic time-dependent, OA-associated cartilage degeneration processes (Mononen et al., 2018; Orozco et al., 2018; 2021; Eskelinen et al., 2019; Elahi et al., 2021). The Cartilage Adaptive Reorientation Degeneration (CARED) model, namely, the most recently developed integrated adaptive model, uses a fibril-reinforced poroviscoelastic material model to calculate the strain in collagen fibrils and maximum shear strain as mechanoregulatory parameters that drive collagen degradation and PG depletion, respectively (Elahi et al., 2021; Elahi et al., 2023). Although the CARED model is able to predict degenerative changes in an explant model (osteochondral plug with 3 mm diameter) under different loading conditions, this is not representative for cartilage degeneration in the whole knee joint. To extend the insights on the impact of mechanical loading to cartilage degeneration in vivo, we need to define a multi-scale modeling workflow that includes a whole joint to tissue analysis and that relies on a complex cartilage material model to allow evaluating the impact of the mechanical tissue environment on constituent degeneration.

Whole knee joint FE models in combination with adaptive models have been developed (Mononen et al., 2018; Orozco et al., 2021) to predict collagen fibril degradation in the native knee joint (Mononen et al., 2016) and, more recently, to predict proteoglycan depletion and collagen fibril re-orientation (Heck et al., 2015; Orozco et al., 2018). To the best of our knowledge such algorithms have not been used to study cartilage degenerative processes specifically in early-stage OA cohorts. As such, insight into mechanoregulatory factors that predispose to cartilage degeneration in early disease stages remain undocumented.

This study proposes a unique, multi-scale, patient-specific in silico approach to evaluate the mechanical environment of the human knee joint based on loading conditions derived from MSK modeling during gait in a healthy and an early-stage knee OA (KOA) subject displaying (fast) progression over a 2-year time frame. These subject-specific loading and boundary conditions were estimated based on experimental gait analysis data collected at baseline, prior to a 2-year follow-up, which is processed through a MSK modeling workflow. The proposed modeling pipeline in this study uses a unique MSK modeling-informed approach, where the subject-specific boundary and loading conditions are used to drive an FE model that is then used to evaluate the mechanical environment of cartilage and meniscus tissues within the knee joint. Firstly, mechanoregulatory factors known to underly cartilage constituent damage were evaluated, at baseline, to characterize the mechanical tissue environment prone to cartilage degeneration. More specifically, strain in fibril direction (SFD) and maximum shear strain (MSS) in the articular cartilage were evaluated as mechanoregulatory factors driving collagen degradation and PG depletion, respectively. Secondly, tissue responses in the menisci, a structure known to participate in strain mitigation and protection of the underlying cartilage tissue (Englund, 2009; Netravali et al., 2010; Choi et al., 2014; Verdonk et al., 2016), were analyzed to evaluate their role in altering the mechanical environment of the cartilage tissue in the presence of subject-specific loading conditions. We hypothesize that under subject-specific loading conditions and consistent use of generic knee articular geometries in both MSK and FE modeling workflows, the mechanoregulatory factors are distinct between a healthy and a KOA subject and that the menisci are differentially implicated in load transmission to the cartilage surface. If confirmed, we anticipate that this study will provide proof of concept data showing the potential of these mechanical factors as functional biomarkers to discriminate patients at risk of (fast) disease progression.

 2 Methods

The study workflow, illustrated in Figure 1, consisted of consecutive MSK and FE modeling pipelines. The red arrows indicate two unique methodological features of this study: 1) consistency in the use of the knee articular geometries in both MSK and FE modeling pipelines and 2) MSK modeling-informed subject-specific loading and boundary conditions as inputs to the developed FE model.

Workflow of the study: (A–C) represent the MSK modeling pipeline to estimate subject-specific knee joint kinematics and kinetics (D,E) represent the FE modeling pipeline to estimate cartilage tissue responses. Red arrows indicate unique approaches: consistent knee geometries in MS and FE pipelines, and MSK-workflow derived subject-specific loading and boundary conditions as inputs for FE model.

One healthy control subject and one progressing early KOA subject with similar demographics (Table 1) were selected from a longitudinal historical dataset (Baert et al., 2013; Mahmoudian et al., 2016). The selection of the KOA subjects was based on the classification criteria of Luyten et al. (2012): More specific, a subject with a minimum increase of 1 Kellgren-Lawrence (KL) score (Kellgren and Lawrence, 1957) - evaluated by the same rater blinded for disease or follow-up status - over a two-year follow-up. The patient had no history of prior trauma or surgery with regards to the lower limbs and back. Furthermore, the patient did not have any musculoskeletal disorders (other than KOA) nor did they suffer from any neurological diseases that would affect their balance or coordination during gait during a six-months period preceding the data collection. All procedures were approved by the local ethical committee (S50534) and written informed consent was obtained from the participants.

Patient demographics *Kellgren-Lawrence classification score for a right knee lateral compartment progressor at baseline and after a 2-year follow-up.

Subject type Gender Age (years) Weight (kg) Height (m) Body mass index (BMI) (kg/m2) KL-score baseline* KL-score follow-up*
medial lateral medial lateral
Control Female 72 61 1.56 25 0 0 0 0
KOA Female 70 68 1.56 28 1 1 1 2
2.1 MS-based estimation of patient-specific knee joint kinematics and kinetics 2.1.1 Experimental data collection

Experimental gait analysis data from the aforementioned longitudinal historical dataset (Baert et al., 2013; Mahmoudian et al., 2016; Meireles et al., 2017) was analyzed at baseline (before the two-year follow-up). This included three-dimensional marker position data acquired with a motion capture system (10MX Vicon, 100 Hz) and ground reaction force data collected with in-ground force plates (AMTI, 1000 Hz). First, a static calibration trial was performed (Cappozzo et al., 1995), followed by four trials of over-ground walking at participant’s self-selected speed. The data were then converted using a custom Matlab script (MATLAB R2020b, The Math Works, Inc., Natick, Massachusetts, United States ) as a pre-processing step for the MSK workflow.

 2.1.2 Musculoskeletal modeling

The OpenSim Joint Articular Mechanics (JAM) (Smith et al., 2014) modeling workflow was used to estimate knee joint kinematics, moments and loading. A 3D MSK model with 6 degrees of freedom (DoF) tibiofemoral (TF) and patellofemoral joint kinematics (Lenhart et al., 2015) was used and included the major knee ligaments which were modeled as non-linear spring bundles (Blankevoort and Huiskes, 1991). The state-of-the-art OpenSim-JAM workflow uses standard OpenSim tools (Delp et al., 2007) in combination with a unique knee joint contact model capable of estimating joint secondary coordinates (Anterior-posterior and medial-lateral translations, internal-external and adduction-abduction rotation). As such, the menisci are not explicitly modeled but are included in the parameters of the knee contact model.

The experimental gait analysis data is processed through the MSK OpenSim-JAM workflow. First, the generic model is scaled to the subject’s anthropometry, based on experimental marker data from the static trails. Then an inverse kinematics algorithm was used to estimate joint angles (i.e., primary kinematics) based on a minimization of marker position between the model and experimental data (Lu and O’Connor, 1999). The concurrent optimization of muscle activations and kinematics (COMAK) algorithm (Smith et al., 2018), as part of OpenSim-JAM, was used to simulate the desired secondary tibiofemoral kinematics and contact forces.

The following specific outputs of the OpenSim-JAM were used as boundary conditions for the FE model: TF anterior-posterior translation, TF medial-lateral translation, TF internal-external rotation, knee flexion angle, varus-valgus moment and superior-inferior knee joint contact force. One representative trial for each subject was processed, more details on the outcome of the MS simulations are presented in Supplementary Material.

 2.2 Finite element modeling of cartilage tissue responses 2.2.1 Geometries and mesh

FE models of the human knee joint were created using Abaqus software (Dassault Systèmes, United States), that includes the tibial and femoral cartilage surfaces, (part of the) tibia bone and the menisci based on previous FE models developed by Esrafilian et al. (2020), Halonen et al. (2016) and Bolcos et al. (2018) which were validated against experimental data (Halonen et al., 2017; Esrafilian et al., 2020; Esrafilian, 2021; Esrafilian et al., 2022). The patella was not included in this model as this study focusses on the tibiofemoral cartilage tissue responses.

Identical generic knee joint geometries were used in the FE model and the MSK knee contact model to enforce geometrical consistency between the two workflows. Hexahedral meshes, generated in ANSA (v21.0.1, BETA CAE Systems International AG, Switzerland), were based on the meshes of the Lenhart et al. (2015) knee contact model that were derived from MRI segmentations of a healthy subject. Given that the KOA patient presented with KL score 1 indicative of early OA, no structural changes in the cartilage such as cartilage thinning were assumed to be present. Hence the use of a generic mesh was deemed acceptable. Despite not being incorporated in the MSK knee contact model, the designated mesh of the menisci was taken from the same healthy subject and was used to generate the hexahedral mesh for the FE model. Soft tissues (apart from the ligaments) were modeled as porous 8-node elements were used (C3D8P element type), with a total of 49,810 elements (femoral cartilage = 8,880 elements, menisci = 7,680 elements, tibial cartilage = 33,250 elements). Figure 2 shows the hexahedral meshes of the tibial cartilage, femoral cartilage and menisci. The Tibia bone was modeled as 675 8-node elements (C3D8 element type). A mesh convergence test was previously performed by Halonen et al. (2016) to verify the mesh size and number of elements in the model. To further enhance the model convergence, the superficial element layers of tibial cartilage have been subdivided resulting in a four-layered mesh with two finer superficial layers and two more coarse bottom layers as suggested by Halonen et al. (2016). Furthermore, the FE model included the anterior cruciate ligament (ACL), posterior cruciate ligament (PCL), lateral collateral ligament (LCL) and medial collateral ligament (MCL), using identical insertion points, origins and stiffness as in the MSK model. The ligaments were modeled as non-linear springs (SpringA element type, (Blankevoort and Huiskes, 1991)).

The hexahedral mesh of the soft tissues in the model: tibial cartilage (33,250 elements), femoral cartilage (8,880 elements) and menisci (7,680 elements). All soft tissue parts consist of four layers. The tibial cartilage mesh consists of two finer superficial layers and two coarse bottom layers to facilitate convergence (Halonen et al., 2016). The x-, y-, z-directions represent the anterior-posterior, superior-inferior and medial-lateral direction, respectively.

 2.2.2 Material model

Cartilage was modeled as a fibril-reinforced poroviscoelastic (FRPVE) material including a biphasic (fluid and solid phases) description with a hyperelastic Neo-Hookean solid matrix representing the PGs and a viscoelastic fibril network representing collagen fibrils (Wilson et al., 2004; 2005; Julkunen et al., 2007). The primary collagen fibrils followed the depth-dependent Benninghoff-type arcade architecture and the superficial collagen orientation was indicated by split lines towards the center of area (See Supplementary Figure S2 in the Supplementary Material). Furthermore, the fibril density and the fluid fraction also varied throughout the depth of the tissue (M. E. Mononen et al., 2012; Esrafilian, et al., 2020; Benninghoff, 1925; Goodwin et al., 2004). The effect of these depth-wise material properties on cartilage tissue stresses and strains have been analyzed in a previous study by Halonen et al. (2013). For the menisci, a fibril-reinforced poroelastic (FRPE) material was used, which is analogous to the FRPVE material model apart from the collagen fibrils which were modeled as an elastic fibril network. A linear elastic material model was used for tibia bone and the material parameters of the ligaments were taken from the MSK model (Blankevoort and Huiskes, 1991; Currey, 2004; Lenhart et al., 2015). A more detailed description of the material model and specific material parameters (Wilson et al., 2004; Julkunen et al., 2007; Makris et al., 2011; Dabiri and Li, 2013) is provided in the Supplementary Material. The same material model and material parameters were used for both subjects. Considering the KOA patient’s KL score of 1, indicating early OA, this assumption was considered acceptable.

 2.2.3 Subject-specific loading and boundary conditions

The time-dependent loading and boundary conditions derived from the MSK workflow were applied throughout the stance phase of the gait cycle. This was done by means of a transient analysis, more specifically a soils consolidation analysis, with an implicit FE scheme. The following subject-specific loading and boundary conditions derived from the MS model were implemented at the reference point of the femur in the FE model: Anterior-posterior translation, medial-lateral translation, internal-external rotation, knee flexion angle, varus-valgus moment and the superior-inferior TF joint contact force. This reference point was located between the medial and lateral femoral epicondyles and was tied to the femoral cartilage-bone interface. The origins of the ligaments were also tied to this reference point. The FE model was driven at this femoral reference point and thus the output of the MS model was extracted from the local coordinate system of the proximal tibia as this corresponds to the relative motion of the femur to the tibia.

Femur bone was considered rigid and its effect was included in the FE model with encaster boundary condition on the femoral cartilage internal surface. The tibia bone was modeled as linear elastic for the purpose of facilitating the convergence of the model to a numerical solution. The bottom surface of the tibial bone was encastred. The following frictionless contact types were defined in the model as nodes to surfaces: Tibial cartilage to femoral cartilage, menisci to femoral cartilage, tibial cartilage to menisci and menisci to tibia bone. The assumption of frictionless contact is widely accepted as the standard in knee joint FE modeling, primarily because the main function of articular cartilage is to provide a smooth and frictionless surface for joint articulation and efficient load distribution across the joint surface (Mukherjee et al., 2020).

 2.3 Analysis of FE model results 2.3.1 Evaluation of cartilage tissue response

This study analyzed maximum shear strain (MSS, see Eq. 1) and strain in fibril direction (SFD, see Eq. 2). Thresholds of MSS = 30% and SFD = 10% were chosen as being indicative of PG depletion and collagen degeneration respectively, as implemented in the most recently developed CARED model (Elahi et al., 2021) and based on literature findings (Eskelinen et al., 2019). ε max = max   ɛ p , 1 − ɛ p , 2 , ɛ p , 1 − ɛ p , 3 , ɛ p , 2 − ɛ p , 3 , where ɛp,1, ɛp,2, ɛp,3 are the principal strains of the Green-Lagrangian strain tensor. ε f = ln   F e f , where F is the deformation gradient tensor and ef is the unit vector of the fibril orientation.

First, tibial cartilage tissue responses were evaluated throughout the entire stance phase of the gait cycle. Here, the analysis was restricted to the tibial cartilage volume related to the contact area, this area amounts to 20% of the total tibial cartilage volume as observed in the model, as tissue strains (mostly) outside the contact area between femoral and tibial cartilage and between menisci and tibial cartilage were considerably low and did not contribute to the analysis. The strains were extracted from the centroids of the tibial cartilage elements and element volumes were taken into account. Then, tibial cartilage tissue responses were evaluated at peak superior-inferior knee contact force, as tissue strains reached their absolute maxima at this time-point during the stance phase. Results were evaluated qualitatively based on strain responses at the tibial cartilage surface and quantitatively by means of violin plots of the tissue strains within the cartilage volume. This quantitative analysis was also restricted to the tibial cartilage volume related to the contact area.

 2.3.2 Evaluation of meniscal tissue response

The maximum principal strain (MPS) of the logarithmic strain tensor was analyzed to reflect the overall tissue response of the menisci. The MPS was not directly associated with structural changes in the menisci. The tissue response was evaluated at peak superior-inferior knee contact force to be consistent with the evaluation of the tibial cartilage. Analogous to the evaluation of the tibial cartilage, a qualitative evaluation was done based on the strain at the (superior) surface of the menisci and a quantitative analysis was performed by means of violin plots of the tissue responses within the menisci. Here, the whole volume of the menisci was included in the quantitative analysis as (nearly) the whole surface of the menisci was in contact with either the femoral or tibial cartilage.

 3 Results 3.1 Evaluation of the cartilage tissue response 3.1.1 Strain evolution throughout the stance phase of the gait cycle

During the majority of the stance phase, the MSS (Figure 3A) and SFD (Figure 3B) were higher in the KOA subject than in the control subject. However, neither the median nor the 75th percentile boundary exceeded the degeneration thresholds of MSS and SFD.

Maximum shear strain (A) and strain in fibril direction (B) evaluated throughout the stance phase in both tibial compartment volumes. The graphs represent the median strain and the lower and upper boundaries of the shaded areas represent the 25th and 75th percentile, respectively.

 3.1.2 Strain in fibril direction at peak contact force

For both lateral and medial compartments, the strain in fibril direction (SFD) was more elevated in the KOA compared to the control subject, with the median being 0.82% higher in the lateral and 0.4% higher in the medial compartment respectively (Figure 4A). With the exception of some local strain concentrations, overall, the SFD thus did not exceed the threshold of collagen fibril degeneration (10%). Figure 4B shows a posterior shift in center-of-pressure (COP) in the KOA subject of 2.9 mm in the lateral compartment and 5.1 mm in the medial compartment compared to the control subject.

Strains in fibril direction (SFD) evaluated at peak contact force: (A) Violin plot (depicting the 25th percentile (first quartile), median (second quartile), 75th percentile (third quartile), density trace (smoothed histogram) and whiskers (based on 1.5 interquartile range value)) of the SFD (%) of the medial (top) and lateral (bottom) volumes of the tibial compartment for the control (blue) and KOA (red) subject. (B) SFD (%) in the lateral (left) and medial (right) tibial cartilage of the control (top) and progressor (bottom) subjects. The local strain concentration in the control subject is indicated (B) which corresponds to the bulge at the upper tale of the violin plot (A) of the medial compartment in the control subject. COP indicates the center-of-pressure in each tibial compartment.

 3.1.3 Maximum shear strain at peak contact force

For both compartments, the MSS of the KOA subject was more elevated compared to the control subject (Figure 5), with the median being 3.6% higher in the lateral tibial compartment and 0.7% in the medial tibial compartment respectively (Figure 5A). The MSS did surpass the threshold of degeneration (above 30%) and a larger volume was affected in the KOA subject (Figure 5A).

Maximum shear strain (MSS) evaluated at peak contact force: (A) Violin plot (depicting the 25th percentile (first quartile), median (second quartile), 75th percentile (third quartile), density trace (smoothed histogram) and whiskers (based on 1.5 interquartile range value)) of the MSS (%) of the medial (top) and lateral (bottom) volumes of the tibial compartment for the control (blue) and KOA (red) subject. (B) MSS (%) in the lateral (left) and medial (right) tibial cartilage of the control (top) and progressor (bottom) subject. COP indicates the center-of-pressure in each tibial compartment.

 3.2 Evaluation of the meniscal tissue response: maximum principal strain at peak contact force

The maximum principal strain (MPS) is more elevated in the KOA subject compared to the control subject (Figure 6) The median MPS of the KOA subject is 2.5% and 3.1% higher for lateral and medial meniscus compared to the control, respectively.

Maximum principal strain (MPS) evaluated at peak contact force: (A) Violin plot (depicting the 25th percentile (first quartile), median (second quartile), 75th percentile (third quartile), density trace (smoothed histogram) and whiskers (based on 1.5 interquartile range value)) of the MPS (%) of the medial (top) and lateral (bottom) menisci volumes for the control (blue) and KOA (red) subject. (B) MPS (%) in the lateral (left) and medial (right) menisci of the control (top) and progressor (bottom) subject.

 4 Discussion

This study is the first to compare key mechanoregulatory factors associated with cartilage tissue degeneration between a control subject and an early-stage KOA patient (KL 1) showing disease progression in 2 years. The study focused on two mechanoregulatory factors, namely, strain in the fibril direction (SFD) and maximum shear strain (MSS) and evaluation was performed using subject-specific loading conditions determined during the baseline stance phase of gait.

Using the proposed multi-scale patient-specific in silico approach, the local cartilage strains in terms of SFD and MSS, were found to be higher for the KOA patient compared to the control subject throughout the majority of stance duration and at peak contact force. This shows how these intra-tissue strains during locomotion are distinct between the healthy and early KOA subject already at baseline despite the subjects’ similar demographics. During stance, this distinction is even more clear when evaluating the 75th percentile of the tissue strains that were elevated throughout the gait cycle in the OA subject. This is in line with previous suggestions that not only instantaneous loading that surpasses the threshold of cartilage degeneration, but also the accumulated load over time (cumulative load) during normal ambulatory activities, could induce cartilage degeneration (Mononen et al., 2016). Furthermore, the strain difference at peak contact force between subjects is most prominent in the lateral compartment where after 2 years follow-up, structural progression can be observed in the KOA subject. This could invigorate our hypotheses that the mechanoregulatory factors studied here may serve as functional biomarkers with the potential to discriminate patients at risk of disease progression already in the early phases of the disease process. Indeed, the KOA subject presented with structural lateral compartment OA progression after two-year follow-up as evidenced by an increase in KL score (KL score 1 at baseline to KL score 2 after two-year follow-up).

Importantly, not only the mechanical environment of the cartilage surface is altered, also the maximum principal strain in the menisci is consistently higher in the KOA compared to the control subject, which is likely also due to the posterior shift in the center-of-pressure. This suggests an important role for the meniscus in strain mitigation in early KOA subjects. Any deficiency in meniscal integrity could therefore even further amplify the excessive articular cartilage loading. This is in line with literature where (partial) meniscectomy and meniscal subluxation have been associated with the degeneration process of cartilage (Choi et al., 2014; Verdonk et al., 2016). It is important to note that only one KOA subject was included in this proof-of-concept study. Hence, these findings will need to be validated in larger cohorts.

The altered joint level loading during gait in the patient with KOA impacts the local cartilage micro-environment such that constituent depletion occurs. Indeed, although the SFD was well below the threshold of degeneration (10%), MSS surpassed the threshold (30%) in 3.3% of the cartilage volume of the KOA subject. This would suggest a potential risk of PG depletion to drive the OA progression rather than collagen degeneration, which is highly relevant as it has been suggested that damaged cartilage can suffer from PG depletion under physiological loading conditions (Orozco et al., 2018) and that these changes in PG content can still be reversible which is thought not to be the case for changes in the collagen network (Elahi et al., 2023). An unexpected finding of this study was that 1.8% of the tibial cartilage volume of the control subject also surpassed the MSS threshold of degeneration, suggesting PG depletion. This could be attributed to age or possibly computational limitations as local strain concentrations in the models were found specifically when the load-bearing contact location was at the center of the respective tibial compartments. This is more prominent in the control subject as a posterior shift in contact location was found in the KOA subject.

Finally, in this study, personalized loading conditions were estimated from MSK modeling and used to drive the FE model to study the tissue strain responses. These personalized loading conditions estimated during gait differed between the control subject and the KOA patient, with higher inferior-superior contact force and knee adduction moment, as well as higher knee flexion in the KOA compared to the control subject. Previous studies have shown that the aforementioned characteristics are distinctive for KOA subjects (Hurwitz et al., 2000; Kaufman et al., 2001; Al-Zahrani and Bakheit, 2002; Baliunas et al., 2002; Astephen et al., 2008). Additionally, the posterior shift in center-of-pressure is in-line with previous studies that reported a posterior-lateral shift in center-of-pressure and load-bearing area when comparing KOA subjects to control subjects (Andriacchi and Mündermann, 2006; Meireles et al., 2017). These findings highlight the need to use disease-specific kinematic-kinetic input to the whole knee FE and adaptive workflow modeling, for instance using integrated 3D gait analysis data is used within an integrated MSK and FE modeling workflow to fully capture the altered mechanical environment in the cartilage tissue in the homeostatic joint and a joint susceptible for OA disease progression.

4.1 Limitations and future work

In this study, proof of concept data was obtained from two participants which could serve as a stepping stone for future analysis where more participants are included. It would be particularly interesting to include non-progressive KOA subjects as this would provide more information and likely more insight into the relation between localized tissue strains and cartilage degeneration. It must be noted that executing the FE workflow is a time-consuming task with high computational costs, and it can be challenging to achieve convergence of the model. The advantage of the implemented modeling approach is that it uses generic geometries consistent with the MSK pipeline. This alleviates the need for re-meshing of the geometries, which saves a significant amount of time, however at the expense of losing personalized details on the articular cartilage geometry. Thus, variations in cartilage thickness are not accounted for. Although the use of a generic mesh was deemed acceptable in this study, as no cartilage thinning was assumed to be present in the KOA patient with KL score 1, personalized geometry, in particular thickness, can further enhance the accuracy of strain predictions within the cartilages and menisci.

Although subject-specific loading conditions were applied, further personalization of the model in terms of cartilage thickness and geometry, knee contact with friction or tissue material properties in the FRPVE material model based on the disease state of the subject (healthy, early OA, late-stage OA) needs to be considered and might affect the SFD and MSS measures. Although state-of-the-art, the use of similar material parameters for a healthy and KOA subject is a limitation as it can alter the tissue mechanics significantly. In future work, experimental non-invasive characterization of subject-specific material parameters (e.g., using Raman spectroscopy or MRI) may allow more fine-tuning of the estimation of mechanical parameters driving the mechanoregulatory algorithm. However, in this study, given the consistent implementation of the entire modeling workflow in both subjects and the inclusion of subjects with similar demographics, we are confident in the comparability of the FE model results and the validity of the conclusions regarding the relative differences between the control and KOA subjects. Furthermore, Figure 3, demonstrates that the strains in the KOA subject are higher than those in the control subject for the vast majority of the stance phase indicating a noticeable disparity between the FE model results of the two subjects. As a next step, an additional uncertainty analysis and validation, particularly in a larger cohort of subjects would need to be performed to ensure the accurate prediction of the knee joint’s mechanical environment for each individual.

Altered cartilage tissue response, namely, local strain concentrations, seem to be present at the origin of the collagen split lines where there is a discontinuity in the fibrils. This could be an underlying cause of volume surpassing the MSS threshold in the control subject, suggesting PG depletion. These local strain concentrations can be considered and compensated for in adaptive modeling workflows. Studies by Párraga Quiroga et al. (2017) and Seyed and Makris. (2020) have proposed a non-localization theory, where strains are averaged over a certain volume. This would change the interpretation of these results regarding PG depletion in adaptive modeling workflows.

Additionally, there is an ambiguity regarding the exact thresholds of degeneration, with different studies reporting alternative thresholds as high as 50% for the MSS (Orozco et al., 2018; Orozco et al., 2021). In future studies, analysis of SFD and MSS during gait in a larger cohort of KOA patients presenting with and without disease progression in terms of structural cartilage degeneration will allow further refining of the threshold values of these functional biomarkers allowing identification of patients at risk of tissue degeneration. Furthermore, it should be noted that the relevance of the defined thresholds of the mechanoregulatory factors (MSS and SFD) for in vivo conditions is not clear, given that the current thresholds were selected based on outcomes from in vitro experiments (Eskelinen et al., 2019; Elahi et al., 2021). Further basic research should warrant the development of new technology to determine these thresholds in vivo and also to allow further finetuning, e.g., based on disease-specific mechanisms.

The focus of this study was to compare the effects of subject-specific loading and boundary conditions on the predicted mechanical environment of cartilage and meniscus in a healthy and early-stage KOA subject. In future studies, a sensitivity analysis can be conducted to specifically evaluate the influence of various subject-specific loading and boundary conditions or a combination of these conditions, which serve as an input to the FE model (such as TF anterior-posterior translation, TF medial-lateral translation, TF internal-external rotation, knee flexion angle, varus-valgus moment, and superior-inferior knee joint contact force), on the resulting tissue strains.

 5 Conclusion

Patient-specific gait data at baseline were processed through a multi-scale in silico modeling approach to analyze key mechanoregulatory factors for cartilage degeneration. These mechanoregulatory factors were able to distinguish between a healthy and progressive early-stage KOA subject and even more, able to identify the tibial compartment at risk of progression. Furthermore, evaluated tissue responses in the menisci can be related to specific OA loading conditions. Therefore, the mechanoregulatory factors analyzed could potentially serve as functional biomarkers allowing discrimination of patients at risk of (fast) disease progression. As such, this study presents a first step towards in silico-informed screening in the context of musculoskeletal disorders. Following further validation in larger cohorts, this could not only aid in earlier diagnosis, personalized treatment provision and prediction of progression as well as monitoring for patients with KOA, but may also provide novel fundamental insights into the intricate interplay between mechanical tissue loading during locomotion, the local cartilage micro-environment in vivo and consequent cartilage tissue degeneration.

 Data availability statement

The data analyzed in this study is subject to the following licenses/restrictions: The data was originally collected for other studies. Requests to access these datasets should be directed to IJ, ilse.jonkers@kuleuven.be.

 Ethics statement

The studies involving human participants were reviewed and approved by Approved by the local ethical committee at UZ LEUVEN (S50534). The patients/participants provided their written informed consent to participate in this study.

 Author contributions

IM and BK were responsible for data processing. IJ, SE, AE, and IM contributed to the study design and workflow. SV was responsible for longitudinal follow-up of gait and OA level of subjects. All authors contributed to the article and approved the submitted version.

 Funding

This work was supported by KU Leuven/FWO Happy Joints project (C14/18/077 and G045320N) and by the FWO, fundament fellowship (11M9422N). It is also part of the EOS excellence program Joint-Against OA (G0F8218N).

 Conflict of interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

 Publisher’s note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

 Supplementary material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fbioe.2023.1214693/full#supplementary-material

References Adouni M. Shirazi-Adl A. (2014). Evaluation of knee joint muscle forces and tissue stresses-strains during gait in severe OA versus normal subjects. J. Orthop. Res. 32 (1), 69–78. 10.1002/jor.22472 Al-Zahrani K. S. Bakheit A. M. O. (2002). A study of the gait characteristics of patients with chronic osteoarthritis of the knee. Disabil. Rehabilitation 24 (5), 275–280. 10.1080/09638280110087098 Andriacchi T. P. Mündermann A. (2006). The role of ambulatory mechanics in the initiation and progression of knee osteoarthritis. Curr. Opin. Rheumatology 18 (5), 514–518. 10.1097/01.bor.0000240365.16842.4e Astephen J. L. Deluzio K. J. Caldwell G. E. Dunbar M. J. (2008). Biomechanical changes at the hip, knee, and ankle joints during gait are associated with knee osteoarthritis severity. J. Orthop. Res. 26 (3), 332–341. 10.1002/jor.20496 Baert I. A. C. Jonkers I. Staes F. Luyten F. P. Truijen S. Verschueren S. M. (2013). Gait characteristics and lower limb muscle strength in women with early and established knee osteoarthritis. Clin. Biomech. 28 (1), 40–47. 10.1016/j.clinbiomech.2012.10.007 Baliunas A. J. Hurwitz D. Ryals A. Karrar A. Case J. Block J. (2002). Increased knee joint loads during walking are present in subjects with knee osteoarthritis. Osteoarthr. Cartil. 10 (7), 573–579. 10.1053/JOCA.2002.0797 Benninghoff A. (1925). Form und Bau der Gelenkknorpel in ihren Beziehungen zur Funktion - Erste Mitteilung: Die modellierenden und formerhaltenden Faktoren des Knorpelreliefs. Z. für Anat. Entwicklungsgeschichte 76 (1–3), 43–63. 10.1007/BF02134417 Blankevoort L. Huiskes R. (1991). Ligament-bone interaction in a three-dimensional model of the knee. J. Biomechanical Eng. 113 (3), 263–269. 10.1115/1.2894883 Bolcos P. O. Mononen M. E. Mohammadi A. Ebrahimi M. Tanaka M. S. Samaan M. A. (2018). Comparison between kinetic and kinetic-kinematic driven knee joint finite element models. Sci. Rep. 8 (1), 17351. 10.1038/s41598-018-35628-5 Cappozzo A. Catani F. Della Croce U. Leardini A. (1995). Position and orientation in space of bones during movement: Anatomical frame definition and determination. Clin. Biomech. 10 (4), 171–178. 10.1016/0268-0033(95)91394-T Choi Y. R. Kim J. H. Chung J. H. Lee D. H. Ryu K. J. Ha D. H. (2014). The association between meniscal subluxation and cartilage degeneration. Eur. J. Orthop. Surg. Traumatology 24 (1), 79–84. 10.1007/s00590-012-1144-3 Currey J. (2004). Incompatible mechanical properties in compact bone. J. Theor. Biol. 231 (4), 569–580. 10.1016/j.jtbi.2004.07.013 Dabiri Y. Li L. P. (2013). Influences of the depth-dependent material inhomogeneity of articular cartilage on the fluid pressurization in the human knee. Med. Eng. Phys. 35 (11), 1591–1598. 10.1016/J.MEDENGPHY.2013.05.005 Delp S. L. Anderson F. C. Arnold A. S. Loan P. Habib A. John C. T. (2007). OpenSim: Open-source software to create and analyze dynamic simulations of movement. IEEE Trans. Biomed. Eng. 54 (11), 1940–1950. 10.1109/TBME.2007.901024 Elahi S. A. Castro-Viñuelas R. Tanska P. Korhonen R. K. Lories R. Famaey N. (2023). Contribution of collagen degradation and proteoglycan depletion to cartilage degeneration in primary and secondary osteoarthritis: An in silico study. Osteoarthr. Cartil. 31, 741. 10.1016/j.joca.2023.01.004 Elahi S. A. Tanska P. Korhonen R. K. Lories R. Famaey N. Jonkers I. (2021). An in silico framework of cartilage degeneration that integrates fibril reorientation and degradation along with altered hydration and fixed charge density loss. Front. Bioeng. Biotechnol. 9, 680257. 10.3389/fbioe.2021.680257 Englund M. (2009). The role of the meniscus in osteoarthritis genesis. America: Medical Clinics of North America, 37–43. 10.1016/j.mcna.2008.08.005 Erdemir A. Besier T. F. Halloran J. P. Imhauser C. W. Laz P. J. Morrison T. M. (2019). Deciphering the “art” in modeling and simulation of the knee joint: Overall strategy. J. Biomechanical Eng. 141 (7), 0710021–07100210. 10.1115/1.4043346 Eskelinen A. S. A. Mononen M. E. Venäläinen M. S. Korhonen R. K. Tanska P. (2019). Maximum shear strain-based algorithm can predict proteoglycan loss in damaged articular cartilage. Biomechanics Model. Mechanobiol. 18 (3), 753–778. 10.1007/s10237-018-01113-1 Esrafilian A. (2021). Multiscale computational modeling of the knee joint: Development, validation, and application of musculoskeletal finite element modeling methods’, Multiscale computational modeling of the knee joint: Development, validation, and application of musculoskeletal finite element modeling methods , 428. China: Publications of the University of Eastern Finland. Dissertations in Forestry and Natural Sciences. Esrafilian A. Stenroth L. Mononen M. E. Tanska P. Avela J. Korhonen R. K. (2020). EMG-assisted muscle force driven finite element model of the knee joint with fibril-reinforced poroelastic cartilages and menisci. Sci. Rep. 10 (1), 3026–3116. 10.1038/s41598-020-59602-2 Esrafilian A. Stenroth L. Mononen M. E. Tanska P. Van Rossom S. Lloyd D. G. (2021). 12 degrees of freedom muscle force driven fibril-reinforced poroviscoelastic finite element model of the knee joint. IEEE Trans. Neural Syst. Rehabilitation Eng. 29, 123–133. 10.1109/TNSRE.2020.3037411 Esrafilian A. Stenroth L. Mononen M. E. Vartiainen P. Tanska P. Karjalainen P. A. (2022). An EMG-assisted muscle-force driven finite element analysis pipeline to investigate joint- and tissue-level mechanical responses in functional activities: Towards a rapid assessment toolbox. IEEE Trans. Biomed. Eng. 69 (9), 2860–2871. 10.1109/TBME.2022.3156018 Goodwin D. W. Wadghiri Y. Z. Zhu H. Vinton C. J. Smith E. D. Dunn J. F. (2004). Macroscopic structure of articular cartilage of the tibial plateau: Influence of a characteristic matrix architecture on MRI appearance. Am. J. Roentgenol. 182 (2), 311–318. 10.2214/ajr.182.2.1820311 Halonen K. S. Dzialo C. M. Mannisi M. Venäläinen M. S. de Zee M. Andersen M. S. (2017). Workflow assessing the effect of gait alterations on stresses in the medial tibial cartilage - combined musculoskeletal modelling and finite element analysis. Sci. Rep. 7 (1), 17396–17414. 10.1038/s41598-017-17228-x Halonen K. S. Mononen M. E. Jurvelin J. S. Töyräs J. Kłodowski A. Kulmala J. P. (2016). Importance of patella, quadriceps forces, and depthwise cartilage structure on knee joint motion and cartilage response during gait. J. Biomechanical Eng. 138 (7). 10.1115/1.4033516 Halonen K. S. Mononen M. E. Jurvelin J. S. Töyräs J. Korhonen R. K. (2013). Importance of depth-wise distribution of collagen and proteoglycans in articular cartilage-A 3D finite element study of stresses and strains in human knee joint. J. Biomechanics 46 (6), 1184–1192. 10.1016/j.jbiomech.2012.12.025 Halonen K. S. Mononen M. Jurvelin J. Töyräs J. Salo J. Korhonen R. (2014). Deformation of articular cartilage during static loading of a knee joint - experimental and finite element analysis. J. Biomechanics 47 (10), 2467–2474. 10.1016/j.jbiomech.2014.04.013 Heck T. A. M. Wilson W. Foolen J. Cilingir A. Ito K. van Donkelaar C. (2015). A tissue adaptation model based on strain-dependent collagen degradation and contact-guided cell traction. J. Biomechanics 48 (5), 823–831. 10.1016/j.jbiomech.2014.12.023 Heinegård D. Saxne T. (2011). The role of the cartilage matrix in osteoarthritis. Nat. Rev. Rheumatol. 7, 50–56. 10.1038/nrrheum.2010.198 Hurwitz D. E. Ryals A. R. Block J. A. Sharma L. Schnitzer T. J. Andriacchi T. P. (2000). Knee pain and joint loading in subjects with osteoarthritis of the knee. J. Orthop. Res. 18 (4), 572–579. 10.1002/JOR.1100180409 Julkunen P. Kiviranta P. Wilson W. Jurvelin J. S. Korhonen R. K. (2007). Characterization of articular cartilage by combining microscopic analysis with a fibril-reinforced finite-element model. J. Biomechanics 40 (8), 1862–1870. 10.1016/j.jbiomech.2006.07.026 Katz J. N. Arant K. R. Loeser R. F. (2021). Diagnosis and treatment of hip and knee osteoarthritis: A review. JAMA - J. Am. Med. Assoc. NIH Public Access 325, 568–578. 10.1001/jama.2020.22171 Kaufman K. R. (2001). “Gait characteristics of patients with knee osteoarthritis,” in Journal of Biomechanics (Germany: Elsevier), 907–915. 10.1016/S0021-9290(01)00036-7 Kellgren J. H. Lawrence J. S. (1957). Radiological assessment of osteo-arthrosis. Ann. rheumatic Dis. 16 (4), 494–502. 10.1136/ard.16.4.494 Lenhart R. L. Kaiser J. Smith C. R. Thelen D. G. (2015). Prediction and validation of load-dependent behavior of the tibiofemoral and patellofemoral joints during movement. Ann. Biomed. Eng. 43 (11), 2675–2685. 10.1007/s10439-015-1326-3 Lloyd D. G. Besier T. F. (2003). An EMG-driven musculoskeletal model to estimate muscle forces and knee joint moments in vivo . J. Biomechanics 36 (6), 765–776. 10.1016/S0021-9290(03)00010-1 Lu T. W. O’Connor J. J. (1999). Bone position estimation from skin marker co-ordinates using global optimisation with joint constraints. J. Biomechanics 32 (2), 129–134. 10.1016/S0021-9290(98)00158-4 Luyten F. P. Denti M. Filardo G. Kon E. Engebretsen L. (2012). Definition and classification of early osteoarthritis of the knee. Knee Surg. Sports Traumatol. Arthrosc. 20 (3), 401–406. 10.1007/s00167-011-1743-2 Mahmoudian A. van Dieen J. H. Bruijn S. M. Baert I. A. Faber G. S. Luyten F. P. (2016). Varus thrust in women with early medial knee osteoarthritis and its relation with the external knee adduction moment. Clin. Biomech. 39, 109–114. 10.1016/j.clinbiomech.2016.10.006 Makris E. A. Hadidi P. Athanasiou K. A. (2011). The knee meniscus: Structure-function, pathophysiology, current repair techniques, and prospects for regeneration. Biomaterials 32 (30), 7411–7431. 10.1016/j.biomaterials.2011.06.037 Martel-Pelletier J. Barr A. J. Cicuttini F. M. Conaghan P. G. Cooper C. Goldring M. B. (2016). Osteoarthritis. Nat. Rev. Dis. Prim. 2 (1), 16072–16118. 10.1038/nrdp.2016.72 Meireles S. Wesseling M. Smith C. R. Thelen D. G. Verschueren S. Jonkers I. (2017). Medial knee loading is altered in subjects with early osteoarthritis during gait but not during step-up-and-over task. PLoS ONE 12 (11), e0187583. 10.1371/journal.pone.0187583 Mononen M. E. Tanska P. Isaksson H. Korhonen R. K. (2016). A novel method to simulate the progression of collagen degeneration of cartilage in the knee: Data from the osteoarthritis initiative. Sci. Rep. 6 (1), 21415–21514. 10.1038/srep21415 Mononen M. E. Mikkola M. Julkunen P. Ojala R. Nieminen M. Jurvelin J. (2012). Effect of superficial collagen patterns and fibrillation of femoral articular cartilage on knee joint mechanics—a 3D finite element analysis. J. Biomechanics 45 (3), 579–587. 10.1016/J.JBIOMECH.2011.11.003 Mononen M. E. Tanska P. Isaksson H. Korhonen R. K. (2018). New algorithm for simulation of proteoglycan loss and collagen degeneration in the knee joint: Data from the osteoarthritis initiative. J. Orthop. Res. 36 (6), 1673–1683. 10.1002/jor.23811 Mukherjee S. Nazemi M. Jonkers I. Geris L. (2020). Use of computational modeling to study joint degeneration: A review. Front. Bioeng. Biotechnol. 8, 93. 10.3389/fbioe.2020.00093 Netravali N. A. Giori N. J. Andriacchi T. P. (2010). Partial medial meniscectomy and rotational differences at the knee during walking. J. Biomechanics 43 (15), 2948–2953. 10.1016/j.jbiomech.2010.07.013 Orozco G. A. Bolcos P. Mohammadi A. Tanaka M. S. Yang M. Link T. M. (2021). Prediction of local fixed charge density loss in cartilage following ACL injury and reconstruction: A computational proof‐of‐concept study with MRI follow‐up. J. Orthop. Res. 39 (5), 1064–1081. 10.1002/jor.24797 Orozco G. A. Tanska P. Florea C. Grodzinsky A. J. Korhonen R. K. (2018) ‘A novel mechanobiological model can predict how physiologically relevant dynamic loading causes proteoglycan loss in mechanically injured articular cartilage’, Sci. Rep., 8(1), pp. 15599–15616. 10.1038/s41598-018-33759-3 Párraga Quiroga J. M. Wilson W. Ito K. van Donkelaar C. C. (2017). The effect of loading rate on the development of early damage in articular cartilage. Biomechanics Model. Mechanobiol. 16 (1), 263–273. 10.1007/s10237-016-0815-0 Sinusas K. (2012). Osteoarthritis: Diagnosis and treatment. Am. Fam. Physician 85 (1), 49–56. Smith C. R. Lenhart R. Kaiser J. Vignos M. Thelen D. (2014). Influence of ligament properties on tibiofemoral mechanics in walking. J. Knee Surg. 29 (2), 099–106. 10.1055/s-0035-1558858 Smith C. R. Won Choi K. Negrut D. Thelen D. G. (2018). Efficient computation of cartilage contact pressures within dynamic simulations of movement. Comput. Methods Biomechanics Biomed. Eng. Imaging Vis. 6 (5), 491–498. 10.1080/21681163.2016.1172346 Sophia Fox A. J. Bedi A. Rodeo S. A. (2009). The basic science of articular cartilage: Structure, composition, and function. Sports Health 1 (6), 461–468. 10.1177/1941738109350438 Speirs A. D. Beaulé P. E. Ferguson S. J. Frei H. (2014). Stress distribution and consolidation in cartilage constituents is influenced by cyclic loading and osteoarthritic degeneration. J. Biomechanics 47 (10), 2348–2353. 10.1016/j.jbiomech.2014.04.031 Sun H. B. (2010). Mechanical loading, cartilage degradation, and arthritis. Ann. N. Y. Acad. Sci. 1211 (1), 37–50. 10.1111/j.1749-6632.2010.05808.x Syed N. Makris U. E. (2020). “Osteoarthritis,” in Pain management for clinicians (Cham: Springer International Publishing), 175–180. 10.1007/978-3-030-39982-5_7 Tanska P. Julkunen P. Korhonen R. K. (2018). A computational algorithm to simulate disorganization of collagen network in injured articular cartilage. Biomechanics Model. Mechanobiol. 17 (3), 689–699. 10.1007/s10237-017-0986-3 Verdonk R. Madry H. Shabshin N. Dirisamer F. Peretti G. M. Pujol N. (2016). The role of meniscal tissue in joint protection in early osteoarthritis. Knee Surg. Sports Traumatol. Arthrosc. 24 (6), 1763–1774. 10.1007/s00167-016-4069-2 Wilson W. van Donkelaar C. van Rietbergen B. Huiskes R. (2005). A fibril-reinforced poroviscoelastic swelling model for articular cartilage. J. Biomechanics 38 (6), 1195–1204. 10.1016/j.jbiomech.2004.07.003 Wilson W. van Donkelaar C. van Rietbergen B. Ito K. Huiskes R. (2004). Stresses in the local collagen network of articular cartilage: A poroviscoelastic fibril-reinforced finite element study. J. Biomechanics 37 (3), 357–366. 10.1016/S0021-9290(03)00267-7