The effect of fibrillar degradation on the mechanics of articular cartilage: a computational model

Abstract

The pathogenesis and pathophysiological underpinnings of cartilage degradation are not well understood. Either mechanically or enzymatically mediated degeneration at the fibril level can lead to acute focal injuries that will, overtime, cause significant cartilage degradation. Understanding the relationship between external loading and the basic molecular structure of cartilage requires establishing a connection between the fibril-level defects and its aggregate effect on cartilage. In this work, we provide a multiscale constitutive model of cartilage to elucidate the effect of two plausible fibril degradation mechanisms on the aggregate tissue: tropocollagen crosslink failure (β) and a generalized surface degradation (δ). Using our model, the mechanics of aggregate tissue shows differed yield stress and post-yield behavior after crosslink failure and surface degradation compared to intact cartilage, and the tissue-level aggregate behaviors are different from the fibrillar behaviors observed in the molecular dynamics simulations. We also compared the effect of fibrillar defects in terms of crosslink failure and surface degradation in different layers of cartilage within the macroscale tissue construct during a simulated nanoindentation test. Although the mechanical properties of cartilage tissue were largely contingent upon the mechanical properties of the fibril, the macroscale mechanics of cartilage tissue showed ~ 10% variation in yield strain (tissue yield strain: ~ 27 to ~ 37%) compared to fibrillar yield strain (fibrillar yield strain: ~ 16 to ~ 26%) for crosslink failure and ~ 7% difference for the surface degradation (yield strain variations at the tissue: ~ 30 to ~ 37% and fibril: ~ 24 to ~ 26%) at the superficial layer. The yield strain was further delayed in middle layers at least up to 30% irrespective of the failure mechanisms. The cartilage tissue appeared to withstand more strain than the fibrils. The degeneration mechanisms of fibril differentially influenced the aggregate mechanics of cartilage, and the deviation may be attributed to fiber–matrix interplay, depth-dependent fiber orientation and fibrillar defects with different degradation mechanisms. The understanding of the aggregate stress–strain behavior of cartilage tissue, cartilage degradation and its underlying biomechanical factors is important for developing engineering approaches and therapeutic interventions for cartilage pathologies.

Appendices

Appendix A

1.1 Multiscale modeling of cartilage strain energy

In this study, the evolution of plastic stress in the fibril is triggered by the effective fibril yield stress derived from fibril strain energy (Tang et al. 2009). The yield stress is connected to the fibril yield condition as obtained from MDS (i), the plastic strain rate (ii) and the flow resistance (iii) of the tissues as follows:

$$\left\{ {\begin{array}{*{20}l} {\sigma_{y}^{\text{eff}} = \frac{4}{3}\bar{I}_{4e} \frac{{\partial W_{\text{fl}} }}{{\partial \bar{I}_{4e} }} = \sigma_{{_{y(\beta )} }}^{\text{fl}} \,{\text{and}}\,\sigma_{{_{y(\delta )} }}^{\text{fl}} } \hfill & {(i)} \hfill \\ {\dot{\gamma } = \dot{\gamma }_{o} \left| {\frac{{\sigma_{\text{eff}} }}{\alpha \left( t \right)}} \right|^{{{1 \mathord{\left/ {\vphantom {1 p}} \right. \kern-0pt} p}}} sig\left( {\sigma_{y}^{\text{eff}} } \right)} \hfill & {(ii)} \hfill \\ {\alpha \left( t \right) = \int {h\dot{\gamma }f\left( \alpha \right)} } \hfill & {(iii)} \hfill \\ \end{array} } \right.$$

(6)

These relationships help describe the hierarchical coupling between nanoscale collagen degradation and material properties of fibrils (Fig. 6). The plastic parameters, $\dot{\gamma } = 0.01/\text{s}$ and p = 0.05, are fixed due to wide agreement in the literature with respect to cartilage’s plastic rate of deformation (Gasser and Holzapfel 2002). A collagen fiber, which is modeled as the fiber-reinforced composite, includes the descriptions of fibril and incompressible neo-Hookean matrix. The elastic strain energy of the fiber can eventually be expressed under extension as

$$\left\{ {\begin{array}{*{20}l} {W_{fbt} \left( {\bar{I}_{4} ,\bar{I}_{4e} } \right) = v_{\text{fl}} W_{\text{fl}} \left( {\bar{I}_{1ef} ,\bar{I}_{4e} } \right) + v_{ml} \left( {\frac{{\mu_{fm} }}{2}\left( {\bar{I}_{1f} - 3} \right)} \right)} \hfill \\ {W_{fs} \left( {\bar{I}_{1f} ,\bar{I}_{4} ,\bar{I}_{4e} } \right) = \frac{1}{2}\mu^{efffb} \left( {\bar{I}_{4e} } \right)\left( {\bar{I}_{1fb} - \bar{I}_{1f} } \right)} \hfill \\ \end{array} } \right.$$

(7)

The total elastic strain energy density of the fiber is therefore expressed as

$$W_{fb} \left( {\bar{I}_{1f} ,\bar{I}_{4} ,\bar{I}_{4e} } \right) = W_{fbt} \left( {\bar{I}_{4} ,\bar{I}_{4e} } \right) + W_{fs} \left( {\bar{I}_{1f} ,\bar{I}_{4} ,\bar{I}_{4e} } \right)$$

(8)

The framework is hierarchically used to model the soft tissue (cartilage) as the fiber-reinforced composite material as well, and the corresponding axial and shear strain energy are written as

$$\left\{ {\begin{array}{*{20}l} {W_{tt} \left( {\bar{I}_{4} ,\bar{I}_{4e} } \right) = v_{f} W_{fb} \left( {\bar{I}_{1f} ,\bar{I}_{4} ,\bar{I}_{4e} } \right) + v_{m} \left( {\frac{{\mu_{m} }}{2}\left( {\bar{I}_{1} - 3} \right)} \right)} \hfill \\ {W_{ts} \left( {\bar{I}_{1f} ,\bar{I}_{4} ,\bar{I}_{4e} } \right) = \frac{1}{2}\mu^{\text{eff}} \left( {\bar{I}_{4e} } \right)\left( {\bar{I}_{1} - \bar{I}_{1f} } \right)} \hfill \\ \end{array} } \right.$$

(9)

The total strain energy of the tissue is defined by

$$W_{t} \left( {\bar{I}_{1f} ,\bar{I}_{4} ,\bar{I}_{4e} } \right) = W_{tt} \left( {\bar{I}_{4} ,\bar{I}_{4e} } \right) + W_{ts} \left( {\bar{I}_{1f} ,\bar{I}_{4} ,\bar{I}_{4e} } \right)$$

(10)

The noticeable feature of the strain energy function is its dependence on both elastic stretch and the total deformation. This is due to the assumption that the plastic deformation occurs only in the fiber while the matrix materials always remain elastic. The total strain energy for the fiber (W_fb) and tissue (W_t) were therefore formulated mathematically combining axial and shear strains for the both, and eventually, the total stress σ^t was expressed with fibrillar σ^f and nonfibrillar $\sigma^{\text{nf}}$ stress tensors as shown in Eq. (4) in the main article. The plastic stress in Eq. (4) becomes dominant when the effective stress, $\sigma_{y}^{\text{eff}}$, is more than the fibril yield strength, which in this case varies depending on the crosslink failure and surface degradation. The elastic part of the deformation gradient is determined by ${\mathbf{F}}_{e} = {\mathbf{FF}}_{p}^{ - 1}$, where the plastic deformation gradient is $f\left( {\lambda_{\text{fp}} } \right)$. The stretch, λ_fp, is unity up to the yield point and only comes into play beyond the fiber yield strength. The plastic stretch depends on the plastic strain rate $\dot{\gamma }$, which ultimately controls the plastic part of the deformation gradient (Fig. 7).

Fig. 7

Multiscale hierarchical organization of articular cartilage

Appendix B

2.1 Nonlinear optimization of degenerated fibril

The nonlinear optimization function lsqnonlin in MATLAB (The MathWorks Inc., Natick, MA, USA) was employed with the trust region reflective algorithm, which minimized an objective function $f\left( {\mathbf{x}} \right)$ in a least square sense $\left( {\sum\nolimits_{n} {\left[ {f\left( {\mathbf{x}} \right)} \right]^{2} } } \right)$. The optimal objective function provides a good compromise between σ ^MDS_fl (represents either $\left( {\sigma_{{^{fl} }} \left( \beta \right)} \right)$ and $\left( {\sigma_{{^{fl} }} \left( \delta \right)} \right)$ depending on degradation mechanisms, crosslink failure and generalized surface degradation, respectively) and $\sigma_{\text{fl}}^{\text{SEF}}$ for the fibril along the fiber direction,

$$f\left( {\mathbf{x}} \right) = \frac{{\sigma_{\text{fl}}^{\text{SEF}} \left( {\mathbf{x}} \right) - \sigma_{\text{fl}}^{\text{MDS}} }}{{\sigma_{\text{fl}}^{\text{MDS}} }}$$

(11)

where ${\mathbf{x}}$ is a vector of the unknown fibril parameters $\left( {\mu_{o} ,I_{o} ,a_{1} ,a_{2} ,a_{3} } \right)$ with each parameter confined by reasonable bounds as provided in earlier studies (Adouni and Dhaher 2016; Tang et al. 2009). To ensure unbiased estimate of the fibril parameters $\left( {\mu_{o} ,I_{o} ,a_{1} ,a_{2} ,a_{3} } \right)$, multiple sets of the five fibril parameters selected from a plausible range of values were then used (MultiStart MATLAB function) as the initial inputs to the optimization procedure. The multiple outcomes of the optimization process were then averaged to represent the best fit of the fibril parameters used to characterize the fibril continuum model employed in the subsequent FEA simulations.

Figure 8a displays the fibril stress–strain plot of MD data (solid line) and its corresponding fitted curve (dotted line) for the TC crosslink failure (β), varying from native fibril (100% crosslink) to no crosslink. Figure 8b shows the stress–strain curves (MD data and fitted curve) of a fibril at different levels of surface degradation (δ) varying from an intact fibril.

Fig. 8

Curve fitting by nonlinear optimization to MD simulation data to estimate input fibril parameters for a tropocollagen crosslink failure and b generalized surface degradation

For the fitted data, the goodness-of-fit (GOF) values—coefficient of determination, R² (presented as mean ± standard error (SE))—were found 0.983 ± 0.04 for the crosslink failure and 0.978 ± 0.02 for the surface degradation. The fitted curves exhibit acceptable fits to the MDS data for the two degeneration mechanisms.

Appendix C

3.1 Estimation of surface degradation

The enzymatic surface degradation can be estimated based on the fraction of surface occupied by the adsorbate (enzyme). The accumulation of the adsorbate on a surface (adsorbent or substrate) is therefore known as adsorption. In the experiment carried out by Laasanen et al. (2003), collagenase type VII (C 0773, Sigma Chemical Co., St. Louis, MO, USA) was applied to degrade the surface of collagen-II in which surface degradation was estimated considering the binding of collagenase molecules on the collagen surface (Fig. 9). The incubation time (in 37 °C, 5% CO₂ atmosphere) for collagenase (30 U/ml)-treated samples was 44 h (Laasanen et al. 2003). Because of unavailability of specific parametric values, the surface coverage (degradation) was calculated using Langmuir adsorption isotherm for a range of physical and mechanical properties widely used for collagen and collagenase enzyme. The Langmuir adsorption isotherm provides one of the simplest and most direct methods to quantify the adsorption process. Since the Langmuir isotherm model typically well suited with isotherm data from protein/enzyme adsorption studies, the model is often used to estimate the protein binding affinity (Metzmacher et al. 2007a, b; Wilson et al. 2005).

Fig. 9

Schematic of the estimation of surface degradation

3.2 Surface coverage

Surface coverage as per Langmuir Isotherm equation is

$$\theta = \frac{{k_{\text{a}} C}}{{1 + k_{\text{a}} C}}$$

(12)

where θ is the fractional surface coverage, C is the molar concentration of the solution (here collagenase VII), $k_{\text{a}}$ is the association constant and $k_{\text{d}}$ is the dissociation or binding constant. $k_{\text{a}}$ and $k_{\text{d}}$ are related as $k_{\text{a}} = \frac{1}{{k_{\text{d}} }}$.

The Collagen Digestion Unit (CDU) of the collagenase varies between 1000 and 3000 CDU/mg, and its molecular weight is between 68 and 125 kDa (Sigma-Aldrich; Webb 1992). For the ranges of CDU and molecular weight, the concentration, C, of collagenase (30 U/ml) was estimated in the range of $C = 0.08 \times 10^{ - 6} \sim 0.44 \times 10^{ - 6} \,{\text{mol}}/{\text{L}}$.

The binding constant, $k_{\text{d}}$, between collagenase VII and collagen varies widely (O(10⁻⁵) to O(10⁻⁷)) whether determined experimentally and computationally (Addi et al. 2016; Bella et al. 1994; Matsushita et al. 1998, 2001; Toyoshima et al. 2001). However, in most of the reported literature, $k_{\text{d}}$ was estimated to be in the order of O(10⁻⁶) (Bella et al. 1994; Evans 1981; Matsushita et al. 1998, 2001; Toyoshima et al. 2001; Wilson et al. 2003) and was used to calculate the surface coverage using Eq. (12). Considering $k_{\text{d}} = 4.99 \times 10^{ - 6}$ (Wilson et al. 2003), the fraction of binding sites on collagen occupied by collagenase was estimated to $\theta_{ \hbox{min} } \cong 0.016$ and $\theta_{\hbox{max} } \cong 0.081$.

3.3 Estimation of collagen fibrils in cartilage

The number of collagen fibril needs to estimate for the cartilage tissue used in the experiment (Laasanen et al. 2003). The following assumptions are rationally considered to calculate the amount of collagen fibrils in cartilage tissue:

• 1. Superficial layer or zone is typically 15% of articular cartilage tissue height $\left( {h_{\text{tissue}} } \right)$(Adouni et al. 2012; Hollander et al. 1994; Shirazi et al. 2008; Shirazi and Shirazi-Adl 2008).

• 2. Amount of collagen fibers in the superficial layer is approximately 15% (Adouni et al. 2012; Shirazi and Shirazi-Adl 2009; Shirazi et al. 2008).

• 3. Fibers are idealized as perfect cylinders of tightly packed monomeric fibrils

• 4. The fiber diameter is larger than the fibril diameter, $d_{\text{fb}} \gg d_{\text{fl}}$.

• 5. Diameter of a fiber $\left( {d_{\text{fb}} } \right)$ in cartilage generally in the range of $15.2 \pm 8.3 \sim 29.2 \pm 5.4\,{\text{nm}}$ (Halberg et al. 1988; Holmes and Kadler 2006; Moskowitz 2007; Mwenifumbo et al. 2007; Silver and Siperko 2003; Watanabe et al. 1994)

• 6. Length of the fiber is assumed to 1 µm (Domene et al. 2016; Gautieri et al. 2011; Liu et al. 2015)

An idealized fiber of tightly packed fibrils is shown in Fig. 10. This idealization is valid as long as the fiber diameter, $d_{\text{fb}}$, is much larger than the diameter of the fibril, $d_{\text{fl}}$. When the collagen fibrils come in contact with the enzyme solution, the enzyme diffuses into the fibers and binds to specific sites on fibrils located at the surface of the fibers. Due to the size and organization of the fibrils, the enzyme molecules cannot penetrate the tightly packed (crosslinked) fibrils that make up an individual fiber. The fibrils are also assumed to be incompressible and nonoverlapping (small circles in Fig. 10b).

$${\text{Volume}}\,{\text{of}}\,{\text{sample}}\,{\text{cartilage}}\,{\text{plug}},\quad V_{\text{tissue}} = \left. {\pi r^{2} h} \right|_{\text{tissue}} \approx 17.6\,{\text{mm}}^{3}$$

$${\text{Volume}}\,{\text{of}}\,{\text{the}}\,{\text{fibers}}\,{\text{in}}\,{\text{the}}\,{\text{superficial}}\,{\text{layer}}\,{\text{of}}\,{\text{the}}\,{\text{cartilage}}\,{\text{plug}} \approx 0.4\,{\text{mm}}^{3}$$

Fig. 10

Schematics of a fiber (a) and its cross section (b) perpendicular to the axis of the fiber [based on (Tzafriri et al. 2002)]

Total number of collagen fibers is assumed to be $n_{\text{fb}}$ that contribute to the total fiber volume. However, the fiber dimeter is variable depending on its young or mature state, and consequently, the diameters of fibers in the superficial layer vary widely (Moskowitz 2007; Mwenifumbo et al. 2007; Silver and Siperko 2003). Taking into account the minimum and maximum mean values of the fiber diameter, the minimum fiber volume, $V_{\text{fb}}^{\hbox{min} }$, and maximum fiber volume, $V_{\text{fb}}^{\hbox{max} }$, of a single collagen fiber have been estimated. Thus, based on the $V_{fb}^{\hbox{min} }$ and $V_{\text{fb}}^{\hbox{max} }$, the minimum and maximum number of collagen fibers was estimated to $n_{\text{fb}}^{\hbox{min} } = 5.913 \times 10^{11}$ and $n_{\text{fb}}^{\hbox{max} } = 2.216 \times 10^{12}$, respectively.

3.4 Estimation of tropocollagen molecules on the surface

During the surface degradation process, enzymes bind and cleave the tropocollagen molecules on the surface (along the circumference of the fibril) and are schematically shown in Fig. 11. The TC molecules were assumed to be arranged tightly around the circumference of the fibrils. The dimensions of each TC molecule are $d_{\text{TC}} = 1.5\,{\text{nm}}$, $l_{\text{TC}} = 300\,{\text{nm}}$, with a (D-periodic) gap of $67\,{\text{nm}}$ between the TC molecules (Chen et al. 1995; Franchi et al. 2007; Graham et al. 2004; Shoulders and Raines 2009).

Fig. 11

Schematic of a fiber cross section perpendicular to the axis of the fibril showing TC molecules on the surface

Considering the minimum and maximum circumference of the fibers, staggered organization of TC molecules with the periodicity of $67\,{\text{nm}}$ and (collagen type II) fiber length of $1\,\upmu{\text{m}}$ (Domene et al. 2016; Gautieri et al. 2011; Liu et al. 2015), the minimum and maximum number of TC molecules around the circumferential was estimated to $n_{{{\text{TC}}_{\text{surf}} }}^{\hbox{min} } = 5.68 \times 10^{13}$ and $n_{{{\text{TC}}_{\text{surf}} }}^{\hbox{max} } = 3.89 \times 10^{14}$.

3.5 Estimation of surface degradation

The collagenase type VII used in the Laasanen et al. experiments is bacterial collagenase, which acts both as a collagenase and a gelatinase. Collagenases can only cleave at one site on collagen molecules, but gelatinases can cleave the rest of the sites on the collagen once collagenase has made an initial cut. The bacterial collagenase can cut in multiple places on collagen and was accounted here for estimating the number of binding sites (Fig. 12). Thus, in the current estimation, we considered up to three binding sites/TC (Piluso et al. 2017; Riley and Herman 2005; Xu et al. 2000; Zderic 1995). Hence, the total number of collagenases on cartilage surface was estimated in the range of $n_{\text{Collagenase}}^{\text{surface}} \approx 2.726 \times 10^{12} \sim 9.453 \times 10^{13}$ molecules.

Fig. 12

Surface degradation due the cleavage of molecules by enzyme

To the best of our knowledge (with extensive literature search), two different rates of collagen cleavage by the bacterial collagenase have been found and considered here in the estimation process. These are: (i) 22 molecules of collagen are degraded per molecule of collagenase per hour (Barrett et al. 2012; Welgus et al. 1980) and (ii) a unit of collagenase activity is defined to produce 10% cleavage of collagen in 2.5 h at 37^o C (Zderic 1995).

Considering cleavage rate as per (i), the number of degraded molecules for 44 h of incubation time (Laasanen et al. 2003) was in the range of $n_{\text{collagen}}^{{{\text{degraded}}\,{\text{mol}}}} \approx 4.67 \times 10^{16} \sim 2.565 \times 10^{17}$ molecules. With this, the minimum and maximum surface degradation was estimated to 0.01% and 0.201%, respectively. According to cleavage rate as per (ii), the number of degraded molecules for the same incubation time (Laasanen et al. 2003) was in the range of $n_{\text{collagen}}^{{{\text{degraded}}\,{\text{mol}}}} \approx 3.0 \times 10^{15} \sim 2.054 \times 10^{16}$ molecules. With this, the minimum and maximum surface degradation was estimated to be 0.013% and 3.51%, respectively.

Although Langmuir adsorption isotherm is often considered to estimate the protein binding affinity, the Langmuir model does not take into account the interactions between enzymes—especially when the enzymes (considered as voluminous objects) are charged exhibiting enzyme-to-enzyme interactions and, in a context of a joint, exhibiting interactions with other molecules in the synovial fluid. Thus, the Langmuir model neglects any type of chemical interaction with the surroundings except the surface. Moreover, a fixed number of adsorption sites are typically approximated on a static surface under equilibrium conditions as per the Langmuir model. But, experimentally the surface of a collagen fibril is dynamic in nature since the number of adsorption sites changes with time as the degradation progresses, and thus, equilibrium is difficult to achieve. The Langmuir model, however, phenomenologically seems to capture the behavior of surface degradation and the associated adsorption sites that not only give a gross estimation of the amount of surface degradation but provide an insight into the degradation mechanisms as well.