Mechanobiological model for simulation of injured cartilage degradation via pro-inflammatory cytokines and mechanical stimulus

Abstract

Post-traumatic osteoarthritis (PTOA) is associated with cartilage degradation, ultimately leading to disability and decrease of quality of life. Two key mechanisms have been suggested to occur in PTOA: tissue inflammation and abnormal biomechanical loading. Both mechanisms have been suggested to result in loss of cartilage proteoglycans, the source of tissue fixed charge density (FCD). In order to predict the simultaneous effect of these degrading mechanisms on FCD content, a computational model has been developed. We simulated spatial and temporal changes of FCD content in injured cartilage using a novel finite element model that incorporates (1) diffusion of the pro-inflammatory cytokine interleukin-1 into tissue, and (2) the effect of excessive levels of shear strain near chondral defects during physiologically relevant loading. Cytokine-induced biochemical cartilage explant degradation occurs near the sides, top, and lesion, consistent with the literature. In turn, biomechanically-driven FCD loss is predicted near the lesion, in accordance with experimental findings: regions near lesions showed significantly more FCD depletion compared to regions away from lesions (p<0.01). Combined biochemical and biomechanical degradation is found near the free surfaces and especially near the lesion, and the corresponding bulk FCD loss agrees with experiments. We suggest that the presence of lesions plays a role in cytokine diffusion-driven degradation, and also predisposes cartilage for further biomechanical degradation. Models considering both these cartilage degradation pathways concomitantly are promising in silico tools for predicting disease progression, recognizing lesions at high risk, simulating treatments, and ultimately optimizing treatments to postpone the development of PTOA.

About the Authors:

Atte S. A. Eskelinen

Roles Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Visualization, Writing – original draft

* E-mail: [email protected]

Affiliation: Department of Applied Physics, University of Eastern Finland, Finland

ORCID logo http://orcid.org/0000-0003-1516-3083

Petri Tanska

Roles Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Software, Supervision, Writing – review & editing

Affiliation: Department of Applied Physics, University of Eastern Finland, Finland

ORCID logo http://orcid.org/0000-0002-9684-6902

Cristina Florea

Roles Conceptualization, Data curation, Formal analysis, Investigation, Writing – review & editing

Affiliations Department of Applied Physics, University of Eastern Finland, Finland, Departments of Biological Engineering, Electrical Engineering and Computer Science and Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, United States of America

ORCID logo http://orcid.org/0000-0001-7719-2600

Gustavo A. Orozco

Roles Conceptualization, Data curation, Investigation, Methodology, Software, Writing – review & editing

Affiliation: Department of Applied Physics, University of Eastern Finland, Finland

Petro Julkunen

Roles Conceptualization, Funding acquisition, Methodology, Supervision, Writing – review & editing

Affiliations Department of Applied Physics, University of Eastern Finland, Finland, Department of Clinical Neurophysiology, Kuopio University Hospital, Finland

ORCID logo http://orcid.org/0000-0002-0461-1058

Alan J. Grodzinsky

Roles Conceptualization, Methodology, Resources, Supervision, Writing – review & editing

Affiliation: Departments of Biological Engineering, Electrical Engineering and Computer Science and Mechanical Engineering, Massachusetts Institute of Technology, Cambridge, United States of America

Rami K. Korhonen

Roles Conceptualization, Funding acquisition, Investigation, Methodology, Project administration, Resources, Supervision, Writing – review & editing

Affiliation: Department of Applied Physics, University of Eastern Finland, Finland

ORCID logo http://orcid.org/0000-0002-3486-7855

Introduction

Acute joint insult can debilitate the functioning of articular cartilage and lead to pain, joint stiffness and disability in patients. Moreover, trauma such as anterior cruciate ligament (ACL) rupture [1,2] may be associated with cartilage injury and increased susceptibility to cartilage degeneration which can culminate in post-traumatic osteoarthritis (PTOA) [3–5]. Two major mechanisms have been suggested to play a role in PTOA progression. The first mechanism involves joint inflammation [3,6,7] and the resulting diffusion of pro-inflammatory cytokines from the synovial fluid into cartilage [8–11], leading to proteolysis of cartilage matrix. The second mechanism involves biomechanical factors [12] including induction of chondral defects and aberrant loading characteristics in the knee joint milieu [13–15], leading to elevated shear strains near chondral lesions [16–18]. While these mechanisms have been recognized for years, current healthcare approaches lack effective tools to identify patients with increased risk of PTOA development. Prevention or even postponing the onset of the disease and surgical interventions would be desirable. Thus, we urgently need novel tools to predict disease progression.

Computational modeling is a critically important and cost-efficient approach which can serve as a predictive tool of possible PTOA progression. Finite element (FE) modeling has been successfully used to assess cartilage damage at both tissue [10,17,19] and joint levels [20,21]. Furthermore, subject-specific numerical models can be helpful in optimization of personalized intervention strategies, rehabilitation procedures, and surgeries. They could also be salutary in development of immunomodulating therapies [22,23] and disease-modifying osteoarthritis drugs [9,24,25]. When validated with experimental data, mechanobiological models can provide key quantitative and qualitative insights into biochemical and biomechanical adaptive processes in cartilage. However, computational models combining both inflammatory and loading mechanisms are scarce [26].

Inflammation plays a multifaceted role post-injury. The innate inflammatory response is an inherent part of host defense mechanisms aimed at maintaining tissue integrity and homeostasis. When tissue homeostasis is disturbed via inflammation, synovial macrophages [27,28] and elevated concentrations of inflammatory messenger molecules (such as alarmins, chemokines, and cytokines) orchestrate signaling cascades that promote degradation of the extracellular matrix (ECM) of cartilage [3,29,30]. Cytokines originate primarily from cells in the synovial lining; subsequent increase in cytokine signaling can lead to recruitment of additional immune cells to the inflamed area, and the newly recruited cells produce even more cytokines. These small proteins diffuse into cartilage [8] and bind to cell receptors on chondrocyte (cartilage cell) surface, thus regulating the balance of anabolic and catabolic gene expression in cartilage [31]. Some cytokines in cartilage are pro-anabolic (anti-inflammatory), such as transforming growth factor beta (TGF-β), interleukin (IL)-1 receptor antagonist, IL-4, IL-10, and erythropoietin [22,28,32]. Other cytokines are pro-catabolic (pro-inflammatory), such as tumor necrosis factor alpha (TNFα), IL-1, IL-6, IL-8, IL-12, and IL-15 [3,8,23,33,34]. Previous experiments [4,35–37] suggest that chondrocytes affected by pro-inflammatory cytokines upregulate production of proteolytic aggrecanases ADAMTS-4,5 (a-disintegrin and metalloproteinase with thrombospondin motifs-4,5) and matrix metalloproteinases including the collagenases MMP-1,13. These proteases cleave matrix proteoglycans (PGs) and collagen fibrils, respectively. However, tissue inhibitors of metalloproteinases (TIMPs) can slow these catabolic processes by ADAMTS (specifically, TIMP-3 inhibits ADAMTS-4,5 [38]) and MMPs [39]. To predict the experimentally observed cartilage matrix loss [40], Kar and coworkers [9,10] developed a numerical model including the diffusion of IL-1 into cartilage and the subsequent degeneration of tissue aggrecan (the main PG in articular cartilage) and collagen. Other in silico studies [22,26,41] have also sought to cast new light on inflammation and complex biochemically-induced degradation mechanisms of cartilage.

Mechanical loading is also an essential regulator of cartilage homeostasis [42]. Previous experimental studies have proposed that the harsh loading environment of the knee joint [43], unfavorable loading frequency [44], and presence of chondral injuries [45] can predispose cartilage to biomechanical degradation. In efforts to replicate experimental findings, computational models have focused on the effects of mechanical loading-modulated tissue stresses (such as maximum principal stress) [20,46], strains [17,19,47], fiber stretches [41], volumetric deformations [48], interstitial fluid velocities [49], and altered tissue swelling properties [50] on induction of ECM degeneration. These models have included collagen fibril damage [41,46], PG or related fixed charge density (FCD) loss [17,49], tissue softening [47,51], and increase of tissue permeability [51].

Taking a step towards tools estimating PTOA progression, in this study we provide novel mechanobiological insights into biochemically and biomechanically-driven bovine cartilage degradation shortly after traumatic injury. Building on the previous research by Kar et al. [9,10] and Orozco et al. [49], we present a mechanobiological fibril-reinforced porohyperelastic swelling model of cartilage degeneration which considers biochemical (diffusion of pro-inflammatory IL-1 into tissue and subsequent release of ADAMTS) and biomechanical (elevated levels of the maximum shear strain during physiologically relevant dynamic loading) degradation mechanisms separately and simultaneously. We simulate spatial and temporal decreases in FCD and compare the numerical predictions with previous experimental results on cartilage plugs harvested from 1-2-week-old calves [35,36,40,49,52]. To the best of our knowledge, this is the first computational study involving the known extent of chondral lesions and including simultaneously both biochemical and biomechanical degradation pathways. We concentrate here on IL-1 due to its known role in cartilage inflammation [10,28] and the comparative lack of quantitative information on the dynamic temporal patterns of simultaneously acting inflammatory cytokines (such as IL-6, TNFα) and subsequent cartilage adaptation. Furthermore, we here include only FCD loss and not yet collagen degradation [10,41] or reorientation [53,54], since the reported experiments [49] were performed on a relatively short timescale (12 days); during this short time, substantial collagen degradation has not yet occurred and thus collagen content was assumed to remain unaltered as suggested in earlier animal model experiments [55]. In addition, as one of the early changes in PTOA, FCD loss was presumed to occur prior to collagen damage [10,40,56], consistent with observations by Li et al. [40] and Kar et al. [10] that collagen loss was initiated only after 12 days of IL-1 culture.

