Effect of Degradation and Osteoarthritis on the Viscoelastic Properties of Human Knee Articular Cartilage: An Experimental Study and Constitutive Modeling

: Articular cartilage, as a hydrated soft tissue which covers diarthrodial joints, has a pivotal role in the musculoskeletal system. Osteoarthritis is the most common degenerative disease that affects most individuals over the age of 55. This disease affects the elasticity, lubrication mechanism, damping function, and energy absorption capability of articular cartilage. In order to investigate the effect of osteoarthritis on the performance of articular cartilage, the mechanical behavior of human knee articular cartilage was experimentally investigated. Progressive cyclic deformation was applied beyond the physiological range to facilitate degradation of the tissue. The relaxation response of the damaged tissue was modeled by means of a fractional-order nonlinear viscoelastic model in the framework of ﬁnite deformations. It is shown that the proposed fractional model well reproduces the tissue’s mechanical behavior using a low number of parameters. Alteration of the model parameters was also investigated throughout the progression of tissue damage. This helps predict the mechanical behavior of the osteoarthritic tissue based on the level of previous damage. It is concluded that, with progression of osteoarthritis, the articular cartilage loses its viscoelastic properties such as damping and energy absorption capacity. This is also accompanied by a loss of stiffness which deteriorates more rapidly than viscosity does throughout the evolution of tissue damage. These results are thought to be signiﬁcant in better understanding the degradation of articular cartilage and the progression of OA, as well as in the design of artiﬁcial articular cartilages.


