|Home | About | Journals | Submit | Contact Us | Français|
Articular cartilage is the load bearing soft tissue that covers the contacting surfaces of long bones in articulating joints. Healthy cartilage allows for smooth joint motion, while damaged cartilage prohibits normal function in debilitating joint diseases such as osteoarthritis. Knowledge of cartilage mechanical function through the progression of osteoarthritis, and in response to innovative regeneration treatments, requires a comprehensive understanding of the molecular nature of interacting extracellular matrix constituents and interstitial fluid. The objectives of this study were therefore to (1) examine the timescale of cartilage stress-relaxation using different mechanistic models and (2) develop and apply a novel (termed “sticky”) polymer mechanics model to cartilage stress-relaxation based on temporary binding of constituent macromolecules. Using data from calf cartilage samples, we found that different models captured distinct timescales of cartilage stress-relaxation: monodisperse polymer reptation best described the first second of relaxation, sticky polymer mechanics best described data from ~1-100 seconds of relaxation, and a model of inviscid fluid flow through a porous elastic matrix best described data from 100 seconds to equilibrium. Further support for the sticky polymer model was observed using experimental data where cartilage stress-relaxation was measured in either low or high salt concentration. These data suggest that a complete understanding of cartilage mechanics, especially in the short time scales immediately following loading, requires appreciation of both fluid flow and the polymeric behavior of the extracellular matrix.
Articular cartilage is a complex macromolecular biopolymer [1-4]. The tissue constituents include type II collagen, the aggregating proteoglycan (aggrecan), and numerous other molecules such as decorin, superficial zone protein, and numerous collagens. Further, the hydrated tissue has a distinct zonal structure with defined orientations of collagen and precise localization of proteoglycans. Macroscopically, the cartilage ultrastructure gives rise to a unique deformable and load-bearing solid with low friction and wear characteristics at the articular surface. Unfortunately, disruption of the tissue structure during the course of degenerative diseases, such as osteoarthritis, results in altered mechanics and severe tissue wear.
Historically, cartilage viscoelasticity has been separated into flow-dependent and flow-independent contributions [5, 6]. The flow-dependent contribution results from fluid flowing through a permeable solid matrix and occurs on a long timescale as a result of the low permeability of cartilage . The flow-independent portion of cartilage relaxation includes viscoelasticity not resulting from fluid flow and is often modeled phenomenologically using Fung's quasi-linear viscoelastic model [8, 9]. Previously, we demonstrated that the relaxation function derived for simple (monodisperse) polymeric solutions (i.e. reptation), where all polymers are assumed to have the same molecular weight, was able to model the early stress relaxation of articular cartilage in unconfined compression . Although a generally-accepted mechanistic explanation of flow-independent viscoelasticity remains to be determined, this work remains  as one candidate mechanism.
While our previous work modeled cartilage stress relaxation by monodisperse reptation, it is well-known that cartilage macromolecules are polydisperse (e.g. ) and likely exhibit complex molecular interactions. Considering the complex nature of cartilage structure, it is possible that a polydisperse polymeric solution, where the polymers are assumed to have a distribution of molecular weights, and thus a more complex set of interactions, provides a better description of relaxation data. Additionally, the incorporation of fluid flow in the numerical modeling is essential considering the flow-dependent contributions to viscoelasticity for the hydrated tissue. However the extent to which any of these mechanistic models describes the time course of relaxation is unknown. Thus, our knowledge of cartilage viscoelasticity, especially for timescales immediately following loading, is largely incomplete.
The purpose of this paper is to investigate the relative contributions of monodisperse and polydisperse polymeric models, as well as fluid flow, to cartilage viscoelasticity. Herein, we first compare monodisperse and polydisperse models for short-term stress relaxation, then combine the polymer model with the well known KLM fluid flow model [5, 12] to separate relaxation caused by flow-dependent and flow-independent mechanisms. The objectives of this study were to (1) examine the timescale of cartilage stress-relaxation using different mechanistic models and (2) develop and apply a novel (termed “sticky”) polymer mechanics model, based on temporary binding of constituent macromolecules, to cartilage stress-relaxation. We further tested the sticky polymer model using experimental measurements of cartilage stress relaxation at both low and high salt concentrations, under the hypothesis that high salt concentration would screen temporary electrostatic bonds.
Many biological molecules are long relative to their diameters. A large length-to-diameter ratio is characteristic of structural molecules in connective tissues  and many engineering polymers [14-17]. Mechanical properties in a monodisperse system, or one where all polymers are assumed to have the same molecular weight, are governed by molecules that are long, flexible, and continuously changing conformation . Reptation describes this constant motion as “every [polymer] chain, at a given instant, is confined within a ‘tube’ as it cannot intersect the neighboring chains. The chain thus moves inside the tube like a snake.”  (Figure 1). Stress relaxation based on reptation theory assumes that the stress relaxes with the movement of the polymer chains in their tubes . After creation of the tube, the chain moves from the original tube by diffusion. For a stress relaxation experiment it is assumed that the stress is proportional to the fraction of the chain remaining in the original tube. For a monodisperse, reptating system, the fraction of the chain remaining in the initial tube (equivalent to the relaxation function) is
where τd is the characteristic disengagement time for the chain to escape from its tube, p is an odd counting variable and t represents time. The theoretical development of the relaxation function provides an estimate of the disengagement time τd in terms of measures of the molecular structure. However this is not relevant to the current study and was discussed previously .
Cartilage is a heterogeneous material made of multiple polydisperse polymers with various binding interactions [19-21]. Subsequently, we demonstrate how three additional complexities in a polymer model separately result in a stretched exponential stress relaxation function: . The stretching parameter (β) depends on the polymer mechanism being modeled and may provide insight into active polymer mechanisms in articular cartilage.
The polymer relaxation function (Equation 1) was derived assuming that all of the molecules are linear and of the same length (i.e. a monodisperse system). In a natural system such as cartilage there are different-sized biopolymers interacting to result in macroscale mechanical properties . The issue of a polydisperse reptating polymer system was examined by de Gennes who determined the relaxation function for a system with an exponential molecular weight distribution 
where, f is an average relaxation rate, t is time, and β=x/(x+1) where x>0 is the weight distribution exponent. De Gennes' derivation was for broad classes of reptating systems and demonstrated that for exponential weight distributions that a stretched exponential (or Kolrausch-Williams-Watt, abbreviated KWW) relaxation function is expected for temperatures above the glass transition temperature. In the original derivation, the exponent β was predicted to be temperature-independent with a range between 0.25 and 0.33 for the specific weight distributions considered.
Stretched exponential behavior also can arise from mechanisms other than reptation of a polydisperse system. Another source of stretched exponential behavior occurs when there are barriers to molecular motion (e.g. steric interference or temporary bonding). For example, Edwards and Vilgis  examined the case of long rod molecules that interfere with each others' motions and demonstrated that stress relaxation behaved as , where σ0 is the initial stress and τ is a time constant. This is a less-general stretched exponential expression than derived by de Gennes, but has the same form with β=0.5.
Motivated by the work of Edwards and Vilgis, we examined the case of a linear polymer having NT temporary binding sites with probability p of being unbound at any particular time and spaced such that the length of the free segment is Lm when m consecutive sites are unbound. This theory was developed because of both the well-known non-covalent interactions of cartilage matrix molecules  and the physical properties of cartilage proteoglycans, which were originally defined as “mucopolysaccharides”  because of their high viscosity. The expected numbers Nm of free lengths of polymer of length Lm are (see Electronic Supplement, Appendix 1)
For NT of a reasonable size, there will always be at least some bound sites, and, therefore, relaxation will not be by simple reptation, depending upon the probability of binding, p. However, each free section of the polymer is mobile, and relaxation of stress may occur by motion of the free sections. The characteristic relaxation function for a free segment will be assumed as: Ψ(t,Lm), and thus the overall relaxation function for the polymer chain is formed as a weighted sum over the free lengths Lm as
In Equation 3, the number of sites of length Lm is seen to be proportional to the probability of the sequence “0111…1110” with “m” 1's as would be expected. In Eq. 4 the relaxation function in the sum is weighted by the probability of the sequence “111…111” with “m” 1's. Assuming Debye (single exponential) relaxation of stress for each free segment m with a time constant proportional to Lm as: Ψ(t,Lm) = exp (−at/Lm), with constant a, results in the following “sticky” relaxation kernel:
This can be fit well by a stretched exponential function with 0.58<β<0.69 for many combinations of parameters (Electronic Supplement, Appendix 2).
In summary, the different complex polymer models result in distinct predictions for the stretched exponential stretching parameter β. The polydispersity model of de Gennes resulted in a prediction for the stretched exponential exponent of 0.25<β<0.33, the Edwards and Vilgis “steric interference” model makes the prediction β = 0.5, and our model for temporary bonding results in the range 0.58<β<0.69. With these disparate results it is appropriate to test whether experimental cartilage stress relaxation data is consistent with any of the models.
The flow of interstitial fluid in cartilage caused by pressure gradients describes flow-dependent stress relaxation. For unconfined compression of a linear poroelastic cylinder of radius r filled with an inviscid fluid, the relaxation kernel is
where the coefficients of the series are An(ν) = ((1-ν)(1-2ν)/(1+ν))(1/[(1-ν)2αn2-(1-2ν)]), ν = Poisson's Ratio (between 0 and 0.5), αn = characteristic equation roots (described subsequently), Ha=Es(((1-ν)/(1+ν)(1-2ν)), Es = equilibrium Young's modulus of the material after complete relaxation, and k = permeability. The characteristic equation is J1(x)-((1-ν) x J0(x)/(1-2ν))=0, where J0(x) and J1(x) are Bessel functions and αn are the roots of the equation. This relaxation kernel is well known as the result of solving the KLM constitutive model for unconfined compression . Herein, we will use this relaxation kernel in combination with polymer fluid relaxation kernels to assess the mechanisms that contribute to the time course of stress relaxation in cartilage explants.
Benchtop and in silico experiments were performed to investigate the relative contributions of monodisperse and polydisperse polymeric models, as well as fluid flow, to cartilage viscoelasticity. Monodisperse and polydisperse (i.e. stretched exponential) models were first evaluated using short-term stress relaxation data on a timescale (less that 60 seconds) that is considerably shorter than expected for fluid flow [7, 27]. Second, an elastic-matrix—inviscid-fluid-flow—stretched-exponential model evaluated full-term stress relaxation data (to 1800 seconds) to determine (1) the time scales where fluid flow and polymer dynamics were dominant and (2) whether the stretching parameter (β) of the model was consistent with any of the specific theoretical polymer mechanisms. Finally, we evaluated the ability of the sticky polymer model to predict results under either low or high ionic concentration .
Five mm diameter osteochondral plugs were aseptically harvested within 6 hours of slaughter from load bearing regions of the lateral femoral condyles of bovine calf stifle joints using a cork borer. The subchondral bone was removed, and the cartilage plugs for the short-time tests were equilibrated at 37° C for 5 days in chemically-defined media composed of DMEM/F-12 supplemented with 0.1% (v/v) BSA, 100 units/ml penicillin, 100 μg/ml streptomycin, and 50 μg/ml ascorbate-2-phosphate . On the day of short-time testing, tissue thickness (4.39±0.0254 mm) was measured with a dial indicator (Harbor Freight, Camarillo, CA Model 623-0VGA). All samples were the same thickness within the precision (25.4 microns) of our measuring equipment, which resulted in less than 3% error in the level of applied strain. Mechanical testing was performed in unconfined compression (Enduratec ELF 3200, Bose Electroforce, Eden Prarie, MN) using polished stainless steel platens. For the long-time experiments, samples were equilibrated under a 5N preload before application of a 5% nominal compression for 1800 s. For the short-time experiments, a 20% compression was applied for 1 minute.
Stress-relaxation tests were performed in unconfined compression. A 20% nominal compressive strain was applied with polished stainless steel platens for 1 minute in a media bath using an Instron 8511 Materials Testing System. Short-time (60 s) stress-relaxation data were obtained from twelve independent samples taken from twelve separate joints. Strain was applied in approximately 0.1 seconds, and force was recorded at 200 Hz for 60 seconds following the application of strain. The relaxation portion of the data (stress data following the peak stress) was used to determine the applicability of the relaxation models, described below. Next, data were fit to two different models, a monodisperse relaxation function
where p=1,3… is an odd index counting the terms used to approximate the infinite sum of the reptation relaxation function, τd is the characteristic disengagement time, and the term S1 is the elastic (long term) portion of the response, and a stretched exponential function
where S0 and S1 are fitting constants (with other terms defined above). Therefore, Equations 7 and 8 approximate the relaxation of cartilage as the combination of an elastic matrix and a viscoelastic polymer. Models were fit using MATLAB (The Mathworks, Natick, MA) to minimize the sum-squared-error (SSE) between the measured stress values and those predicted from the reptation or stretched exponential models. The normalized averaged error (NAE) was then computed to represent the average error per data point, normalized to the maximum stress, resulting from failure of the model to accurately represent the actual data.
The 5 mm diameter samples were collected in the same manner as described above and equilibrated in the same chemically-defined medium for 30 minutes after harvest and prior to testing. A total of 10 samples were tested in this experiment, coming from the medial and lateral patellofemoral grooves midway between the proximal and distal boundaries from 5 independent joints. Testing was performed using polished stainless steel platens in unconfined compression with an Enduratec ELF 3200 electromechanical testing system. The specimens were submerged in culture medium during the compression. A 5N compressive preload was applied for 10 minutes after which the initial height of the sample was measured. Subsequently, a 5% compression was applied at a nominal rate of 5 mm/s (load was completed in 0.2-0.3 s). Force data were sampled at 200 Hz for the first minute of relaxation, and 60 Hz for the following 29 minutes. From comparison of the monodisperse and stretched exponential models it was determined that the latter was better for modeling the full 60 seconds of relaxation. We therefore decided to determine the relative contributions of the stretched exponential, elastic matrix and inviscid fluid flow models to 1800 seconds of relaxation data.
The inviscid fluid-elastic model was composed as the KLM relaxation kernel plus a constant term:
where, S∞, B, and ν were fitting constants, and the remaining terms were as defined above. Thus, this is a function that represents elastic-inviscid fluid (i.e. KLM) behavior for unconfined compression of a homogenous cylinder. This function was fit using a brute-force algorithm in C++ to perform a grid search for the optimal set of model parameters. One hundred terms of the series were used in the curve fit. Initial curvefitting of the KLM function to the entire 1800 seconds of relaxation data found that the KLM function was only able to fit the latter portion of the data. Therefore, we subsequently fit the function to the later portion of the data in a sequential fashion (e.g., for one data set we fit the function to the data from N seconds to 1800 seconds where N=930, 744, 596, 477, …, 33 seconds) until the average error of the fit increased to more than ten times the error of the best fit. This process defined the timescale of relaxation that was accurately modeled by the KLM function. After fitting the KLM function to the data, a stretched exponential was fit to the portion of the stress that was unexplained by fluid flow or the elastic term as:
where, SKWW > 0, τKWW > 0, and β (0 ≤ β < 1) are fitting terms. Brute-force fitting was performed using C++ as described above. Together, Equations 9 and 10 model the relaxation of cartilage as a combination of a linear elastic matrix, inviscid fluid flow, and a polymer fluid. The fitting approach maximized the contribution of the elastic and inviscid fluid flow contribution, in order to determine the minimal possible contribution of the polymer fluid (stretched exponential) model. The mean error of the KLM portion of the fit (i.e., the fit of Equation 10) was calculated as the root-mean-square error between the data and the model. The coefficient of determination (COD) for the total fit was calculated from the force variance of the data (FV) and the residual force variance (RFV) as: COD=1-RFV/FV. Residual force variance is defined as the variance of the residual error of the total fit.
The sticky polymer model predicts that stress-relaxation will proceed faster when the probability of temporary binding is decreased. We sought to validate this model using cartilage stress-relaxation data obtained in either low or high ionic strength solutions, assuming that the high ion concentration would minimize temporary bonds associated with the anionic glycosaminoglycans chains of aggrecan. Cartilage explants were subjected to repeated stress-relaxation tests . The first test occurred in low-ionic strength solution (either 0.15M NaCl or 0.075M CaCl2, n=8 samples). The second test occurred after equilibration in high concentration solution (either 1M NaCl or 0.5M CaCl2).
To fit the temporary binding model to the stress-relaxation data, we used nonlinear optimization in MATLAB. A cost function was constructed using the squared residuals between the data and the model of Eq. 14:
In Eq. 14, p represents the probability that each of the NT binding sites is unbound, i is a counting variable for the sum, lo represents the monomer length, and So, t, and S∞ have been defined above. Note that the model time constant in Equation 11 reflects the hypothesis that aggrecan will temporarily bind to a relatively less mobile component of the cartilage matrix (e.g. collagen). At each data point, the cost function was weighted by the signal-to-noise ratio. To fit the data, we used two methods. First, we used an unbiased method in which all fitting variables (So, S∞, lo, and p) were freely adjustable by the minimization routine for both the first and second datasets. Second, we used a biased method where the results of the unbiased method on the first (low-concentration) dataset were used to define sample-specific lo values, which were then used in a 3-parameter (A, B, and p) fit of each high ion concentration dataset.
Both monodisperse and polydisperse models provided good visual fits to the 60 second data (Figure 2). However, the polydisperse model had a smaller mean SSE than the monodisperse model for the entire data range (3.5±1.6 versus 465.5±60.9 MPa2 p<0.05, t-test), indicating that the stretched exponential is a better model for 60 seconds of relaxation. Sequential fitting of the sixty seconds of data from the twelve specimens found different results between the two relaxation functions. Consistent with previous results , the reptation model showed a minimum error when a short period of data was fitted (0.176±0.061 s; Figure 2B). Additionally, the reptation model had a smaller normalized error than the stretched exponential for times less than approximately one second. The stretched exponential did not show a minimum within the sixty seconds of relaxation (Figure 2B). Therefore, we chose to use the stretched exponential function to test the relative contribution of polymer mechanics and inviscid fluid flow over long term (1800 s) relaxation experiment.
Rapid loading was sufficiently fast to avoid relaxation during compression, as indicated by the linear portion of the relaxation test (r2>0.99) (Figure 3). Further, for all start times, the model (Equation 10) provided good fits to the data (r2>0.975 (Figure 5). However, the KLM portion of the fit showed a rapidly increasing error for start times shorter than ~100 seconds (Figure 4A). Permeability estimated from the KLM portion of the fit (Figure 4D) was dependent upon the fit start time. For start times ~100 seconds and longer, the permeability was similar to literature values . The rapid increase in average error, the presence of an “optimum” start time of approximately 100 seconds, and the convergence of the fitted permeability values for start times greater than 100 seconds, indicated that the KLM model (inviscid fluid—elastic matrix) was able to model the temporal majority (from ~100-1800 s) of the relaxation (Figure 4). However, stress relaxation prior to 100 seconds was best fit by the stretched exponential portion of the fit (Table 1, Figure 5). The average quality of the fit (Table 1) was excellent for the total model (coefficient of determination 0.9998±0.0002) with average error for the KLM portion of the same magnitude as the underlying noise in the load cell (~30 mN). The average time constant for the KLM portion of the curve was significantly longer (p>0.05) than that of the polymer component. The average stretching parameter, β, of the fits (Table 1) was 0.65±0.04, within the range of the sticky model.
Curvefitting the temporary binding model found larger fitted values of unbound probability (p) at higher salt concentrations, providing experimental support for the temporary binding mechanism represented by the model. For samples tested in 0.15 M NaCl p was smaller than samples tested in 1M NaCl (for both biased and unbiased methods p<0.02, Figure 6A-C). For samples tested in 0.075 M CaCl2 p was smaller than samples tested in 0.5M CaCl2 (both p<0.01, Figure 6D-F).
Cartilage is composed of heterogeneous polydisperse biopolymers with multiple complex interactions. To determine the molecular mechanisms of cartilage viscoelasticity, models must be based on specific molecular interactions. In addition, rapid initial loading is necessary to capture the fast relaxation modes that occur on timescales shorter than 100 seconds. In this study, we investigated the relative contributions of monodisperse and polydisperse polymeric models, as well as fluid flow, to cartilage viscoelasticity. Accordingly, the objectives of this study were to (1) examine the timescale of cartilage stress-relaxation using different mechanistic models, and (2) develop and apply a novel (termed “sticky”) polymer mechanics model, based on temporary binding of constituent macromolecules, to cartilage stress-relaxation. These different models necessarily have different numbers of free parameters which reflect the differences in the underlying theories and may result in empirical fitting capabilities that differ between models. More precise separation of the model timescales may be possible using additional methods such as Akaike's Information Criterion .
Stress relaxation in articular cartilage is determined by distinct mechanisms, each associated with different time scales and molecular dynamics. At the earliest times (up to ~0.2 s, Figure 2), the monodisperse reptation polymer model provided the best fit to the data, consistent with our previous research . Interestingly, it is possible that inclusion of constraint-release  or contour length fluctuation  concepts may improve these good fits. This observation supports the hypothesis that for short times after loading molecular motion of cartilage constituents is similar to reptation.
The subsequent timescale (~0.2-100 s) was best described by the stretched exponential model. Most of the stress dissipates during this regime, and one possible mechanism underlying the stretched exponential relaxation is the temporary binding of macromolecules (Equation 11). The finding that curvefits of this model result in larger unbound probabilities for high-salt stress-relaxation provides experimental support for temporary binding mechanisms in cartilage viscoelasticity. The mechanism of temporary binding may be ionic crosslinks between anionic glycosaminoglycans  although other salt-induced phenomenon including osmotic pressure may be relevant. Importantly, this model is consistent with other polymer mechanics treatments of temporary binding [35, 36], although further experimental verification, including the possible use of spectroscopic analysis to quantify temporary binding, would further support the relevance of this mechanism to cartilage mechanics.
Finally, the slowest timescale of cartilage stress-relaxation accounted for the temporal majority (~100-1800 s) and was best described by the fluid-flow model. The results of the curvefits using a combined polymer and KLM approach (1) fit the data well (average coefficient of determination 0.9998), (2) provided estimates of the permeability consistent with literature values, and (3) provided insight about which polymer mechanisms might be active in the medium term.
The observation that stress relaxation data separate into a “fluid” and “polymer” regime at approximately 100 seconds is consistent with the experimental results of Wong et al  who found that for stress relaxation in unconfined compression after rapid loading that apparent stress was proportional to specimen volume only after ~60 seconds of relaxation. Theory for a linear elastic porous material interacting with an inviscid fluid (KLM) finds that stress is linearly related to the specimen volume for unconfined compression of a cylinder. Our analysis (Figure 4) demonstrates that KLM becomes less precise when fitting data prior to 100 seconds of relaxation. In fact, we recapitulate the results of Wong et al (one minute cutoff point) if we use the optimum total coefficient of determination to define the beginning of the fluid-flow timescales (Figure 4B). The consistent findings of this study and Wong et al strongly support the idea that fluid flow occurs on a slow timescale.
The polymer relaxation approach we present here complements other research using finite element models. For example, the paper of Li et al  finds that mesoscopic molecular structures (e.g., Benninghoff arcade-like structures) in cartilage result in a tension-compression anisotropy of total relaxation. However, that study uses phenomenological viscoelasticity theory to describe the behavior of the collagen matrix and was not intended to address the underlying molecular nature of viscoelasticity. In contrast, our current study presents candidate mechanisms that, when further studied, could provide additional insight into the molecular nature of flow-independent viscoelasticity in cartilage.
We found that the time course of stress relaxation after rapid loading in articular cartilage is dominated by distinct mechanisms. The first regime (~0-0.2 s) was best fit by a monodisperse polymer model (reptation). The second regime (~0.2-100 s) was best described by a polydisperse polymer model (stretched exponential). One candidate mechanism underlying this model is temporary binding of macromolecules, and was presented herein as a “sticky” polymer model. The third regime (~100-1800 s) accounted for the temporal majority of the stress-relaxation and was consistent with fluid flow through a porous elastic matrix. Neither polymer dynamics alone nor fluid flow alone were sufficient; both models were necessary to describe the experimental data, indicating that both mechanisms are likely critical in cartilage viscoelasticity.
The authors thank Dr. A. Hari Reddi for kind support in the development of cartilage culture methods and for access to his laboratory facility, and Mr. Roger Zauel for implementation of the C++ curvefitting routines. Code is available from the senior author upon request.
Conflict of Interest: The authors have no conflict of interest.
Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.