See Table 1 for a list of abbreviations.

[Figure omitted. See PDF.]

Table 1. List of abbreviations.

https://doi.org/10.1371/journal.pcbi.1007998.t001

Materials and methods

Experiments

Previous experimental research has assessed cytokine-mediated and biomechanically-driven cartilage degeneration [36,40,49,52,57,58]. These studies involved evaluation of (1) loss of sulphated glycosaminoglycans (GAGs) to culture medium via the dimethylmethylene blue (DMMB) assay, or (2) loss of FCD via decrease in optical density obtained with digital densitometry (DD) of Safranin-O stained sections [59,60]. Both methods are indicative of PG matrix damage. In the current study, we quantified average FCD loss near and further away from cartilage lesions based on the data of Orozco et al. [49] to (1) characterize how the presence of lesions affects early matrix damage and PG release, and (2) to compare biomechanically-induced FCD loss with our computational estimates. Due to lack of data in the literature, we did not focus on localized FCD loss in injured samples cultured with cytokines. See S1 Supplementary Material (Subsections S1.1–1.4) for additional details on the previously reported experimental data used for the present modeling analyses.

Computational mechanobiological models

We developed mechanobiological models to predict the loss of cartilage FCD content over time caused by inflammatory and shear strain mechanisms. In the biochemical model (Fig 1A), FCD loss was associated with diffusion of the cytokine IL-1 into cartilage and subsequent increased levels of ADAMTS proteases. In the biomechanical model (Fig 1B), the decrease in the FCD content was driven by excessive levels of the maximum shear strain. Lastly, we implemented both of these mechanisms simultaneously in a combined biochemical and biomechanical model (Fig 1C). The finite element (FE) model simulations were conducted using COMSOL Multiphysics (version 5.3a, Burlington, MA, USA) for the biochemical model, ABAQUS (2018, Dassault Systèmes, Providence, RI, USA) for the biomechanical model, and both software packages simultaneously for the combined model. Furthermore, a custom-made interface was developed in MATLAB (R2017b, The MathWorks, Inc., Natick, MA, USA) for handling the flow of PG matrix degradation data between COMSOL and ABAQUS. The cartilage disk thickness and radius for the 2D FE model were h = 1 mm and r = 1.5 mm, respectively. The cartilage lesion geometry (depth 144 μm, width 316 μm) was adopted from the study by Orozco et al. [49] in which the FE mesh generated in ABAQUS was based on segmentation of a histological slice from a representative injured cartilage disk. The output of all the models was the regional average FCD loss (decrease in FCD levels with respect to initial levels) over time for the regions (1) within 0.1 mm (±10%) from the lesion, and (2) 0.15 mm deep (±10%, from the surface) and 0.45 mm wide (±10%) away from the lesion at a location in the middle between the lesion and disk edges, as performed in the analysis of experimental biomechanical FCD depletion (see S1 Supplementary Material Subsection S1.2). Moreover, we created an intact reference model to analyze the FCD loss from the middle of the disk in a same-sized region as away from the lesion in the injury model. In all the models, the possible minimum FCD concentration was set to 10% of the initial minimum FCD value (found from the superficial zone) to avoid convergence issues. The most essential model inputs, outputs, and equations defining the computational FCD loss are recapitulated in Fig 2. For more details, readers are referred to Kar et al. [10], Orozco et al. [49], and the S1 Supplementary Material (Subsections S1.5–1.9, S1 Table).

[Figure omitted. See PDF.]

Fig 1. Cartilage degeneration algorithm.

Outline of the proposed proteoglycan degeneration algorithm based on the loss of fixed charge density (FCD) via A) diffusion of pro-inflammatory cytokine interleukin-1 into the tissue and subsequent biochemical degradation of cartilage matrix, B) excessive levels of the maximum shear strain especially near lesions and C) combination of these cartilage-degrading mechanisms. In the combined degradation model, both mechanisms are in effect simultaneously during each simulation iteration, followed by an update of non-fibrillar matrix contents for the next iteration. Experimental biomechanical degradation images modified with permission from Orozco et al. [49].

https://doi.org/10.1371/journal.pcbi.1007998.g001

[Figure omitted. See PDF.]

Fig 2. Workflow.

Workflow for creation of A) biochemical, B) biomechanical and C) combined biochemical and biomechanical degradation models. Lesion geometry was based on histological findings after injurious loading. All of the mechanobiological models estimated fixed charge density (FCD) loss near and away from the lesion over time. Histology image obtained with permission from Orozco et al. [49].

https://doi.org/10.1371/journal.pcbi.1007998.g002

Simulation of biochemical degradation

Biochemically-driven aggrecan loss was modeled via diffusion of IL-1 cytokines (1 ng/ml in the culture medium) into tissue and subsequent increase in protease biomolecule concentrations (time integration in COMSOL Multiphysics). In short, the model [10] considered the net degradative effect of ADAMTS-4,5 and TIMP-3 on aggrecan content. The model inputs include effective diffusion coefficients of IL-1 [10,61], intact aggrecan [62–64], and aggrecanases [65], aggrecan biosynthesis rate depending on the local IL-1 concentration [10,66–68], initial aggrecan distribution [69,70], catabolic rate for aggrecan cleavage by aggrecanases [71], and production/degradation of ADAMTS [10,72]. In the COMSOL model, the primary variable of interest was intact aggrecan concentration, the changes of which were interpreted as changes in FCD concentration in ABAQUS (see Subsection Materials and methods: Biochemical fixed charge density loss; ultimately, the biochemical model output was FCD loss). We chose to model biochemical degradation for 21 days because in the model by Kar et al. [10] the basal rate of aggrecan synthesis and the mass transfer coefficient of aggrecan were calibrated within this timeframe with respect to PG loss in young animals [66]. In COMSOL, we selected four-node linear quadrilateral plate elements (element type CQUAD4) with the same node locations as in the ABAQUS mesh. This mesh was imported into COMSOL in a NASTRAN-format file generated with a custom-made MATLAB script.

Cytokine diffusion model

The diffusion of IL-1 into porous biphasic cartilage and proteolytic effects of aggrecanases were modeled with parabolic reaction—diffusion partial differential equations [10,73](1)where Cj is the concentration/amount of constituent j, t is time, Dj is the effective diffusivity of chemical species j and Rj is the corresponding source/sink term, which describes the rate of generation/repair or degradation/consumption of individual species. The species j includes chondrocytes, IL-1, intact aggrecans, and aggrecanases. The source/sink terms Rj for proteolytic aggrecanases were designed in accordance with Michaelis—Menten kinetics. Experimentally observed time delay in secretion of ADAMTS after initiation of biochemical challenge [74] was included in the form of a stimulus equation [10](2)where S1 is a variable describing IL-1-mediated stimulus response of ADAMTS, C* is the concentration of IL-1-IL-1-receptor complexes and α1 is rate constant for stimulus. This stimulus variable S1 was used to define source terms Rj in Eq (1) for ADAMTS.

Initial and boundary conditions in biochemical model

As initial condition, the concentration of IL-1 and aggrecanases were set to 0 ng/ml, and aggrecanase-stimulus (S1 in Eq (2) at time t = 0 d) to 0 mol/m3 inside the tissue. The chondrocyte concentration [75] was uniform and constant throughout the simulation, and no apoptosis nor necrosis was modeled due to lack of cell death rate (dead cells/(s⋅m3)) data in the literature.

For the IL-1 concentration, Dirichlet boundary condition of 1 ng/ml was set on the free surfaces (top, lesion, and lateral) throughout the simulation. We assumed a constant cytokine level outside tissue since during the experiments the culture medium was changed every other day. For the same reason, concentration of aggrecanases was set to 0 ng/ml on the free and bottom surfaces. Mass transfer of intact aggrecan was allowed through the free surfaces. This was handled by Robin boundary conditions (combination of Dirichlet and Neumann boundary conditions) to account for diminished diffusion of chemical species through the cartilage-culture medium interface [76]. No mass transfer nor diffusion was allowed for any of the chemical species through the bottom surface.

Biochemical fixed charge density loss

In the COMSOL model, the primary degrading constituent is intact aggrecan, having concentration Cagg. Intact aggrecan concentration is affected by the effects of aggrecanases and by synthesis of new aggrecans by the chondrocytes. However, to compare the FCD concentration distribution (which is the primary variably of interest in ABAQUS model) between all the mechanobiological models, intact aggrecan concentration in COMSOL was converted to FCD concentration in ABAQUS. Therefore, relative changes in integration point-wise aggrecan concentrations during one simulation iteration i (50 iterations in total corresponded to 21 days of simulated time, so each iteration simulated a time period of 21/50 d ≈ 10 h; see more details from Subsection Materials and methods: Simulation of combined biochemical and biomechanical degradation) were implemented with MATLAB as coefficients of FCD depletion for each integration point(3)where FCDi is the integration point-wise FCD concentration in ABAQUS at iteration i, Cagg,i,end and Cagg,i,start are the integration point-wise aggrecan concentrations in COMSOL at the end and start of iteration i, and is the biochemical degradation rate. The mesh was the same in both COMSOL and ABAQUS. After each iteration, the decreased FCD in each integration point was imported to COMSOL as new integration point-wise aggrecan concentration for the following iteration Cagg,i+1,start(4)where for the biochemical model Cagg,i+1,start = Cagg,i,end. This was done with the COMSOL import tool, which used linear interpolation (see S1 Supplementary Material Subsection S1.6).