Introduction
Articular cartilage (AC), as a hydrated soft tissue which covers diarthrodial joints, has a pivotal role in the musculoskeletal system. It reduces friction at the surface of joints, minimizes contact stresses by distributing loads [1,2] and facilitates transmission of cyclic loads without damage under normal physiological loading. AC is actually a matrix, saturated with water and mobile ions. The matrix contains chondrocytes embedded in the extracellular matrix (ECM). The mechanical properties of AC mostly depend on the ECM, which consists of collagen fibers (primarily type II collagen) and negatively charged proteoglycans (PGs)-responsible for tolerating tensile and compressive loadings, respectively [3,4]-while chondrocytes are responsible for ECM maintenance.
Osteoarthritis (OA) is the most common degenerative disease in individuals over the age of 55 [5,6]. OA is a multifactorial disease; however, its progression is described by depletion of PGs and break down of collagen fibrils, leading to changes in the mechanical properties of the AC. OA affects the elasticity, lubrication mechanism, damping function, and energy absorption capability of the tissue. Osteoarthritic damage is a chemomechanical process where the mechanical properties of the tissue strongly depend on the OA grade [7,8]. Nonetheless, as there is no effective treatment for OA, it is highly important to investigate alteration of the main constituents of the tissue throughout the progression of disease in order to estimate their functionality at various OA grades. This can also be helpful in the development of artificial cartilages.
Many studies have been conducted in order to develop constitutive relations for the mechanical behavior of articular cartilage. The prevailing theories account for the biphasic (solid-fluid) and triphasic (solid-fluid-ionic) structure of the cartilage. Mow et al. demonstrated the ionic nature of the cartilage where synovial fluid contributes to the mechanical response by being drawn into the porous matrix [9]. Other research has focused on interaction between the collagen matrix and the lubricating synovial fluid which permeates the joint capsule. This interaction includes frictional drag and compressibility of the matrix. For more details, one can refer to [10][11][12], among others. Though these models are physiologically comprehensive, they mostly do not fit the experimental results well [13][14][15].
AC exhibits a time-dependent mechanical response, meaning that reaching equilibrium under an external loading may take from a few minutes to hours, depending on the mode of loading [13]. This time-dependent behavior, also called viscoelasticity, is thought to be caused by both the fluid content and solid matrix of the tissue [16]. Stress relaxation is one of the important characteristics of viscoelastic materials. A variety of experimental methods are available to quantify viscoelastic properties. Moreover, viscoelastic mechanical models can be used to describe the viscoelastic behavior (e.g., stress relaxation and creep) of the tissue. Most viscoelastic models make use of spring-dashpot combinations to represent the time-dependent behavior of AC. As AC exhibits a nonlinear mechanical behavior, such models rarely provide a good fit to experimental data unless they include a few springs and dashpots [17,18]. Smyth et al. employed fractional calculus to model the stress-relaxation of AC in small strain regime [19]. They made advantage of fractional calculus to reduce the number of required model parameters.
Researchers have modeled the fibril network alone by means of spring-dashpot combinations [20][21][22]. Wilson et al. proposed a poroviscoelastic fibril-reinforced model which focuses on the morphology of the collagen fibers and their influence on stress and strain distribution [21]. Later, Julkunen et al. showed a profound agreement between the experimental data and Wilson's model in stress-relaxation [23]. Stender et al. proposed an anisotropic constitutive model to account for the development of damage in the AC [24]. They did not consider the viscoelastic behavior as well as the effect of fluid flow in the AC. Hosseini et al. focused on a fibril-reinforced viscoelastic biphasic model to predict the progression of induced damage in AC, as well as the relationship between damage to fibrils and PGs [25]. Korhonen et al. investigated the effects of collagen network degradation, as well as PG depletion, on the stress-relaxation response by applying a poroelastic model [26]. Experimental investigations demonstrate that excessive mechanical loading can lead to degradation of ECM [27][28][29]. Khajehsaeid and Abdollahpour investigated deformationinduced degradation of human knee articular cartilage by applying excessive deformations experimentally. They studied the evolution of the mechanical properties throughout the progression of tissue damage and concluded that modeling the alterations of mechanical properties can be an effective method for predicting degradation of AC under excessive loading [8]. Recent works have proposed degeneration algorithms to study the effect of OA-induced changes on the mechanical behavior of AC [30,31]. They simulate collagen degradation by decreasing the stiffness of the collagen network.
In this work, the stress-relaxation response of human knee articular cartilage was experimentally investigated. To facilitate degradation of the articular cartilage, progressive cyclic deformation was extended beyond the physiological range. The nonlinear relaxation response of the damaged tissue was modelled in the framework of finite deformations. To reduce the number of required material parameters, a fractional-order nonlinear viscoelastic model is proposed. Alteration of the model parameters was also investigated throughout the progression of damage in order to predict the mechanical behavior of the degraded tissue based on the level of previous damage.

Experiments
To conduct experimental observations on the viscoelastic behavior of osteoarthritic cartilage, 14 samples were harvested from removed knee cartilages of 11 patients undergoing total knee arthroplasty; 6 males (60 ± 3 years old) and 5 females (60 ± 2 years old). Specimens with dimensions 30 × 10 × (3 ± 1) mm (as shown in Figure 1) were cut from parts of the femoral cartilage that appeared almost healthy to the orthopedic surgeons. The specimens were stored in phosphate buffered saline (PBS) and processed penicillin at 4 • C. As tests were performed just 15 h after surgery, freezing was not required. Prior to testing, the specimens were kept at 20 • C for 30 min to reach swelling equilibrium. As excessive mechanical loading can lead to breakdown of the cartilage network and depletion of the PGs, we aimed to emulate OA-associated damage via excessive loading of the tissue. It has been demonstrated that alteration of the tissue's constituents may lead to OA, which can be reflected in the mechanical properties of the tissue. Most studies focus on compressive tests of AC; however, in this study, to achieve larger deformations and accumulated damage in testing, degradation of cartilage was investigated through cyclic tensile tests. This provides an opportunity to study the alteration of the viscoelastic properties throughout OA/degradation progression. Therefore, the specimens were subjected to progressive cyclic uniaxial tensile tests with constant speed of 5 mm/min, resulting in a constant strain rate of 8.3 × 10 −3 s −1 .
The tests were conducted using an STM-2 (SANTAM) Servo Electromechanical Universal Tensile Testing Machine. Upon completion of the loading step, the specimens were allowed 5 min to relax prior to unloading. After each load-hold-unload cycle, each specimen was kept in PBS for 20 min to fully recover prior to the next cycle.