Simulation of biomechanical degradation

Biomechanical degradation of PG content was associated with increased levels of intra-tissue shear strain during loading, since shear strains have been suggested to increase near chondral lesions [16–18,77] and may induce cartilage matrix breakdown [19] and chondrocyte apoptosis [78–80]. Ultimately, this leads to loss of FCD content. The biomechanical model solved a highly non-linear problem consecutively (50 times, updating FCD content between iterations) with a transient SOILS analysis in ABAQUS to account for the viscoelastic behavior of the fibril-reinforced porohyperelastic material with Donnan osmotic swelling and chemical expansion (see S1 Supplementary Material Subsection S1.7 for further details on the material model). The structural inputs for the biomechanical model included depth-dependent fluid fraction [81], FCD [49], and collagen distributions [82], and collagen fibril orientation in immature tissue [49,83]. The output of the biomechanical model was FCD loss. Linear axisymmetric continuum elements with pore pressure (element type CPE4P) were implemented for the model in ABAQUS.

Boundary conditions in biomechanical model

Initially, the cartilage was allowed to reach mechanical equilibrium with external salt concentration via swelling [84]. Then the cartilage disk was subjected to the same dynamic loading protocol (modeled for two full loading-unloading cycles (two seconds), since our preliminary testing suggested that a longer simulation did not cause any substantial changes in strains) with an impermeable platen as in the experiments. During loading, the cartilage-platen contact was frictionless (contact type: surface-to-surface; contact behavior: penalty), the explant was allowed to bulge radially (bottom surface was fixed axially, one bottom corner node was also fixed horizontally to prevent the disk from rotating and translating), and fluid flow was allowed through the free surfaces.

Biomechanical fixed charge density loss

The FCD content was decreased in each integration point when the maximum shear strain levels [17,85](5)exceeded the degeneration threshold of = 50%, where εp,k are the principal strains of the Green—Lagrangian strain tensor. This threshold value was chosen to match experimental findings and predicted biomechanical FCD loss as done by Orozco et al. [49]. The integration point-wise maximum shear strain during loading (at the second peak of compressive 1-Hz load during two-second simulation) and the degeneration threshold defined together the biomechanical degeneration rate Dmech [17,51](6)

See also the Subsection Materials and methods: Simulation of combined biochemical and biomechanical degradation for justification of this equation. Thus, the integration point-wise FCD loss for simulation iteration i was defined as(7)

In total, the model was run for 50 iterations to reach near-equilibrium (see Subsection Results: Simulation of biomechanical degradation) in terms of the predicted FCD distribution. The loss of FCD affected the Donnan osmotic swelling pressure gradient Δπ, chemical expansion stress Tc, and ultimately the total stress tensor σtot (see S1 Supplementary Material Subsection S1.7 for definition of these terms). After each iteration, the decreased FCD was imported as a new FCD distribution for the following ABAQUS simulation iteration (cF,0 in Eq. (S15) in S1 Supplementary Material Subsection S1.7).

Simulation of combined biochemical and biomechanical degradation

In a third model, both biochemical and biomechanical degradative mechanisms were acting concomitantly. (See Subsections Materials and methods: Initial and boundary conditions in biochemical model and Materials and methods: Boundary conditions in biomechanical model above for initial and boundary conditions for both of these mechanisms, respectively.) The combined FCD loss was implemented similarly as in Eqs (3) and (7)(8)and the integration point-wise FCD loss in ABAQUS was used to solve the intact aggrecan distribution to be imported to COMSOL with Eq (4). Similar approach (multiplication of the biochemical and biomechanical degradation effects) has also been previously used [26].

Simulation time of 21 days in COMSOL was chosen to correspond to 50 iterations in ABAQUS. 21 days was used originally for the biochemical model calibration [10] and 50 iterations was found to result in near-equilibrium of biomechanically-driven FCD loss (negligible change in FCD loss with more iterations, see Subsection Results: Simulation of biomechanical degradation). To further justify our decision, the biomechanical degeneration rate (Eq 6) [17,51] combined with this time scale and number of iterations captures well the FCD loss near lesions in our experiments (see Subsection Results: Simulation of biomechanical degradation and S1 Supplementary Material Subsection S2.1) and the rapid early loading-associated GAG loss in mechanically loaded samples [83,86–88] (see Subsection Discussion: Biomechanical degradation).

Mesh density

The sensitivity of the predicted FCD loss on the FE mesh was tested with five different meshes and the algorithm including both biochemically and biomechanically-driven degradations. These meshes were increasingly dense near the lesion and the free edges. A mesh with 918 elements was chosen, since increasing the mesh density beyond that resulted in only minor changes in predicted FCD loss within 0.1 mm from the lesion and away from the lesion (S1 Supplementary Material Subsection S1.9, S1 Fig).

Results

Summary of cartilage degradation

Representative samples of experimental biochemically-driven degradation (Fig 3A) demonstrated progressive FCD loss starting from the free surfaces and extending over time to the deeper tissues. Mechanically injured and dynamically loaded samples showed pronounced FCD loss near defects (Fig 3B). Simulated biochemical degradation occurred near the top surface, edges and below the lesion (Fig 4). On the other hand, with biomechanical degradation algorithm we predicted matrix losses only near the lesion (Fig 5). In the model combining both mechanisms, matrix losses occurred near the free surfaces and especially in the vicinity of the lesion (Fig 6).

[Figure omitted. See PDF.]

Fig 3. Previous experimental findings.

A) In vitro experiments with exogenous pro-inflammatory cytokine interleukin-1 (IL-1) challenge (10 ng/ml in culture medium) show marked matrix degradation near the sample surface and edges [40]. B) In vitro experiments with injurious loading (50% at rate 100%/s) following a dynamic loading period (15% strain amplitude, haversine waveform, 1 Hz, 1h on and 5h off cycles) in unconfined compression show lower optical density (~lower fixed charge density, FCD) especially near lesions. The numbers [%] show the average localized FCD loss compared to the day when cartilage was mechanically injured. Biomechanical degradation images modified with permission from Orozco et al. [49].

https://doi.org/10.1371/journal.pcbi.1007998.g003

[Figure omitted. See PDF.]

Fig 4. Simulated biochemical cartilage degradation.

Predicted biochemical fixed charge density (FCD) losses in A) intact reference model and B) injury model. Biochemical degradation shown at time t = 4 d was simulated with 1 ng/ml of exogenous interleukin-1.

https://doi.org/10.1371/journal.pcbi.1007998.g004

[Figure omitted. See PDF.]

Fig 5. Simulated biomechanical cartilage degradation.

Predicted biomechanical fixed charge density (FCD) losses in A) intact reference model and B) injury model. Biomechanical degradation shown at time t = 12 d was simulated with the degeneration threshold of = 50% of the maximum shear strain. This threshold was chosen to match model predictions and experimental findings.

https://doi.org/10.1371/journal.pcbi.1007998.g005

[Figure omitted. See PDF.]

Fig 6. Simulated combined biochemical and biomechanical cartilage degradation.

Predicted combined biochemical and biomechanical fixed charge density (FCD) losses in A) intact reference model and B) injury model. Combined biochemical and biomechanical degradation shown at time t = 4 d was simulated with 1 ng/ml of exogenous interleukin-1 and with the degeneration threshold of = 50% of the maximum shear strain. This threshold was chosen to match model predictions and experimental findings.

https://doi.org/10.1371/journal.pcbi.1007998.g006

Simulation of biochemical degradation

Biochemical degradation occurred near the free surfaces in both intact reference (Fig 4A) and injury (Fig 4B) models. At timepoint t = 4 d, the predicted average FCD losses were 38% for the intact reference model, and 44% both near and away from the lesion for the injury model. Below the lesion, the ECM was degraded deeper in the tissue than closer to the sides. Near the edge of the lesion in the superficial zone, the maximum FCD losses were 58% (t = 4 d) and 97% (t = 12 d) (Figs 4B and 7A). Temporally, the average regional biochemical FCD loss was of sigmoidal shape; first occurring at a relatively slow pace, then attaining the highest rate of degradation at around day 4, and finally reaching the maximum loss set up by our algorithm (10% of the initial minimum FCD was the lowest allowable FCD content) at around day 8 (Fig 7A). The localized peak FCD loss profile with logarithmic-like shape showed a rapid degradation near the edge of the lesion (Fig 7A). However, the overall degradation of cartilage did not cease at day 8 as shown by a steady increase of the degraded area after day 8 (Fig 8A; cartilage was referred as “degraded” if element-wise FCD loss compared to initial FCD content was equal or more than 20%). The degraded area increased rapidly around day 2 (Fig 8A) starting from the free surfaces (see S1 Supplementary Material, Subsection S2.4, S5 Fig, S1 Animation).

[Figure omitted. See PDF.]

Fig 7. Temporal estimations of cartilage degradation.

Average and maximum fixed charge density (FCD) losses over time predicted by A) biochemical, B) biomechanical, and C) combined biochemical and biomechanical model. The dashed line in biomechanical model includes the base degeneration rate parameter, which was chosen according to the experimental FCD loss in freely-swollen cartilage disks (see S1 Supplementary Material, Subsection S1.8). The blue and black errorbars (with 95% confidence intervals) at days 7 and 12 represent the experimental average FCD loss near lesions after dynamic loading and in freely-swollen samples, respectively. The FCD concentration was not allowed to decrease below 10% of the initial minimum FCD concentration in the model due to computational stability.

https://doi.org/10.1371/journal.pcbi.1007998.g007

[Figure omitted. See PDF.]

Fig 8. Simulation of degraded cartilage area.

Temporal estimations of degraded cartilage area with A) biochemical and combined biochemical and biomechanical degradation and B) biomechanical degradation. In these simulations, cartilage was considered”degraded” when FCD loss in elements was equal or greater than 20%. See S1 Supplementary Material (Subsection S2.4, S5 Fig, S1–S4 Animation) for animations showing the propagation of the degraded area over time.

https://doi.org/10.1371/journal.pcbi.1007998.g008

Simulation of biomechanical degradation

The intact reference model showed no FCD loss anywhere (Fig 5A), whereas the injury model showed an average FCD loss of 35% near the lesion at timepoint t = 12 d (Fig 5B). The localized peak FCD losses followed logarithmic-like shape and were 89% (t = 4 d) and 97% (t = 12 d) right below the lesion surface (Figs 5B and 7B). The average regional FCD loss occurred at a fast rate early-on (days 0–5), then leveling off to a slower rate over time (Fig 7B). At day 21, the predicted FCD distribution had reached near-equilibrium (Fig 7B, solid blue line; the FCD loss changed less than 0.15% per iteration after iteration 40). An addition of the base degeneration rate term (indicating spontaneous GAG loss from unloaded samples via passive diffusion, see S1 Supplementary Material Subsection S1.8) increased the predicted average FCD loss (Fig 7B). Then, the model predictions suggested a 6-7-fold (at day 2) and a 3-fold (at day 12) increase in FCD loss in areas near the lesion compared to areas away from the lesion (where base degeneration accounts for all of the FCD loss; Fig 7B). The degraded area increased early (days 0–2) at a fast pace and slowly throughout the rest of the simulation (Fig 8B and S1 Supplementary Material Subsection S2.4, S5 Fig, S2 Animation and S3 Animation).

Simulation of combined biochemical and biomechanical degradation

The intact reference model predicted only biochemical degradation (38% FCD loss at t = 4 d, Fig 6A). However, combined biochemical and biomechanical degradation occurred near the free surfaces and especially in the vicinity of the lesion in the injury model (Fig 6B). The predicted average FCD losses were 55% and 44% near and away from the lesion, respectively. The localized maximum FCD losses were 81% (t = 4 d) and 98% (t = 12 d) below the lesion surface (Figs 6B and 7C). The FCD loss reached deeper layers of cartilage below the lesion compared to the tissues near lateral edges. During the early days of simulation (days 0–6), the predicted FCD loss near the defect was higher than away from the defect, where the FCD loss profile exhibited sigmoidal behavior (Fig 7C). Furthermore, the rate of FCD loss was almost constant near the lesion during the early days and initially higher near the lesion than away from it (Fig 7C). Below the lesion, the logarithmic-shaped maximum FCD loss profile predicted fast degradation (Fig 7C). At around day 8, the average FCD loss near and away the lesion reached equilibrium, but the overall cartilage degeneration continued (Fig 8A). Animations in the S1 Supplementary Material (Subsection S2.4, S5 Fig, S4 Animation) and Fig 8A show the rapid increase of degraded areas near the lesion and subsequently near the free surfaces.

Discussion

Summary

In this study, we have developed a novel mechanobiological FE model capable to predict cartilage degradation via both biochemically and biomechanically-driven mechanisms. We provided spatio-temporal estimations of matrix damage, which are in accordance with our experimental findings (S1 Supplementary Material Subsection S2.1, S2 Fig) and the current literature reports. With the chosen concentration of pro-inflammatory cytokine IL-1 (1 ng/ml, representing moderate inflammation), the ECM was degraded more dramatically over larger areas (near free surfaces and below the lesion) compared to the effect of mechanical loading alone after injury. The shear strain-modulated biomechanical degeneration that occurred only near the lesion was rapid, especially during 1–3 days following injury.

Biochemical degradation

The intact reference model predicted slightly smaller FCD loss (38%) than the injury model both near (44%) and away from (44%) the lesion (Fig 4). This might be associated with the increased surface area available (lesion surfaces) for cytokines to diffuse into the tissue in the injury model. Likewise, this geometrical alteration may explain the increase in FCD loss near the edge of the lesion (peak FCD loss 58%) and just below the lesion compared to regions away from the defect. Therefore, the mere presence of lesions might affect biochemical degradation to some degree.

For comparison of biochemical model and experiments, Kar et al. [10] already provided data showing sound agreement between their model and experimental IL-1-induced degradation of young calf cartilage in terms of time-dependent aggrecan loss to the culture medium [40]. In addition, Patwari et al. [35] reported an average total GAG loss of 280 μg/ml per disk (600% more compared to untreated controls) over 6-day treatment using 10 ng/mL IL-1, i.e., 10-fold higher cytokine concentration compared to that used in our model. Similarly, Stevens et al. [36] found that cartilage disks treated with 10 ng/ml of IL-1 for five days released aggrecan 2.6 times the amount released by control cartilage disks. These findings are supportive of the predictions of our biochemical model using 1 ng/ml of IL-1. FCD losses in bulk of the tissue (see S1 Supplementary Material Subsection S2.2, S3 Fig) were 31% (day 5) and 39% (day 6), while near the lesion they were 56% (day 5) and 67% (day 6). Furthermore, when we decreased the IL-1 levels to physiological levels in synovial fluid (0.01 ng/ml [89]) the predicted FCD loss decreased drastically (at day 4, -70% near the lesion and -83% away from the lesion) compared to treatment with 1 ng/ml of IL-1 (S1 Supplementary Material Subsection S2.3, S4 Fig).

Biomechanical degradation

Biomechanical degradation models predicted localized loss of FCD near the lesion. These predictions were in good agreement with our experimental findings; loss of FCD near lesions was significantly higher than further away from the lesions (S1 Supplementary Material Subsection S2.1, S2 Fig). These experimental data and the biomechanical degradation model together suggest that the presence of lesions themselves is the primary reason for biomechanically-induced loss of FCD since dynamic loading itself did not degrade the ECM of intact cartilage. In support of these model simulations, we found experimentally that significantly more FCD was lost locally near lesions than away from them, while the average FCD concentrations were similar (1) away from lesions, (2) in freely-swollen cartilage disks, and (3) in uninjured dynamically loaded disks (see S1 Supplementary Material Subsection 2.1, S2 Fig). In addition, inclusion of the base degeneration rate (spontaneous degeneration of free swelling samples) led to a better fit with the experimental findings compared to the original injury model (Fig 7B).

Research quantifying the loss of GAG over time associated with the combination of initial impact injury followed by physiologically relevant dynamic loading are scarce [52] and, to the best of our knowledge, the current literature lacks data of localized GAG loss under aforementioned conditions. In contrast, studies using injurious loading alone [83,87,90] report rapid early GAG depletion. Specifically, Thibault et al. [86] found significantly higher GAG release rates during the first four days after loading compared to unloaded controls. Additionally, Rolauffs et al. [83] reported rapid GAG loss (~20%) by 48h after injurious loading (using the same injury-loading protocol as in the current study), comparable to our predictions obtained with biomechanical degradation rate Dmech in Eq (6). Similarly fast loss of GAG was reported by DiMicco et al. [87] and Mohanraj et al. [88] in mechanically injured samples. Mohanraj et al. found early damage in cartilage tissue analog constructs after already 12, 24, 48, and 120 hours after injury (50% or 75% strain amplitude with rate 50%/s). Furthermore, they found a 2-3-fold increase in bulk PG matrix loss in injured constructs compared to uninjured controls by 48h after injury. This was predicted quite well by our biomechanical degradation model with base degeneration at t = 2 d (6-7-fold increase in localized FCD loss near the lesion compared to away from the lesion, see Subsection Results: Simulation of biomechanical degradation, and 2.2-fold increase in bulk FCD loss in the injury model versus intact reference model, see S1 Supplementary Material Subsection S2.2, S3 Fig).

Combined biochemical and biomechanical degradation

The intact reference model combining both biomechanical and biochemical degradation mechanisms exhibited the same FCD loss as the biochemical model (38%, Figs 4A and 6A). The injury model revealed marked depletion of FCD near the lesion compared to regions away from the lesion. Moreover, the FCD loss predicted by the combined degradation model and that by the biochemical model was essentially identical away from the defect. This suggests that biomechanical degradation does not play a marked role away from lesions even in the case when those regions are already affected by biochemical perturbations.