Modeling
In order to study the stress-relaxation response of AC and alteration of its viscoelastic properties, one needs to model the tissue behavior by means of viscoelastic models. The standard linear viscoelastic models include either parallel combination of Maxwell elements or serial combination of Kelvin elements. For example, for a model consisting of n Maxwell elements, as shown in Figure 2, we have: where σ is the total stress applied to the model and ε i e and ε i v are elastic and viscous strains in the i-th element, respectively. In (1), τ ι = η ι /E ι is called the relaxation time for the i-th element, and E 0 stands for the instantaneous stiffness of the model (i.e., the stiffness at extremely high loading rate): The above model exhibits a solid-like behavior if at least one Maxwell element does not include any dashpot, but only a spring with stiffness E ∞ . To model the viscoelastic behavior of AC, one needs to employ at least 3 Maxwell elements needing at least 5 model parameters (3 coefficients and 2 relaxation times associated with the Prony series) to be determined from experimental tests [32]. This number of parameters often cannot be uniquely identified by only one experiment. Therefore, it is advantageous to reduce the number of model parameters as much as possible, provided that the model is still able to reproduce the material behavior accurately enough.
In this study, a fractional-order nonlinear viscoelastic model is proposed to characterize the stress-relaxation behavior of AC by a reduced number of model parameters. The fractional-order viscoelastic models employ fractional dashpots instead of the ordinary dashpots of classic viscoelastic models.
As shown in Figure 3, considering a single chain in (1) and replacing ε v by ε − ε e one obtains a constitutive equation of the standard Maxwell element: In order to obtain a fractional Maxwell element, a fractional dashpot replaces the ordinary one. The differential equation of the fractional dashpot is written as follows: where η is the viscosity associated with the fractional dashpot. The obtained fractional Maxwell element shown in Figure 3 is characterized by the following constitutive equations: In (4) and (5), d α dt α corresponds to a fractional derivative of order α which is defined according to the Reimann-Liouville definition: where Γ is the Gamma function: The stress-relaxation response of the fractional Maxwell element for a unit step deformation ε(t) = u(t) is defined as [33]: where the kernel function E α is the α-order Mittag-Leffler function: As discussed earlier, in order to yield a solid-like behavior, we need to add at least one spring in parallel, as shown in Figure 4. The constitutive equation of the obtained model is written as follows [34]: This approach can be extended to finite deformations as well. Following the method proposed in [35], (1) can be rewritten in the framework of finite deformations as follows: where W i are the strain energy functions associated with the nonlinear springs, σ and ε refer to any conjugate finite stress and strain tensors, Q 1 is a stress-like internal variable, and W 0 = W 1 + W ∞ is the instantaneous stored strain energy. Following a similar approach for the fractional model presented in (10), and using the Cauchy stress σ and the deforma-tion gradient F as conjugate tensors, one obtains the model shown in Figure 5, with the constitutive Equation (12). In modeling the loading and unloading behavior of a viscoelastic material, it is important to choose an appropriate strain energy function capable of reproducing the nonlinear mechanical behavior of the material under study [36][37][38]. As most strain energy functions are defined in terms of invariants of the left Cauchy-Green deformation tensor B = FF T , the differentiations in (12) should be written with respect to B or its invariants equivalently [39,40]: In (13), I k are invariants of B. Separating the effects of deformation and time is a common practice in the theory of viscoelasticity [41,42]. Thus, the general solution of (11)-(13) can be written as follows: where σ ins (ε) is the contribution of the deformation and h(t) tands for the effect of deformation rate (time). Solution of the fractional constitutive equation represented in (10) may not be straightforward for an arbitrary deformation history. It can be analytically determined only for simple deformation histories. For example, the stress-relaxation solution of (10) for a unit step deformation was obtained in [43]: Specializing the above equation for α = 1 one obtains the expression for the stressrelaxation function corresponding to the standard linear solid, noting that Γ(1 + k) = k!: For an arbitrary deformation, the deformation history ε(t) must be introduced to the following relation [43]: where s is the Laplace variable and L −1 is the inverse Laplace transform which can be calculated as: In Mittag-Leffler notation we can write: The same procedure can also be adopted to (12) for finite deformations. In case of finite deformations, one may choose the neo-Hookean energy W = 1 2 µ(I 1 − 3) because of its simplicity to solve for the Cauchy stress in a step deformation, where µ is the shear modulus of the material. Using (13) and (15): where µ 0 and µ ∞ are the instantaneous and the asymptotic shear moduli of the model, respectively. For an isotropic incompressible material (J = 1), subject to uniaxial tension, B can be written as follows: In (21), the relation λ 1 = λ, λ 2 = λ 3 = λ 1/2 was adopted. Due to the stress-free conditions in the transverse directions (σ 22 = σ 33 = 0), the first component of the Cauchy stress is calculated by means of (20) and (21): For a ramp deformation (i.e., constant strain rate), the relation ε(t) = .
εt should be used in conjunction with (17) to determine the corresponding stress response in terms of time. However, if one is interested only in the relaxation stage, the normalized stress σ(t)/σ max may be easier to investigate. According to (14), the normalized relaxation stress will be a function of time only as the deformation is kept constant during the relaxation stage. So we have: (23) In (23), τ, α, and µ ∞ /µ 0 are the parameters associated with the nonlinear model presented in Figure 5. The model tries to reproduce the stress-relaxation data by a minimum number of parameters which can be determined by means of available experimental data.