The maximum loss of FCD at the lesion site predicted by the combined biochemically and biomechanically-driven model was slightly smaller compared to model driven only by biomechanical degradation mechanism (for example, at timepoint t = 4 d the combined and biomechanical models predicted 81% and 89% of peak FCD loss, respectively, see Fig 7C). The reason for this is that the combined model considers both passive diffusion of aggrecan and synthesis of new aggrecan, and their combined effect decreases the net rate of aggrecan loss. In contrast, the biomechanical model does not include aggrecan synthesis and, thus, does not include mechanisms for generation of new FCD.

Experimental literature that includes the known extent of quantified cartilage cracks along with biochemical degradation (e.g., IL-1-induced) and dynamic loading is lacking. However, Li et al. [52] reported 40% total GAG loss after an 8-day treatment with 10% strain amplitude dynamic loading in the combined presence of IL-6, IL-6 soluble receptor (sIL-6R) and TNFα. Similarly, Sui et al. [34] reported 45% total GAG loss after a 6-day treatment with these same inflammatory cytokines and injurious loading just prior to the cytokine culture. These results are in accordance with our numerical predictions in bulk of the tissue (see S1 Supplementary Material Subsection S2.2, S3 Fig) with only 1 ng/ml of IL-1 after 6 days (42% FCD loss) and 8 days (61%) of simulation. Furthermore, the predicted average FCD losses near the lesion were 77% (day 6) and 95% (day 8) (Fig 7C).

Limitations

Several limitations in this study are worth mentioning. Characterizing the temporal dynamics associated with combinations of inflammatory cytokines following mechanical injury is a difficult task. Due to lack of quantitative data about the interplay of different cytokines, we modeled only the net effect of IL-1, ADAMTS-4,5, and TIMP-3. In the end, a model including the whole range of pro- and anti-inflammatory cytokines with their individual activation characteristics would most likely be computationally very demanding. However, such a model might not predict substantially different spatio-temporal matrix degradation compared to the models used in this study. Lack of quantitative experimental data also makes it difficult to include cytokine-mechanical loading-interactions into our models. Although some evidence exists that injurious loading stimulates a few selected signaling pathways also stimulated by IL-1 and TNFα [29,91], modeling the dynamics of these pathways and their subsequent effect on cartilage composition was out of scope of this study. Also, the predictions of combined effect of cytokines and mechanical injury on the localized FCD loss especially at the lesion site were not compared to experimental findings due to lack of data in the literature.

Several model parameters (such as effective diffusion and catalytic rate coefficients of aggrecanases) need calibration and new experimental validation, as also recognized elsewhere [10]. Another drawback of the biochemical model is the exclusion of electric charge-based partitioning (Gibbs–Donnan effect) of IL-1 at the cartilage–medium interface. In fact, an updated version of a cytokine diffusion model [9] included the Nernst–Planck equations which would account for the electric charges of IL-1 (valency of -1 [9]) and Donnan partitioning. We are aware that the charge of cytokines may play a role in cytokine diffusion, but this was beyond the scope of the present study. We presume that Donnan partitioning associated with the negative charge of both IL-1 and GAGs may result in slightly smaller FCD loss if this effect was included in our model.

Another drawback is associated with modeling in 2D rather than 3D. The latter is especially important for future applications with patient-specific PTOA predictions, which in turn would need validation of model parameters in human instead of bovine cartilage. We chose to model degradation of bovine tissue in 2D to develop and investigate the computational concept of combined biochemical and biomechanical degradation mechanisms. At this step of model development, the concept was tested using a simple geometry in tandem with repeatable and relatively inexpensive experimental tissue degradation models. However, we acknowledge that the diffusion of cytokines in 3D cartilage geometries most likely differs from the 2D situation. Moreover, the complex 3D state of strain might also deviate from simple 2D case and thus have an impact on biomechanical damage. In addition, we chose not to include crack propagation since this would be challenging to validate experimentally.

Although chondrocyte apoptosis is appreciable in cartilage samples treated with catabolic cytokines [40] and in highly sheared/compressed regions [79,80], we chose not to include cell death explicitly in our present models. These processes would introduce new time and depth-dependent cell death parameters in need of validation. In the model of this study, decreasing cell concentration (1) decreases local aggrecan biosynthesis (thus, decreases FCD content) but in turn (2) leads to smaller amount of available IL-1 receptors (C* in Eq (2) decreases) and subsequently decreases expression of aggrecanases (thus, increases FCD content). The current model parameters should be critically scrutinized, since now in a preliminary model with decreased initial cell concentration the net effect of these competing processes increases FCD content compared to a model with unaltered amount of cells.

The biomechanical model without the base degeneration reached a near-equilibrium state in terms of FCD loss. We acknowledge that this might not be a real-life scenario; for example, FCD loss might lead to higher fluid velocity [49], altered shear strain threshold ( in Eq (6)), or crack propagation during dynamic loading which in turn might lead to accelerated FCD loss. However, the model with the selected assumptions led to consistent results with the experiments.

In these mechanobiological models, we concentrated only on the non-fibrillar (i.e., proteoglycan) matrix loss and purposely did not include the collagen degradation present in the original biochemical model [10]. Based on earlier animal model experiments [55], the loss of collagen is negligible over short periods of time post-injury. Also, in our preliminary models, we observed that inclusion of the collagen degradation did not have a substantial effect on tissue shear strains and thus on biomechanical degradation.

Future

The functional significance of mechanobiological model predictions is yet to be fully appreciated. In the future, computer models like the ones used in the current study could facilitate prediction of early PTOA progression, recognition of high-risk chondral lesions, and simulation of treatment and rehabilitation outcomes. However, these aims need more experimental data on cartilage matrix loss at early times post-trauma via inflammatory mediators and shear strain-modulated tissue failure. Such data could also provide new information of possible “therapeutic window of opportunity” for effective drug interventions of PTOA [40,92]. Therefore, we plan to conduct new experiments with simultaneous biochemical and biomechanical degradations after injurious loading. This will be done at several timepoints to calibrate the degeneration rate parameters and to validate the mechanobiological algorithm.

The computational models should also be brought from tissue to joint level in order to predict PTOA progression for clinical relevance. This could be possible with synovial fluid or serum data on biochemical mediators [93,94] combined with subject-specific joint models obtained via medical imaging [18]. For example, it would be interesting to include predictions of inflammation and biomechanically-driven degradation of ACL injured patients. Furthermore, besides joint trauma, increased levels of inflammatory mediators are suggested to play a role in the phenotype of osteoarthritis induced by obesity [95].

The current study also demonstrates the high need for fast and robust programming environments or software which could handle both inflammatory and biomechanical aspects simultaneously. This would reduce the strenuous coupling of several solvers and make the approach less complex and faster. Currently, the COMSOL Multiphysics can run the biochemical simulation alone in two minutes, but due to the complex biomechanical material model in ABAQUS the full simulation with both degradation mechanisms takes three days.

Conclusions

Biochemical cues and abnormal mechanical environment are potential contributors linking knee trauma to progressive cartilage degeneration. Our numerical models considered these two cartilage-degrading mechanisms simultaneously for the first time in finite element models of mechanically injured cartilage, providing novel mechanobiological insights into PTOA progression. The model predictions were supported by current literature; cartilage subjected to pro-inflammatory stressors degraded near free surfaces, whereas in injured and dynamically loaded cartilage the FCD content was depleted especially near tissue lesions. Moreover, we suggest that the occurrence of chondral lesions (1) might affect the cytokine diffusion-driven degradation to some extent, and (2) could be the trigger for further accelerated biomechanical damage. In the future, mechanobiological models could evolve to in silico analysis tools that provide clinicians decision-making guidance for treatment of PTOA. This could lead to increased quality of life of patients via forestalling the onset of PTOA.

Supporting information

S1 Supplementary Material. Electronic supplementary material.

The supplementary material containing more detailed information about the underlying experimental findings, biochemical and biomechanical models, finite element mesh density, physiological levels of interleukin-1, and estimations of degraded cartilage area over time.

https://doi.org/10.1371/journal.pcbi.1007998.s001

(DOCX)

S1 Fig. Mesh sensitivity analysis.

Average fixed charge density (FCD) loss at time t = 4 d with combined biochemical and biomechanical degradation with A) 499, B) 697, C) 918, D) 1217, and E) 2029 elements. The mesh with 918 elements was chosen, as increasing the mesh density from this did not yield quantitatively nor qualitatively different predictions for average FCD loss anymore.

https://doi.org/10.1371/journal.pcbi.1007998.s002

(TIF)

S2 Fig. Experimental findings in the absence of cytokine interleukin-1 suggest that cartilage loses more of its fixed charge density content near lesions compared to areas away from lesions.

Quantification of average optical density (OD) and fixed charge density (FCD) loss in freely-swollen control samples, uninjured dynamically loaded samples and injured dynamically loaded samples treated for 7 or 12 days. These biomechanical degradation experiments were carried out by Orozco et al. [49]. In injured samples, ODs were calculated as an average within 0.1 mm (±10%) from a lesion. In regions away from lesions, ODs were calculated as an average from a 0.45 mm (±10%) wide and 15 mm (±10%) thick surface region at the midway between lesion and sample edge. In freely-swollen and uninjured dynamically loaded samples, ODs were calculated as an average from a 0.45 mm (±10%) wide and 15 mm (±10%) thick surface region at the middle of samples. A) At day 0, average ODs between the freely-swollen and injured samples were similar (one-way ANOVA, p = 0.068). B) Average OD decreases significantly over time both away from lesions (one-way ANOVA, p = 5.6⋅10−4) and near lesions (p = 2.6⋅10−7; the figure shows Tukey's honestly significant difference (HSD) test results). C) At days 0, 7 and 12, average ODs were statistically similar between freely-swollen, uninjured dynamically loaded, and away from the lesion -groups (one-way ANOVA, p = 0.105 for day 0, p = 0.213 for day 7, p = 0.416 for day 12). However, within these treatment groups the average ODs decreased in time, especially in the away from the lesion -group which exhibited statistically significant decrease in OD (one-way ANOVA, p = 0.065 for freely-swollen samples, p = 0.075 for uninjured dynamically loaded samples, and p = 5.6⋅10−4 for injured dynamically loaded samples away from lesions; the figure shows Tukey’s HSD test results). D) FCD losses (calculated from average ODs at day 7 and 12 compared to day 0) near lesions were significantly greater than away from lesions (the figure shows dependent samples t-test results). Box plots display values as range (brackets), interquartiles and median (solid bars).

https://doi.org/10.1371/journal.pcbi.1007998.s003

(TIF)

S3 Fig. Average cartilage degradation in the whole explant geometry.

Simulated average bulk fixed charge density (FCD) losses in the whole explant with biochemical, biomechanical (with and without base degeneration, see S1 Supplementary Material Subsection S1.8) and combined biochemical and biomechanical degradation models.

https://doi.org/10.1371/journal.pcbi.1007998.s004

(TIF)

S4 Fig. Cartilage degradation with physiological concentrations of interleukin-1.

A) Predicted biochemically driven fixed charge density (FCD) losses A) over time under moderate inflammation (1 ng/ml of exogenous interleukin-1) and physiological levels (0.01 ng/ml) of pro-inflammatory mediators. Still images at day t = 4 d show B) markedly higher matrix losses with moderate inflammation compared to C) physiological conditions.

https://doi.org/10.1371/journal.pcbi.1007998.s005

(TIF)

S5 Fig. Still image preview of animations of cartilage degradation.

Animations of predicted fixed charge density (FCD) losses over 21 days of biochemical, biomechanical (without or with base degeneration), and combined biochemical and biomechanical degradation.

https://doi.org/10.1371/journal.pcbi.1007998.s006

(TIF)

S1 Table. Compositional and material parameters for the fibril-reinforced porohyperelastic swelling (FRPHES) material model.

z is the normalized distance from the cartilage surface (z = 0) to the bottom (z = 1).

https://doi.org/10.1371/journal.pcbi.1007998.s007

(PNG)

S1 Animation. Animation of simulated biochemical degradation.

https://doi.org/10.1371/journal.pcbi.1007998.s008

(GIF)

S2 Animation. Animation of simulated biomechanical degradation without base degeneration.

https://doi.org/10.1371/journal.pcbi.1007998.s009

(GIF)

S3 Animation. Animation of simulated biomechanical degradation with base degeneration.

https://doi.org/10.1371/journal.pcbi.1007998.s010

(GIF)

S4 Animation. Animation of simulated combined biochemical and biomechanical degradation.

https://doi.org/10.1371/journal.pcbi.1007998.s011

(GIF)

Acknowledgments

The authors appreciate the support of the University of Eastern Finland and the Massachusetts Institute of Technology to conduct this study. CSC–IT center for Science Ltd., Finland, is acknowledged for providing the modeling software. Patrick Grahn, D.Sc. (Tech.), from COMSOL Multiphysics Support is acknowledged for technical support.

Citation: Eskelinen ASA, Tanska P, Florea C, Orozco GA, Julkunen P, Grodzinsky AJ, et al. (2020) Mechanobiological model for simulation of injured cartilage degradation via pro-inflammatory cytokines and mechanical stimulus. PLoS Comput Biol 16(6): e1007998. https://doi.org/10.1371/journal.pcbi.1007998

References

1. Nelson F, Billinghurst R, Pidoux R, Reiner A, Langworthy M, McDermott M, et al. Early post-traumatic osteoarthritis-like changes in human articular cartilage following rupture of the anterior cruciate ligament. Osteoarthr Cartil. 2006;14: 114–119. pmid:16242972

2. Kiapour AM, Murray MM. Basic science of anterior cruciate ligament injury and repair. Bone Joint Res. 2014;3: 20–31. pmid:24497504

3. Lieberthal J, Sambamurthy N, Scanzello C. Inflammation in Joint Injury and Post-Traumatic Osteoarthritis. Osteoarthr Cartil. 2015;23: 1825–1834. pmid:26521728

4. Carbone A, Rodeo S. Review of current understanding of post-traumatic osteoarthritis resulting from sports injuries. J Orthop Res. 2017;35: 397–405. pmid:27306867

5. Mahmoudian A, Van Assche D, Herzog W, Luyten FP. Towards secondary prevention of early knee osteoarthritis. RMD Open. 2018;4: 1–12. pmid:30167325

6. Berenbaum F. Osteoarthritis as an inflammatory disease (osteoarthritis is not osteoarthrosis!). Osteoarthr Cartil. 2013;21: 16–21. pmid:23194896

7. Andriacchi T, Favre J, Erhart-Hledik J, Chu C. A Systems View of Risk Factors for Knee Osteoarthritis Reveals Insights into the Pathogenesis of the Disease. Ann Biomed Eng. 2015;43: 376–387. pmid:25224078

8. Loeser RF, Goldring SR, Scanzello CR, Goldring MB. Osteoarthritis: A disease of the joint as an organ. Arthritis Rheum. 2012;64: 1697–1707. pmid:22392533

9. Kar S, Smith DW, Gardiner BS, Grodzinsky AJ. Systems based study of the therapeutic potential of small charged molecules for the inhibition of IL-1 mediated cartilage degradation. PLoS One. 2016;11: 1–38. pmid:27977731

10. Kar S, Smith DW, Gardiner BS, Li Y, Wang Y, Grodzinsky AJ. Modeling IL-1 induced degradation of articular cartilage. Arch Biochem Biophys. 2016;594: 37–53. pmid:26874194

11. Krishnan Y, Grodzinsky AJ. Cartilage diseases. Matrix Biol. 2018;71–72: 51–69. pmid:29803938

12. Guilak F. Biomechanical factors in osteoarthritis. Best Pract Res Clin Rheumatol. 2011;25: 815–823. pmid:22265263

13. Mündermann A, Dyrby CO, Andriacchi TP. Secondary gait changes in patients with medial compartment knee osteoarthritis: Increased load at the ankle, knee, and hip during walking. Arthritis Rheum. 2005;52: 2835–2844. pmid:16145666

14. Harding GT, Hubley-Kozey CL, Dunbar MJ, Stanish WD, Astephen Wilson JL. Body mass index affects knee joint mechanics during gait differently with and without moderate knee osteoarthritis. Osteoarthr Cartil. 2012;20: 1234–1242. pmid:22902710

15. Meireles S, Wesseling M, Smith CR, Thelen DG, Verschueren S, Jonkers I. Medial knee loading is altered in subjects with early osteoarthritis during gait but not during step-up-and-over task. PLoS One. 2017;12: e0187583. pmid:29117248

16. Zevenbergen L, Gsell W, Chan D, Sloten J, Himmelreich U, Neu C, et al. Functional assessment of strains around a full-thickness and critical sized articular cartilage defect under compressive loading using MRI. Osteoarthr Cartil. 2018;26: 1710–1721. pmid:30195045

17. Eskelinen A, Mononen M, Venäläinen M, Korhonen R, Tanska P. Maximum shear strain-based algorithm can predict proteoglycan loss in damaged articular cartilage. Biomech Model Mechanobiol. 2019;18: 753–778. pmid:30631999

18. Myller KAH, Korhonen RK, Töyräs J, Salo J, Jurvelin JS, Venäläinen MS. Computational evaluation of altered biomechanics related to articular cartilage lesions observed in vivo. J Orthop Res. 2019; 1–30. pmid:30839123

19. Párraga Quiroga J, Wilson W, Ito K, van Donkelaar C. The effect of loading rate on the development of early damage in articular cartilage. Biomech Model Mechanobiol. 2017;16: 263–273. pmid:27514541

20. Liukkonen M, Mononen M, Klets O, Arokoski J, Saarakkala S, Korhonen R. Simulation of subject-specific progression of knee osteoarthritis and comparison to experimental follow-up data: Data from the osteoarthritis initiative. Sci Rep. 2017;7: 9177. pmid:28835668

21. Mononen M, Liukkonen M, Korhonen R. Utilizing Atlas-Based Modeling to Predict Knee Joint Cartilage Degeneration: Data from the Osteoarthritis Initiative. Ann Biomed Eng. 2019;47: 813–825. pmid:30547410

22. Graham JM, Ayati BP, Ding L, Ramakrishnan PS, Martin JA. Reaction-Diffusion-Delay Model for EPO/TNF-α Interaction in articular cartilage lesion abatement. Biol Direct. 2012;7: 9. pmid:22353555