Results
Deformation histories of the specimens are shown for the female and male specimens in Figure 6. At strain amplitudes of less than (average) 12.5% for female specimens and 15% for male specimens, no significant change was observed in the stress-strain behavior, revealing no permanent damage in the tissue. In order to study the impact of degradation on the viscoelastic properties of AC, only the cycles with strain amplitude more than 12.5% for females and 15% for males are reported in this study. In this study, we focus on the stress-relaxation response of AC and the alteration of its viscoelastic properties. The (average) stress-relaxation results are shown in Figure 7. In order to focus on the early stage of relaxation, only the initial 200s of the relaxation tests are shown. The stress-relaxation response of the AC tissue was reproduced by means of the fractional-order viscoelastic model presented in (13), (20), and (23); the results are shown in Figures 8 and 9 for the female and male specimens, respectively.
The optimization toolbox of MATLAB was utilized to identify the model parameters by minimizing an appropriate cost function. The obtained parameters are reported in Tables 1 and 2 for the female and male specimens, respectively.  The stress-relaxation response of the AC tissue was reproduced by means of the fractional-order viscoelastic model presented in (13), (20), and (23); the results are shown in Figures 8 and 9 for the female and male specimens, respectively.  The optimization toolbox of MATLAB was utilized to identify the model parameters by minimizing an appropriate cost function. The obtained parameters are reported in Tables 1 and 2 for the female and male specimens, respectively. Using the proposed fractional-order nonlinear viscoelastic model, the material parameters were identified. Figure 10 show alteration of the parameters with the maximum applied strain amplitude, where it was observed that the evolution of parameters was almost linear through the degradation process. The OA grade regions in Figure 10 are defined based on the amount of maximum I 1 experienced previously by the tissue [8]. Choosing another measure for the deformation threshold might shift the areas slightly. Two predicted curves are also provided in Figures 8 and 9 for the tissue with prior 35% and 45% deformation, respectively. It should be noted that, within the physiological range (shown as the safe region in the figure), the material parameters are unaltered as the tissue does not experience unrecoverable microstructural changes or permanent alteration in its mechanical properties. In this study, we focused on the stress-relaxation response of AC and the alteration of its viscoelastic properties, while the stress-strain response of osteoarthritic AC was thoroughly studied in an authors' recent work where tissue degradation, as well as alteration of the mechanical properties, were linked to OA [8].