23. Guilak F, Pferdehirt L, Ross AK, Choi Y, Collins KH, Nims RJ, et al. Designer stem cells: Genome engineering and the next generation of cell‐based therapies. J Orthop Res. 2019; jor.24304. pmid:30977548

24. Bajpayee AG, De la Vega RE, Scheu M, Varady NH, Yannatos IA, Brown LA, et al. Sustained intra-cartilage delivery of low dose dexamethasone using a cationic carrier for treatment of post traumatic osteoarthritis. Eur Cell Mater. 2017;34: 341–364. pmid:29205258

25. Grodzinsky AJ, Wang Y, Kakar S, Vrahas MS, Evans CH. Intra-articular dexamethasone to inhibit the development of post-traumatic osteoarthritis. J Orthop Res. 2017;35: 406–411. pmid:27176565

26. Kapitanov G, Wang X, Ayati B, Brouillette M, Martin J. Linking Cellular and Mechanical Processes in Articular Cartilage Lesion Formation: A Mathematical Model. Front Bioeng Biotechnol. 2016;4: 1–22.

27. Murakami S, Muneta T, Furuya K, Saito I, Miyasaka N, Yamamoto H. Immunohistologic Analysis of Synovium in Infrapatellar Fat Pad After Anterior Cruciate Ligament Injury. Am J Sports Med. 1995;23: 763–768. pmid:8600748

28. Kennedy A, Fearon U, Veale DJ, Godson C. Macrophages in synovial inflammation. Front Immunol. 2011;2: 1–9.

29. Goldring MB, Otero M, Tsuchimochi K, Ijiri K, Li Y. Defining the roles of inflammatory and anabolic cytokines in cartilage metabolism. Ann Rheum Dis. 2008;67: 75–82. pmid:19022820

30. Boehme KA, Rolauffs B. Onset and Progression of Human Osteoarthritis—Can Growth Factors, Inflammatory Cytokines, or Differential miRNA Expression Concomitantly Induce Proliferation, ECM Degradation, and Inflammation in Articular Cartilage? Int J Mol Sci. 2018;19: 2282. pmid:30081513

31. Bevill SL, Boyer KA, Andriacchi TP. The Regional Sensitivity of Chondrocyte Gene Expression to Coactive Mechanical Load and Exogenous TNF-α Stimuli. J Biomech Eng. 2014;136: 091005. pmid:24976081

32. Kothari P, Pestana R, Mesraoua R, Elchaki R, Khan KMF, Dannenberg AJ, et al. IL-6-Mediated Induction of Matrix Metalloproteinase-9 Is Modulated by JAK-Dependent IL-10 Expression in Macrophages. J Immunol. 2014;192: 349–357. pmid:24285838

33. Lotz M, Terkeltaub R, Villiger PM. Cartilage and joint inflammation. Regulation of IL-8 expression by human articular chondrocytes. J Immunol. 1992;148: 466–473. pmid:1729366

34. Sui Y, Lee JH, DiMicco MA, Vanderploeg EJ, Blake SM, Hung HH, et al. Mechanical injury potentiates proteoglycan catabolism induced by interleukin-6 with soluble interleukin-6 receptor and tumor necrosis factor α in immature bovine and adult human articular cartilage. Arthritis Rheum. 2009;60: 2985–2996. pmid:19790045

35. Patwari P, Gao G, Lee JH, Grodzinsky AJ, Sandy JD. Analysis of ADAMTS4 and MT4-MMP indicates that both are involved in aggrecanolysis in interleukin-1-treated bovine cartilage. Osteoarthr Cartil. 2005;13: 269–277. pmid:15780640

36. Stevens AL, Wishnok JS, White FM, Grodzinsky AJ, Tannenbaum SR. Mechanical Injury and Cytokines Cause Loss of Cartilage Integrity and Upregulate Proteins Associated with Catabolism, Immunity, Inflammation, and Repair. Mol Cell Proteomics. 2009;8: 1475–1489. pmid:19196708

37. Troeberg L, Nagase H. Proteases involved in cartilage matrix degradation in osteoarthritis. Biochim Biophys Acta—Proteins Proteomics. 2012;1824: 133–145. pmid:21777704

38. Kashiwagi M, Tortorella M, Nagase H, Brew K. TIMP-3 Is a Potent Inhibitor of Aggrecanase 1 (ADAM-TS4) and aggrecanase 2 (ADAM-TS5). J Biol Chem. 2001;276: 12501–12504. pmid:11278243

39. Su S, Grover J, Roughley PJ, DiBattista JA, Martel-Pelletier J, Pelletier JP, et al. Expression of the tissue inhibitor of metalloproteinases (TIMP) gene family in normal and osteoarthritic joints. Rheumatol Int. 1999;18: 183–191. pmid:10399794

40. Li Y, Wang Y, Chubinskaya S, Schoeberl B, Florine E, Kopesky P, et al. Effects of insulin-like growth factor-1 and dexamethasone on cytokine-challenged cartilage: Relevance to post-traumatic osteoarthritis. Osteoarthr Cartil. 2015;23: 266–274. pmid:25450855

41. Faisal TR, Adouni M, Dhaher YY. The effect of fibrillar degradation on the mechanics of articular cartilage: a computational model. Biomech Model Mechanobiol. 2019. pmid:30604303

42. Felson DT. Osteoarthritis as a disease of mechanics. Osteoarthr Cartil. 2013;21: 10–15. pmid:23041436

43. Mahmoudian A, van Dieen JH, Bruijn SM, Baert IA, Faber GS, Luyten FP, et al. Varus thrust in women with early medial knee osteoarthritis and its relation with the external knee adduction moment. Clin Biomech. 2016;39: 109–114. pmid:27744006

44. Sadeghi H, Lawless BM, Espino DM, Shepherd DET. Effect of frequency on crack growth in articular cartilage. J Mech Behav Biomed Mater. 2018;77: 40–46. pmid:28888932

45. Morel V, Berutto C, Quinn TM. Effects of damage in the articular surface on the cartilage response to injurious compression in vitro. J Biomech. 2006;39: 924–930. pmid:16488230

46. Mononen M, Tanska P, Isaksson H, Korhonen R. A novel method to simulate the progression of collagen degeneration of cartilage in the knee: Data from the osteoarthritis initiative. Sci Rep. 2016;6: 21415. pmid:26906749

47. Hosseini S, Wilson W, Ito K, van Donkelaar C. A numerical model to study mechanically induced initiation and progression of damage in articular cartilage. Osteoarthr Cartil. 2014;22: 95–103. pmid:24185112

48. Stender M, Regueiro R, Klisch S, Ferguson V. An Equilibrium Constitutive Model of Anisotropic Cartilage Damage to Elucidate Mechanisms of Damage Initiation and Progression. J Biomech Eng. 2015;137: 081010. pmid:26043366

49. Orozco GA, Tanska P, Florea C, Grodzinsky AJ, Korhonen RK. A novel mechanobiological model can predict how physiologically relevant dynamic loading causes proteoglycan loss in mechanically injured articular cartilage. Sci Rep. 2018;8: 1–16.

50. Manzano S, Manzano R, Doblaré M, Doweidar MH. Altered swelling and ion fluxes in articular cartilage as a biomarker in osteoarthritis and joint immobilization: A computational analysis. J R Soc Interface. 2014;12: 20141090. pmid:25392400

51. Mononen M, Tanska P, Isaksson H, Korhonen R. New algorithm for simulation of proteoglycan loss and collagen degeneration in the knee joint: Data from the osteoarthritis initiative. J Orthop Res. 2018;36: 1673–1683. pmid:29150953

52. Li Y, Frank E, Wang Y, Chubinskaya S, Huang H-H, Grodzinsky A. Moderate Dynamic Compression Inhibits Pro-Catabolic Response of Cartilage to Mechanical Injury, TNF-α and IL-6, but Accentuates Degradation Above a Strain Threshold. Osteoarthr Cartil. 2013;21. pmid:24007885

53. Wilson W, Driessen NJB, van Donkelaar CC, Ito K. Prediction of collagen orientation in articular cartilage by a collagen remodeling algorithm. Osteoarthr Cartil. 2006;14: 1196–1202. pmid:16797194

54. Tanska P, Julkunen P, Korhonen R. A computational algorithm to simulate disorganization of collagen network in injured articular cartilage. Biomech Model Mechanobiol. 2018;17: 689–699. pmid:29177932

55. Mäkelä JTA, Rezaeian ZS, Mikkonen S, Madden R, Han SK, Jurvelin JS, et al. Site-dependent changes in structure and function of lapine articular cartilage 4 weeks after anterior cruciate ligament transection. Osteoarthr Cartil. 2014;22: 869–878. pmid:24769230

56. Grenier S, Bhargava M, Torzilli P. An in vitro model for the pathological degradation of articular cartilage in osteoarthritis. J Biomech. 2014;47: 645–652. pmid:24360770

57. Quinn TM, Maung AA, Grodzinsky AJ, Hunziker EB, Sandy JD. Physical and biological regulation of proteoglycan turnover around chondrocytes in cartilage explants. Implications for tissue degradation and repair. Ann N Y Acad Sci. 1999;878: 420–441. pmid:10415746

58. Patwari P, Cheng DM, Cole AA, Kuettner KE, Grodzinsky AJ. Analysis of the relationship between peak stress and proteoglycan loss following injurious compression of human post-mortem knee and ankle cartilage. Biomech Model Mechanobiol. 2007;6: 83–89. pmid:16715319

59. Kiviranta I, Jurvelin J, Säämänen A, Helminen H. Microspectrophotometric quantitation of glycosaminoglycans in articular cartilage sections stained with Safranin O. Histochemistry. 1985;82: 249–255. pmid:2581923

60. Király K, Lapveteläinen T, Arokoski J, Törrönen K, Módis L, Kiviranta I, et al. Application of selected cationic dyes for the semiquantitative estimation of glycosaminoglycans in histological sections of articular cartilage by microspectrophotometry. Histochem J. 1996;28: 577–590. pmid:8894661

61. Dileo MV., Kellum JA, Federspiel WJ. A simple mathematical model of cytokine capture using a hemoadsorption device. Ann Biomed Eng. 2009;37: 222–229. pmid:18949559

62. Comper WD, Williams RP. Hydrodynamics of concentrated proteoglycan solutions. J Biol Chem. 1987;262: 13464–13471. pmid:3654623

63. DiMicco MA, Sah RL. Dependence of cartilage matrix composition on biosynthesis, diffusion, and reaction. Transp Porous Media. 2003;50: 57–73.

64. Catt CJ, Schuurman W, Sengers BG, van Weeren PR, Dhert WJA, Please CP, et al. Mathematical modelling of tissue formation in chondrocyte filter cultures. Eur Cells Mater. 2011;22: 377–392. pmid:22179936

65. Vempati P, Mac Gabhann F, Popel AS. Quantifying the proteolytic release of extracellular matrix-sequestered VEGF with a computational model. PLoS One. 2010;5. pmid:20686621

66. Maroudas A. Glycosaminoglycan turn-over in articular cartilage. Philos Trans R Soc Lond B Biol Sci. 1975;271: 293–313. pmid:239424

67. Häuselmann H, Oppliger L, Michel B, Stefanovic-Racic M, Evans C. Nitric oxide and proteoglycan biosynthesis by human articular chondrocytes in alginate culture. FEBS Lett. 1994;352: 361–364. pmid:7926002

68. Häuselmann H, Flechtenmacher J, Michal L, Thonar E, Shinmei M, Kuettner K, et al. The superficial layer of human articular cartilage is more susceptible to interleukin-1-induced damage than the deeper layers. Arthritis Rheum. 1996;39: 478–488. pmid:8607897

69. Nagase H, Kashiwagi M. Aggrecanases and cartilage matrix degradation. Arthritis Res Ther. 2003;5: 94–103. pmid:12718749

70. Klein TJ, Chaudhry M, Bae WC, Sah RL. Depth-dependent biomechanical and biochemical properties of fetal, newborn, and tissue-engineered articular cartilage. J Biomech. 2007;40: 182–190. pmid:16387310

71. Wittwer AJ, Hills RL, Keith RH, Munie GE, Arner EC, Anglin CP, et al. Substrate-dependent inhibition kinetics of an active site-directed inhibitor of ADAMTS-4 (aggrecanase 1). Biochemistry. 2007;46: 6393–6401. pmid:17487981

72. Yamamoto K, Owen K, Parker AE, Scilabra SD, Dudhia J, Strickland DK, et al. Low density lipoprotein receptor-related protein 1 (LRP1)-mediated endocytic clearance of a disintegrin and metalloproteinase with thrombospondin motifs-4 (ADAMTS-4): Functional differences of non-catalytic domains of ADAMTS-4 and ADAMTS-5 in LRP1 binding. J Biol Chem. 2014;289: 6462–6474. pmid:24474687

73. Soh S, Byrska M, Kandere-Grzybowska K, Grzybowski B. Reaction -Diffusion Systems in Intracellular Molecular Transport and Control. Angew Chem Int Ed Engl. 2010;49: 4170–4198. pmid:20518023

74. Chan PS, Caron JP, Orth MW. Short-term gene expression changes in cartilage explants stimulated with interleukin 1β plus glucosamine and chondroitin sulfate. J Rheumatol. 2006;33: 1329–1340. pmid:16821268

75. Jadin KD, Wong BL, Bae WC, Li KW, Williamson AK, Schumacher BL, et al. Depth-varying density and organization of chondrocytes in immature and mature bovine articular cartilage assessed by 3D imaging and analysis. J Histochem Cytochem. 2005;53: 1109–1119. pmid:15879579

76. Smith D, Gardiner B, Davidson J, Grodzinsky A. Computational model for the analysis of cartilage and cartilage tissue constructs. J Tissue Eng Regen Med. 2016;10: 334–347. pmid:23784936

77. Gratz K, Wong B, Bae W, Sah R. The Effects of Focal Articular Defects on Cartilage Contact Mechanics. J Orthop Res. 2009;27: 584–592. pmid:18979528

78. Quinn T, Allen R, Schalet B, Perumbuli P, Hunziker E. Matrix and cell injury due to sub-impact loading of adult bovine articular cartilage explants: Effects of strain rate and peak stress. J Orthop Res. 2001;19: 242–249. pmid:11347697

79. Patwari P, Gaschen V, James I, Berger E, Blake S, Lark M, et al. Ultrastructural quantification of cell death after injurious compression of bovine calf articular cartilage. Osteoarthr Cartil. 2004;12: 245–252. pmid:14972342

80. Bonnevie E, Delco M, Bartell L, Jasty N, Cohen I, Fortier L, et al. Microscale frictional strains determine chondrocyte fate in loaded cartilage. J Biomech. 2018;74: 72–78. pmid:29729853

81. Wilson W, van Donkelaar C, van Rietbergen B, Huiskes R. A fibril-reinforced poroviscoelastic swelling model for articular cartilage. J Biomech. 2005;38: 1195–1204. pmid:15863103

82. Saarakkala S, Julkunen P. Specificity of fourier transform infrared (FTIR) microspectroscopy to estimate depth-wise proteoglycan content in normal and osteoarthritic human articular cartilage. Cartilage. 2010;1: 262–269. pmid:26069557

83. Rolauffs B, Muehleman C, Li J, Kurz B, Kuettner K, Frank E, et al. Vulnerability of the superficial zone of immature articular cartilage to compressive injury. Arthritis Rheum. 2010;62: 3016–3027. pmid:20556809

84. Wilson W, Van Donkelaar CC, Van Rietbergen B, Huiskes R. A fibril-reinforced poroviscoelastic swelling model for articular cartilage. J Biomech. 2005;38: 1195–1204. pmid:15863103

85. Sibole S, Erdemir A. Chondrocyte deformations as a function of tibiofemoral joint loading predicted by a generalized high-throughput pipeline of multi-scale simulations. PLoS One. 2012;7: e37538. pmid:22649535

86. Thibault M, Poole A, Buschmann M. Cyclic compression of cartilage/bone explants in vitro leads to physical weakening, mechanical breakdown of collagen and release of matrix fragments. J Orthop Res. 2002;20: 1265–1273. pmid:12472239

87. DiMicco M, Patwari P, Siparsky P, Kumar S, Pratta M, Lark M, et al. Mechanisms and Kinetics of Glycosaminoglycan Release Following In Vitro Cartilage Injury. Arthritis Rheum. 2004;50: 840–848. pmid:15022326

88. Mohanraj B, Meloni G, Mauck R, Dodge G. A high-throughput model of post-traumatic osteoarthritis using engineered cartilage tissue analogs. Osteoarthr Cartil. 2014;22: 1282–1290. pmid:24999113

89. McNulty AL, Rothfusz NE, Leddy HA, Guilak F. Synovial fluid concentrations and relative potency of interleukin-1 alpha and beta in cartilage and meniscus degradation. J Orthop Res. 2013;31: 1039–1045. pmid:23483596

90. Kurz B, Jin M, Patwari P, Cheng D, Lark M, Grodzinsky A. Biosynthetic response and mechanical properties of articular cartilage after injurious compression. J Orthop Res. 2001;19: 1140–1146. pmid:11781016

91. Deschner J, Hofman CR, Piesco NP, Agarwal S. Signal transduction by mechanical strain in chondrocytes. Curr Opin Clin Nutr Metab Care. 2003;6: 289–293. pmid:12690261

92. Lotz MK. New developments in osteoarthritis. Posttraumatic osteoarthritis: Pathogenesis and pharmacological treatment options. Arthritis Res Ther. 2010;12. pmid:20602810

93. Catterall JB, Stabler T V, Flannery CR, Kraus VB. Changes in serum and synovial fluid biomarkers after acute injury (NCT00332254). Arthritis Res Ther. 2010;12: R229. pmid:21194441

94. Amano K, Huebner JL, Stabler T V, Tanaka M, Mcculloch CE, Lobach I, et al. Synovial Fluid Profile at the Time of Anterior Cruciate Ligament Reconstruction and Its Association With Cartilage Matrix Composition 3 Years After Surgery. Am J Sports Med. 2018;46: 890–899. pmid:29364702

95. Collins K, Hart D, Seerattan R, Reimer R, Herzog W. High-fat/high-sucrose diet-induced obesity results in joint-specific development of osteoarthritis-like degeneration in a rat model. Bone Joint Res. 2018;7: 274–281. pmid:29922445