Discussion
It was observed that the proposed fractional-order nonlinear viscoelastic model well reproduces the AC stress-relaxation response by means of only three parameters. Alteration of the parameters with the maximum applied strain amplitude was also investigated. Beyond the physiological threshold, the excessive mechanical loading leads to unrecoverable degradation of the tissue [8]. The tissue degradation affects its mechanical properties which can be related to the grade of OA [44]. It should be noted that the physiological deformation threshold is an average value for the examined samples which not only depends on the sex, age, ethnicity, and location of the specimen, but may also depend on the deformation rate [45].
The results shown in Figure 10 help predict behavior of the OA tissue, meaning that once the physiological threshold is exceeded, the behavior of a degraded tissue is predictable based on the maximum applied previous deformation. This is thoroughly discussed in [8] where a relation was proposed between the maximum deformation of the tissue and its damage. This allows prediction of the amount of damage and grade of OA associated with a mechanical loading beyond the physiological range. Figure 10 also shows that the relaxation time (τ) decreases with increase in the maximum deformation ever experienced, which is in fact related to damage of the tissue. This can be explained by a significant decrease in the PG content (decrease in the number of PG-collagen bonds) and the consequent availability of more space in the ECM for fluids to flow out of matrix more easily and rapidly, resulting in a lower relaxation time [46]. In addition, some studies have connected the reduction in the relaxation time to the destruction of the collagen network according to the theory of polymer dynamics. June et al. found that, after enzymatic digestion by bacterial collagenase and hyaluronidase, changes in the stress relaxation of AC become faster [47,48]. Thus, we can conclude that, through the progression of OA, the AC tissue relaxes more rapidly after a loading cycle, and this is assumed to be a sign of loss of viscoelastic properties such as damping and energy absorption capacity.
In general, the viscoelastic behavior of AC is associated with a flow-independent mechanism, e.g., resistance of the collagen fibrils in tension and PGs in compression, as well as a flow-dependent mechanism, e.g., the frictional drag force of interstitial fluid flow. Elasticity of the tissue mostly comes from resistance of the collagen fibrils, as well as PGs. This is why, for a degraded tissue, the mechanical behavior approaches a fluid-like behavior, which is exhibited in the increase in α [49][50][51][52]. The loss of stiffness is also reflected in the decrease in µ ∞ /µ 0 . This reveals that the tissue elasticity deteriorates more rapidly than its viscosity does throughout the damage evolution/OA progression.
The present study only considers low loading rates. The rate of loading may also affect the results. This is assumed to not only alter the physiological threshold but also the microstructural mechanisms of degradation. This could be the subject of a future study.

Conclusions
In this work, the stress-relaxation behavior of human knee articular cartilage was experimentally investigated through progressive cyclic deformation beyond the physiological range to facilitate degradation of the tissue. A fractional-order nonlinear viscoelastic model was proposed for the damaged tissue in the framework of finite deformations. The proposed model well reproduced the tissue's stress-relaxation behavior using a low number of parameters. Alteration of the model parameters was also investigated which helps predict behavior of the osteoarthritic tissue. This means that once the physiological threshold is exceeded, the mechanical behavior of the degraded tissue is predictable based on the level of previous damage. It was concluded that, with progression of OA, the AC tissue loses its viscoelastic properties such as damping and energy absorption capacity. This is also accompanied by a loss of stiffness, which deteriorates more rapidly than viscosity does throughout the damage evolution. These results are thought to be significant in better understanding the degradation of articular cartilage and the progression of osteoarthritis, as well as in the design of artificial articular cartilages  Informed Consent Statement: Informed consent was obtained from all subjects involved in the study.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author.

Conflicts of Interest:
The authors declare no conflict of interest.