Home | About | Journals | Submit | Contact Us | Français |

**|**HHS Author Manuscripts**|**PMC2966345

Formats

Article sections

Authors

Related links

J Biomech Eng. Author manuscript; available in PMC 2010 October 29.

Published in final edited form as:

PMCID: PMC2966345

NIHMSID: NIHMS241446

Gregory C. Thomas,^{1} Anna Asanbaeva,^{2} Pasquale Vena,^{3} Robert L. Sah,^{2} and Stephen M. Klisch^{1}

Corresponding author: Stephen M. Klisch, Ph.D., Associate Professor, Mechanical Engineering Department, California Polytechnic State University, San Luis Obispo, CA 93407, +001(805) 756-1308; FAX +001 (805) 756-1137, Email: ude.yloplac@hcsilks

The publisher's final edited version of this article is available at J Biomech Eng

See other articles in PMC that cite the published article.

A constituent based nonlinear viscoelastic (VE) model was modified from a previous study (Vena et al., *J Biomech Eng*, 2006, 128:449–457) to incorporate a glycosaminoglycan (GAG) - collagen (COL) stress balance using compressible elastic stress constitutive equations specific to articular cartilage (AC). For uniaxial loading of a mixture of quasi-linear VE constituents, time constant and relaxation ratio equations are derived to highlight how a mixture of constituents with distinct quasi-linear VE properties is one mechanism that produces a nonlinear VE tissue. Uniaxial tension experiments were performed with newborn bovine AC specimens before and after ~55 and ~85% GAG depletion treatment with guanidine. Experimental tissue VE parameters were calculated directly from stress relaxation data while intrinsic COL VE parameters were calculated by curve-fitting the data with the nonlinear VE model with intrinsic GAG viscoelasticity neglected. Select tissue and intrinsic COL VE parameters were significantly different between control and experimental groups and correlated with GAG content, suggesting that GAG-COL interactions exist to modulate tissue and COL mechanical properties. Comparison of the results from this and other studies that subjected more mature AC tissue to GAG depletion treatment suggests that GAGs interact with the COL network in a manner that may be beneficial for rapid volumetric expansion during developmental growth while protecting cells from excessive matrix strains. Furthermore, the underlying GAG-COL interactions appear to diminish as the tissue matures, indicating a distinctive remodeling response during developmental growth.

Articular cartilage (AC) covers the ends of articulating bones in synovial joints and functions to transmit loads during joint motion with minimal friction. Approximately 60% to 80% of AC wet weight is interstitial fluid [1] and the remaining weight is a porous solid matrix (SM) composed of glycosaminoglycans (GAGs), collagens (COLs), and other non-collagenous proteins [2]. GAGs provide a fixed negative charge that causes the SM to swell and resist compressive loads and a predominantly type II COL crosslinked network resists tensile and shear loads [3].

The long term goals of this study include improving structure-function relations between biochemical, molecular, and mechanical properties that may be used in continuum mechanics models of AC growth and remodeling [4, 5]. A major obstacle is that AC mechanical properties are difficult to model at the continuum level due to their complexity; they vary with location [2] and depth from the articular surface [6], exhibit anisotropy with respect to directions relative to anatomical split-line direction [7], and exhibit strong tension-compression asymmetry [8, 9]. AC tissue likely experiences finite, multi-dimensional strains when subject to typical loads *in vivo* [10, 11]. These observations suggest the need for finite deformation (i.e. nonlinear) anisotropic stress constitutive equations.

Furthermore, AC generally exhibits both fluid flow dependent and independent (i.e. intrinsic) viscoelastic (VE) behavior [12–17]. In regards to uniaxial tension (UT) stress relaxation, all VE models have a relaxation response that can be described by a time constant, τ, and a relaxation ratio, *R*. Linear and nonlinear VE models may be characterized by τ and *R* parameters that are strain independent and strain dependent, respectively, while in the quasi-linear VE model [18] *R*, but not τ, is strain-dependent [19].^{1} AC tissue VE models have employed linear [16, 20, 21], quasi-linear [12, 14], and nonlinear [17, 22] viscoelasticity yet none of these formulations have used anisotropic elastic stress equations for finite deformations. Recently, [23] introduced a constituent based VE model employing a mixture of quasi-linear VE constituents with anisotropic elastic stress equations for finite deformations. That study showed that a mixture of quasi-linear VE constituents leads to a nonlinear tissue VE response, suggesting one mechanism that causes nonlinear VE behavior. The first objective of this study was to modify the model in [23] to incorporate, and thereby allow study of, AC viscoelasticity.

Previous studies have suggested that VE parameters may change following GAG depletion for AC [24–29] and other biological tissues such as lung tissue [30], temporomandibular joint disc tissue [31], tendon tissue [32], and mitral valve tissue [33]. While most GAG depletion studies with AC involved the use of mature tissue, a recent study [29] suggested that GAG-COL interactions regulate tissue VE properties in a manner dependent on the tissue’s maturational stage. In order to illustrate how a constituent based VE model may be used to improve delineation of AC structure-function relations in tissue growth and remodeling, a second objective of this study was to quantify the VE response for control and GAG depleted immature (i.e. newborn) bovine AC tested in UT. Including comparison of the results with other studies with more mature AC tissue, this objective serves as a first step towards developing a better understanding of how AC remodels during developmental growth to modulate tissue VE properties.

In UT, the effects of flow-dependent viscoelasticity and intrinsic GAG viscoelasticity may be neglected as suggested in several previous studies [22, 34–37].^{2} Consequently, use of such a simplified model with UT experiments facilitates the investigation of several GAG-COL interactions on both tissue and COL VE properties. Here, *indirect GAG interactions* are defined as those that result from a reduction in GAG swelling pressure and a resultant elastic relaxation of the prestressed COL network arising from a GAG-COL stress balance, and *direct GAG interactions* are defined as those that arise due to both covalent and non-covalent molecular interactions between GAGs and the COL network.

The specific aims are to: (1) develop a constituent based nonlinear VE model specific to AC; (2) perform UT stress relaxation experiments for both fresh and GAG depleted immature AC specimens; (3) quantify both experimental tissue and intrinsic COL VE parameters and assess whether significant changes occur due to GAG depletion and if correlations with GAG and COL contents exist.

In this work the AC SM is assumed to occupy a reference configuration *κ _{0}* at time

$$\mathbf{C}={\mathbf{F}}^{T}\mathbf{F}$$

(1)

where **F** is the deformation gradient tensor and *T* the transpose operator. The elastic second Piola-Kirchhoff stress **S** is derived from a strain energy function *W* as

$$\mathbf{S}=\mathit{2}\frac{W\mathbf{C}}{}$$

(2)

and is related to Cauchy stress **T** or first Piola-Kirchhoff stress **P** using

$$J\mathbf{T}={\mathbf{\text{FSF}}}^{T}={\mathbf{\text{PF}}}^{T}.$$

(3)

where the Jacobian *J* is the determinant of **F** i.e. *J* = *det*(**F**).

A GAG-COL stress balance is assumed such that SM stress **S** is the sum of the GAG (**S**^{GAG}) and COL (**S**^{COL}) stresses at every time *t*

$$\mathbf{S}\mathit{(}t\mathit{)}={\mathbf{S}}^{\mathit{\text{GAG}}}\mathit{(}t\mathit{)}+{\mathbf{S}}^{\mathit{\text{COL}}}\mathit{(}t\mathit{)}.$$

(4)

^{3}

Full details of the elastic stress equations are presented in the Appendix. Briefly, the elastic GAG stress ^{e}**S**^{GAG} is based on a two compartmental (i.e. extra- and intra-fibrillar water compartments) isotropic swelling stress model originally proposed in [38] and used in [5, 39] to derive a continuum level GAG stress equation with material constants α_{1} and α_{2}. The elastic COL stress ^{e}**S**^{COL} is based on [40] and contains both an isotropic component ^{e}**S**^{O} with material constant μ and an anisotropic bimodular fiber-reinforced component ^{e}**S**^{BIM} with material constants *γ _{1}, γ_{2}, γ_{3}*, δ which are active (i.e. nonzero) only when their corresponding fibers are in tension.

The GAG and COL constituents are allowed to have distinct reference configurations, ${\kappa}_{\mathit{0}}^{\mathit{\text{GAG}}}\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}{\kappa}_{\mathit{0}}^{\mathit{\text{COL}}}$ (Figure 1). Each constituent is then prescribed an initial deformation, ${\mathbf{F}}_{\mathit{0}}^{\mathit{\text{GAG}}}\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}{\mathbf{F}}_{\mathit{0}}^{\mathit{\text{COL}}}$ to map it to the SM reference configuration κ_{0}. The SM occupies a stress free equilibrium reference configuration at time *t=0*, i.e. **S***(0)*=^{e}**S***(0)* = **0**. However, the GAG isotropic swelling stress model does not have a stress free configuration. Here, the GAG constituent has the same configuration in ${\kappa}_{\mathit{0}}^{\mathit{\text{GAG}}}$ and κ_{0}, which is equivalent to ${\mathbf{F}}_{\mathit{0}}^{\mathit{\text{GAG}}}=\mathbf{I}$ and from (A.1) ^{e}**S**^{GAG}*(0)* = −α_{1}**I**.^{4} To balance this swelling stress, the COL network must support a tensile prestress produced by an initial COL deformation gradient tensor, ${\mathbf{F}}_{\mathit{0}}^{\mathit{\text{COL}}}$, which yields the initial collagen stress tensor, ^{e}**S**^{COL}*(0)*=α_{1}**I**. ${\mathbf{F}}_{\mathit{0}}^{\mathit{\text{COL}}}$ can be determined from equations (A.3) and (A.7). After a deformation **F** is applied, the constituent stresses are calculated relative to their respective reference configurations using

$${\mathbf{F}}^{\mathit{\text{GAG}}}={\mathbf{\text{FF}}}_{\mathit{0}}^{\mathit{\text{GAG}}},{\mathbf{F}}^{\mathit{\text{COL}}}={\mathbf{\text{FF}}}_{\mathit{0}}^{\mathit{\text{COL}}}.$$

(5)

The stress balance hypothesis. The GAG and COL constituents have their own reference configurations, ${\kappa}_{\mathit{0}}^{\mathit{\text{GAG}}}\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}{\kappa}_{\mathit{0}}^{\mathit{\text{COL}}}$. The GAG constituent has the same configuration in ${\kappa}_{\mathit{0}}^{P}$ and κ_{0}, which is equivalent to ${\mathbf{F}}_{\mathit{0}}^{\mathit{\text{GAG}}}=\mathbf{I}$ and ^{e}**S**^{GAG} **...**

To implement constituent based viscoelasticity, the method of [23] for each constituent stress in (4) is implemented as^{5}

$$\mathbf{S}\mathit{(}t\mathit{)}={}^{e}\mathbf{S}^{\mathit{\text{GAG}}}\mathit{(}\mathit{0}\mathit{)}+{\int}_{\mathit{0}}^{t}{G}^{\mathit{\text{GAG}}}\mathit{(}t-\tau \mathit{)}\frac{d{}^{e}\mathbf{S}^{\mathit{\text{GAG}}}}{d\tau}d\tau +{}^{e}\mathbf{S}^{\mathit{\text{COL}}}\mathit{(}\mathit{0}\mathit{)}+{\int}_{\mathit{0}}^{t}{G}^{\mathit{\text{COL}}}(t-\tau )\frac{d{}^{e}\mathbf{S}^{\mathit{\text{COL}}}}{d\tau}d\tau .$$

(6)

Two different relaxation functions, *G ^{GAG}* and

$${G}^{\mathit{\text{GAG}}}\mathit{(}t\mathit{)}=\mathit{1}+{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{GAG}}}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}\mathit{(}-t/{\tau}_{i}^{\mathit{\text{GAG}}}\mathit{)},{G}^{\mathit{\text{COL}}}\mathit{(}t\mathit{)}=\mathit{1}+{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{COL}}}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}\mathit{(}-t/{\tau}_{i}^{\mathit{\text{COL}}}\mathit{)}}}$$

(7)

where *τ _{i}* and

In this study *τ _{i}* and

$$S\mathit{(}t\mathit{)}={}^{e}S\mathit{(}\mathit{0}\mathit{)}+{\displaystyle {\int}_{\mathit{0}}^{t}G\mathit{(}t-\tau \mathit{)}\frac{d{}^{e}S}{d\tau}d\tau .}$$

(8)

Consider the stress response to a step increase in strain and, consequently, a step increase in elastic stress *Δ ^{e}S* at time

$$S\mathit{(}t\mathit{)}={}^{e}S\mathit{(}\mathit{0}\mathit{)}+\Delta {}^{e}S+\Delta {}^{e}S{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}\left(\frac{-\mathit{(}t-{t}_{\mathit{\text{step}}}\mathit{)}}{{\tau}_{i}}\right).}$$

(9)

^{6}

The tissue time constant τ is defined as the amount of time *t* it takes for the tissue stress to relax by 63.2% from its peak value *S(t _{step})* following a step increase in strain and can be calculated from

$$S\mathit{(}\tau \mathit{)}={S}_{\mathit{\text{equil}}}+\mathit{0.368}\mathit{(}S\mathit{(}{t}_{\mathit{\text{step}}}\mathit{)}-{S}_{\mathit{\text{equil}}}\mathit{)}$$

(10)

where the equilibrium stress is

$${S}_{\mathit{\text{equil}}}={}^{e}S^{\mathit{\text{GAG}}}\mathit{(}\mathit{0}\mathit{)}+\Delta {}^{e}S^{\mathit{\text{GAG}}}+{}^{e}S^{\mathit{\text{COL}}}\mathit{(}\mathit{0}\mathit{)}+\Delta {}^{e}S^{\mathit{\text{COL}}}.$$

(11)

The peak stress is obtained upon evaluating (9) at *t=t _{step}* for each constituent

$$S\mathit{(}{t}_{\mathit{\text{step}}}\mathit{)}={}^{e}S^{\mathit{\text{GAG}}}\mathit{(}\mathit{0}\mathit{)}+\Delta {}^{e}S^{\mathit{\text{GAG}}}\mathit{(}\mathit{1}+{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{GAG}}}\mathit{)}+{}^{e}S^{\mathit{\text{COL}}}\mathit{(}\mathit{0}\mathit{)}+\Delta {}^{e}S^{\mathit{\text{COL}}}\mathit{(}\mathit{1}+{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{COL}}}\mathit{)}}}$$

(12)

$$S\mathit{(}\tau \mathit{)}={S}_{\mathit{\text{equil}}}+\mathit{0.368}\mathit{(}\Delta {}^{e}S^{\mathit{\text{GAG}}}\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{GAG}}}+\Delta {}^{e}S^{\mathit{\text{GAG}}}\phantom{\rule{thinmathspace}{0ex}}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{COL}}}\mathit{)}.}}$$

(13)

Evaluating (9) at *t=τ* for each constituent and using (11) one obtains

$$S\mathit{(}\tau \mathit{)}={S}_{\mathit{\text{equil}}}+\Delta {}^{e}S^{\mathit{\text{GAG}}}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{GAG}}}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}\left(\frac{-\mathit{(}\tau -{t}_{\mathit{\text{step}}}\mathit{)}}{{\tau}_{i}^{\mathit{\text{GAG}}}}\right)+\Delta {}^{e}S^{\mathit{\text{COL}}}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{COL}}}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}\left(\frac{-\mathit{(}\tau -{t}_{\mathit{\text{step}}}\mathit{)}}{{\tau}_{i}^{\mathit{\text{COL}}}}\right)}}$$

(14)

Equating (13–14) allows the calculation of τ from

$$\mathit{0.368}\mathit{(}\Delta {}^{e}S^{\mathit{\text{GAG}}}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{GAG}}}+\Delta {}^{e}S^{\mathit{\text{COL}}}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{COL}}}\mathit{)}=}}\Delta {}^{e}S^{\mathit{\text{GAG}}}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{GAG}}}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}\left(\frac{-\mathit{(}\tau -{t}_{\mathit{\text{step}}}\mathit{)}}{{\tau}_{i}^{\mathit{\text{COL}}}}\right)+\Delta {}^{e}S^{\mathit{\text{COL}}}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{COL}}}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}\left(\frac{-\mathit{(}\tau -{t}_{\mathit{\text{step}}}\mathit{)}}{{\tau}_{i}^{\mathit{\text{COL}}}}\right)}}$$

(15)

For a single constituent (15) reduces to

$$\mathit{0.368}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}={\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}\left(\frac{-\mathit{(}\tau -{t}_{\mathit{\text{step}}}\mathit{)}}{{\tau}_{i}}\right)}}$$

(16)

which reveals that when the relaxation function parameters are constants, τ is strain-independent, a characteristic of quasi-linear viscoelasticity. Conversely, for constituent based viscoelasticity (15) reveals that τ depends on the applied UT strain (because of the dependence on stress in (15)) a characteristic of nonlinear viscoelasticity.

The relaxation ratio *R* is defined as equilibrium stress divided by peak stress

$$R=\frac{{S}_{\mathit{\text{equil}}}}{S\mathit{(}{t}_{\mathit{\text{step}}}\mathit{)}}.$$

(17)

$$R=\frac{{}^{e}{S}^{\mathit{\text{GAG}}}\mathit{(}\mathit{0}\mathit{)}+\Delta {}^{e}S^{\mathit{\text{GAG}}}+{}^{e}S^{\mathit{\text{COL}}}\mathit{(}\mathit{0}\mathit{)}+\Delta {}^{e}S^{\mathit{\text{COL}}}}{{}^{e}{S}^{\mathit{\text{GAG}}}\mathit{(}\mathit{0}\mathit{)}+\Delta {}^{e}S^{\mathit{\text{GAG}}}\left(\mathit{1}+{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{GAG}}}}\right)+{}^{e}S^{\mathit{\text{COL}}}\mathit{(}\mathit{0}\mathit{)}+\Delta {}^{e}S^{\mathit{\text{COL}}}\left(\mathit{1}+{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{COL}}}}\right)}$$

(18)

For a single constituent (18) reduces to

$$R=\frac{{}^{e}S\mathit{(}\mathit{0}\mathit{)}+\Delta {}^{e}S}{{}^{e}S\mathit{(}\mathit{0}\mathit{)}+\Delta {}^{e}S\left(\mathit{1}+{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}}\right)}$$

(19)

which reveals that when the relaxation function parameters are constants, *R* is strain-dependent as it is a function of applied UT strain when the initial strain is not zero, a characteristic of quasi-linear viscoelasticity.^{7} For constituent based viscoelasticity, (18) reveals how *R* additionally depends on the constituent properties.

In the current study, the method of [23] is modified in three ways: (1) the initial equilibrium elastic stresses of the constituents, ^{e}**S**^{GAG}*(0)* and ^{e}**S**^{COL}*(0)*, are nonzero in a SM stress free reference configuration; (2) the first term in the relaxation functions equals *1* instead of a variable “*g*_{∞}” so that the equilibrium elastic stress is obtained as *t* → ∞; and (3) the incompressibility assumption is not used. These modifications result in a slightly different time discretization scheme summarized here.^{8}

To simplify the final result, a recursive expression for time-dependent matrices, **c**^{i}*(t)*, and a scalar, *c ^{COL}*, specific to COL are defined as

$${\mathbf{c}}^{i}\mathit{(}t\mathit{)}=\mathit{\text{exp}}\mathit{(}-\Delta t/{\tau}_{i}^{\mathit{\text{COL}}}\mathit{)}{\mathbf{c}}^{i}\mathit{(}t-\Delta t\mathit{)}+{g}_{i}^{\mathit{\text{COL}}}{\tau}_{i}^{\mathit{\text{COL}}}\mathit{(}\mathit{1}-\mathit{\text{exp}}\mathit{(}-\Delta t/{\tau}_{i}^{\mathit{\text{COL}}}\mathit{)}\mathit{)}\left(\frac{{}^{e}{\mathbf{S}}^{\mathit{\text{COL}}}\mathit{(}t\mathit{)}-{}^{e}{\mathbf{S}}^{\mathit{\text{COL}}}\mathit{(}t-\Delta t\mathit{)}}{\Delta t}\right)$$

(20)

and

$${c}^{\mathit{\text{COL}}}=\mathit{1}+\frac{\mathit{1}}{\Delta t}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{COL}}}{\tau}_{i}^{\mathit{\text{COL}}}\mathit{(}\mathit{1}-\mathit{\text{exp}}\mathit{(}-\Delta t/{\tau}_{i}^{\mathit{\text{COL}}}\mathit{)}.}$$

(21)

The time discretization procedure used to calculate **S**^{COL}(*t + Δt*) after a time increment *Δt* is

$${\mathbf{S}}^{\mathit{\text{COL}}}\mathit{(}t+\Delta t\mathit{)}={\mathit{(}{c}^{\mathit{\text{COL}}}\mathit{)}}^{e}{\mathbf{S}}^{\mathit{\text{COL}}}\mathit{(}t+\Delta t\mathit{)}+{\mathit{(}\mathit{1}-{c}^{\mathit{\text{COL}}}\mathit{)}}^{e}{\mathbf{S}}^{\mathit{\text{COL}}}\mathit{(}t\mathit{)}+{\displaystyle \sum _{i=\mathit{1}}^{n}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}\mathit{(}-\Delta t/{\tau}_{i}^{\mathit{\text{COL}}}\mathit{)}{\mathbf{c}}^{i}\mathit{(}t\mathit{)}.}$$

(22)

An analogous expression may be obtained for the GAG constituent.

For UT analysis, GAG viscoelasticity is neglected upon assuming that the AC UT VE response is dominated by intrinsic COL viscoelasticity.^{9} The *UT model* is obtained by setting all GAG relaxation parameters ${g}_{i}^{\mathit{\text{GAG}}}$ equal to zero. As seen below, a different number of Prony series terms for the COL relaxation function are needed for different experimental groups; thus, it is not possible to investigate significant differences for ${g}_{i}^{\mathit{\text{COL}}}\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}{\tau}_{i}^{\mathit{\text{COL}}}$ between groups. Instead, an aggregate COL time constant and an aggregate COL relaxation ratio are defined as functions of ${g}_{i}^{\mathit{\text{COL}}}\phantom{\rule{thinmathspace}{0ex}}\text{and}\phantom{\rule{thinmathspace}{0ex}}{\tau}_{i}^{\mathit{\text{COL}}}$ to render COL viscous properties that may be compared between groups with different number of relaxation function terms.

For the UT model, the amount of time $t={\tau}_{A}^{\mathit{\text{COL}}}$ it takes to relax by 63.2% can be calculated from (15) which reduces upon neglecting intrinsic GAG VE to

$$\mathit{0.368}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{COL}}}={\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{COL}}}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}\left(\frac{-\mathit{(}{\tau}_{A}^{\mathit{\text{COL}}}-{t}_{\mathit{\text{step}}}\mathit{)}}{{\tau}_{i}^{\mathit{\text{COL}}}}\right)}}$$

(23)

where ${\tau}_{A}^{\mathit{\text{COL}}}$ is defined as the *aggregate COL time constant*. The relaxation ratio *R* from (16) reduces upon neglecting intrinsic GAG VE and using a SM stress-free reference configuration to

$$R=\frac{\Delta {}^{e}S^{\mathit{\text{GAG}}}+\Delta {}^{e}S^{\mathit{\text{COL}}}}{\Delta {}^{e}S^{\mathit{\text{GAG}}}+\Delta {}^{e}S^{\mathit{\text{COL}}}\left(\mathit{1}+{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{COL}}}}\right)}$$

(24)

suggesting a definition for an *aggregate COL relaxation ratio* ${R}_{A}^{\mathit{\text{COL}}}$ as

$${R}_{A}^{\mathit{\text{COL}}}=\frac{\mathit{1}}{\mathit{1}+{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}^{\mathit{\text{COL}}}}}.$$

(25)

^{10}

Recall that, relative to a stress- and strain-free configuration, the time constant and relaxation ratio depend on the applied UT strain for a nonlinear VE material but do not depend on the applied UT strain for a quasi-linear VE material. Thus, for the UT model (24) and (23) reveal that the relaxation ratio, but not the time constant, depends on the applied UT strain. Consequently, this model falls somewhere between nonlinear and quasi-linear viscoelasticity.

From an earlier study [41], the UT stress relaxation response was measured on groups of adjacent paired AC explants having 0, ~55, and ~85 % GAG depletion. Newborn (1–3 week old) bovine AC patellofemoral groove explants were harvested from ~0.6 mm below the AC surface and sliced to a thickness of ~0.25 mm. 12 explants were harvested, corresponding to four specimens per group. The control group (no GAG depletion i.e. GD-0) did not receive enzyme treatment. Explants for the ~85% GAG extraction group (GD-85) were rinsed in guanidine HCl (Gnd) to remove ~85% GAG mass. Explants for the ~55% GAG extraction group (GD-55) were rinsed in Gnd to remove ~85% GAG mass and then soaked in a solution of Gnd saturated with cartilage extract to replete GAG mass to ~45% of the average GD-0 GAG mass.

Explants were punched into tapered tensile specimens (length ~4 mm and cross sectional gage area ~0.25 × 0.8 mm) and oriented such that the direction of uniaxial loading was along the anterior-posterior direction, approximately perpendicular to the split line direction [42]. The thickness of each tensile specimen was measured at three locations in the gage region and averaged for cross-sectional area calculations. Specimens were loaded into a testing fixture, pulled to a tare-load of 0.05 N (equivalent to a stress of ~0.2 MPa), allowed to relax, and then pulled to 10% strain under displacement control over a 40 second ramp. Displacement was held at 10% strain during a 900 second stress relaxation period. The measured load and displacement were converted to first Piola-Kirchhoff stress (load / original cross-sectional area) and strain (elongation / initial clamp to clamp distance).

Each tensile specimen and residual tissue obtained during preparation were saved for biochemical analysis. Samples were solubilized with proteinase K [43] and analyzed to quantify content of sulfated GAG [44] and hydroxyproline [45]. Hydroxyproline content was converted to COL content by assuming a mass ratio of COL/hydroxyproline = 7.25 for bovine cartilage [46, 47]. Biochemical parameters were normalized to tissue wet weight (WW) to represent constituent content.

Experimental stress data were used to calculate experimental tissue VE parameters that are independent of the proposed VE model. Equilibrium stress was defined as the stress at the end of the relaxation phase. Experimental Young’s modulus *E ^{exp}* was defined as equilibrium stress divided by applied strain (0.10). The experimental time constant

Experimental stress and biochemical data were used to calculate elastic GAG-COL parameters in elastic constitutive equations using MATLAB. GAG material constants *α _{1}* and

Due to a lack of comprehensive data, the COL elastic parameters (μ, Ф_{+12}, Ф_{+13}, Ф_{+23}) were held constant among all specimens. In this model, μ corresponds to the tissue shear modulus for small strain elasticity theory and was estimated from torsional shear data for newborn bovine AC [48] as 0.11 MPa. Values of (Ф_{+12}, Ф_{+13}, Ф_{+23}) were obtained from [5, 40] as (45 deg., 35 deg., 35 deg.) which produce reasonable predictions for anisotropic and asymmetric Young’s moduli and Poisson’s ratios for newborn bovine and adult human AC. The remaining elastic COL material constants were constrained by γ_{1} = γ_{2} = 2γ_{3} = δ which produce reasonable predictions for mechanical properties of newborn bovine AC [5]. Consequently, there remains only one adjustable COL elastic parameter, γ_{1}, which was calculated on a specimen-specific basis. Specifically, a MATLAB code was developed that solves for the initial COL deformation gradient tensor, ${\mathbf{F}}_{\mathit{0}}^{\mathit{\text{COL}}}$, and γ_{1} that yields a SM stress free reference configuration (i.e. **S***(0)* = **0**), and matches the measured equilibrium stress at 10% strain.

After determining ${\mathbf{F}}_{\mathit{0}}^{\mathit{\text{COL}}}$ and γ_{1}, the VE time discretization procedure (20–22) with two terms in the Prony series for the COL relaxation function were coded in MATLAB to solve the UT boundary-value problem with traction free boundary conditions on the free surfaces using first Piola-Kirchhoff stress. For each specimen, an iterative optimization procedure was implemented to determine COL VE parameters $({\tau}_{\mathit{1}}^{\mathit{\text{COL}}},{g}_{\mathit{1}}^{\mathit{\text{COL}}},{\tau}_{\mathit{2}}^{\mathit{\text{COL}}},{g}_{\mathit{2}}^{\mathit{\text{COL}}})$ that minimize the squared sum of residual errors between experimental and theoretical stress values. A parameter study revealed that the optimal solution for the GD-0 and GD-55 groups did not depend on initial guesses for $({\tau}_{\mathit{1}}^{\mathit{\text{COL}}},{g}_{\mathit{1}}^{\mathit{\text{COL}}},{\tau}_{\mathit{2}}^{\mathit{\text{COL}}},{g}_{\mathit{2}}^{\mathit{\text{COL}}})$; however, the solution did depend on initial guesses for the GD-85 group. Consequently, for GD-85 specimens one term in the Prony series for the COL relaxation function was used.^{11} After optimization, the aggregate time constant ${\tau}_{A}^{\mathit{\text{COL}}}$ and relaxation ratio ${R}_{A}^{\mathit{\text{COL}}}$ were calculated using (23) and (25), respectively.

A parameter study was conducted in order to differentiate between effects of the different GAG-COL interaction mechanisms on the tissue VE response. The MATLAB code was modified to obtain the UT boundary value-problem solution using first Piola-Kirchhoff stress and averaged values of the elastic and viscous parameters for the GD-0 group to identify a baseline case. Additional solutions were obtained to model progressive changes to the GD-0 VE properties due to treatment. First, the effect of an altered GAG-COL stress balance due to reduced GAG content on tissue VE properties were isolated: averaged GAG elastic parameters (*α _{1},α_{2}*) were adjusted to the GD-55 and GD-85 values while the COL VE parameters were kept equal to the GD-0 values. Second, the incremental effect of direct GAG interactions on COL elasticity was isolated: in addition to (

The effects of GAG depletion treatment on experimental tissue VE parameters (*E ^{exp}, τ^{exp}, R^{exp}*), aggregate VE parameters s $({\tau}_{A}^{\mathit{\text{COL}}},{R}_{A}^{\mathit{\text{COL}}})$, and COL elastic parameter (

Experimental tissue VE parameters (*E ^{exp}, τ^{exp}, R^{exp}*) increased (p<0.01), decreased (p<0.001), and increased (p<0.001), respectively, due to increasing GAG depletion treatment (Table 1; Figure 3). All three (

The effects of GAG depletion treatment on experimental tissue VE parameters (*E*^{exp}, τ^{exp}, R^{exp}); mean + 1 standard deviation values shown. GD-0 = control group with no GAG depletion. GD-55 and GD-85 = experimental groups with ~ 55 and 85% GAG depletion, **...**

Experimental tissue VE parameters (*E*^{exp}, τ^{exp}, R^{exp}) vs. GAG and COL contents. GD-0 = control group with no GAG depletion. GD-55 and GD-85 = experimental groups with ~ 55 and 85% GAG depletion, respectively. Linear regression results only shown **...**

Experimental GAG and COL contents normalized to tissue wet weight (WW), GAG elastic material parameters (*α*_{1},α_{2}), COL VE parameters $({\gamma}_{1},{\tau}_{\mathit{1}}^{\mathit{\text{COL}}},{g}_{\mathit{1}}^{\mathit{\text{COL}}},{\tau}_{\mathit{2}}^{\mathit{\text{COL}}},{g}_{\mathit{2}}^{\mathit{\text{COL}}})$, and curve-fit R^{2} values. Mean + 1 standard deviation **...**

Unique VE model parameters for the GD-0 and GD-55 specimens were obtained with two Prony series terms in the COL relaxation function; as mentioned above, unique solutions for the GD-85 specimens were obtained with one Prony series term (Table 1, Figure 5). The COL elastic parameter *γ _{1}* increased due to increasing GAG depletion treatment (p<0.001) (Table 1; Figure 6). The COL VE parameter ${g}_{\mathit{2}}^{\mathit{\text{COL}}}$ decreased due to increasing GAG depletion treatment (p<0.01) (Table 1; Figure 6), while a significant difference was almost detected for ${\tau}_{\mathit{2}}^{\mathit{\text{COL}}}$ (p=0.06) (Table 1; Figure 6). The aggregate VE relaxation parameter ${R}_{\mathit{2}}^{\mathit{\text{COL}}}$ increased due to increasing GAG depletion treatment (p<0.0001) while a significant difference was almost detected for ${\tau}_{A}^{\mathit{\text{COL}}}$ (p=0.08) (Figure 6). The VE parameters $({\gamma}_{1},{g}_{\mathit{2}}^{\mathit{\text{COL}}},{R}_{A}^{\mathit{\text{COL}}})$ were significantly correlated with GAG content (p<0.01, 0.0001, and 0.001, respectively) (Figure 7). None of the VE parameters were significantly correlated with COL content although positive trends existed for

Specimen-specific curve-fits of the constituent based VE model. GD-0 = control group with no GAG depletion. GD-55 and GD-85 = experimental groups with ~ 55 and 85% GAG depletion, respectively.

The effects of GAG depletion treatment on COL VE parameters $({\gamma}_{1},{\tau}_{\mathit{1}}^{\mathit{\text{COL}}},{g}_{\mathit{1}}^{\mathit{\text{COL}}},{\tau}_{\mathit{2}}^{\mathit{\text{COL}}},{g}_{\mathit{2}}^{\mathit{\text{COL}}},{\tau}_{A}^{\mathit{\text{COL}}},{R}_{A}^{\mathit{\text{COL}}})$; mean + 1 standard deviation values shown. GD-0 = control group with no GAG depletion. GD-55 and GD-85 = experimental **...**

COL VE parameters $({\gamma}_{1},{\tau}_{\mathit{1}}^{\mathit{\text{COL}}},{g}_{\mathit{1}}^{\mathit{\text{COL}}},{\tau}_{\mathit{2}}^{\mathit{\text{COL}}},{g}_{\mathit{2}}^{\mathit{\text{COL}}},{\tau}_{A}^{\mathit{\text{COL}}},{R}_{A}^{\mathit{\text{COL}}})$ vs. GAG and COL contents. GD-0 = control group with no GAG depletion. GD-55 and GD-85 = experimental groups with ~ 55 and 85% GAG depletion, respectively. Linear **...**

The parameter studies using averaged GAG elastic parameters (*α _{1},α_{2}*) for the GD-55 and GD-85 groups with COL VE parameters $({\gamma}_{1},{\tau}_{\mathit{1}}^{\mathit{\text{COL}}},{g}_{\mathit{1}}^{\mathit{\text{COL}}},{\tau}_{\mathit{2}}^{\mathit{\text{COL}}},{g}_{\mathit{2}}^{\mathit{\text{COL}}})$ for the GD-0 group showed that indirect GAG interactions resulting from the stress balance hypothesis modulate peak and equilibrium stresses but not relaxation behavior (i.e., time constant and relaxation ratio) (Figure 8). The parameter studies using averaged elastic parameters (

A constituent based nonlinear VE model for AC tissue was developed using quasi-linear VE GAG and COL constituents. For uniaxial loading it was shown that this model produces a nonlinear VE response for which both the time constant and relaxation ratio are strain-dependent, as opposed to a quasi-linear VE response for which only the relaxation ratio is strain-dependent. More specifically, the quasi-linear VE relaxation ratio is strain-dependent when the initial strain is not zero but is strain-independent when the initial strain is zero. That conclusion highlights one mechanism that contributes to nonlinear VE properties; i.e. the presence of multiple constituents that contribute to distinct VE parameters. For the UT model that neglected GAG VE for the analysis of UT experiments, it was shown that the time constant is strain-independent while the relaxation ratio is strain-dependent even when the initial strain is zero, in contrast to quasi-linear viscoelasticity. Consequently, the UT model falls somewhere between quasi-linear and nonlinear viscoelasticity.

The results for the experimental tissue VE parameters (*E ^{exp}, τ^{exp}, R^{exp}*) when compared with results from previous studies that included more mature AC tissue tested in UT suggest that there exist GAG-COL interactions that remodel during developmental growth. Specifically,

In addition to the studies mentioned above for AC tissue tested in UT, GAG depleted tissue was tested in unconfined compression from mature (1–2 year old) bovine AC and modeled using small-strain poroviscoelastiy [27]. Those results showed that VE time constants decreased due to GAG depletion, in agreement with the results of this study, and the unconfined compression elastic modulus decreased due to GAG depletion, in contrast with the results of this study. Again, this difference is likely due to the distinct maturational stages of the specimens.

Collectively, these results suggest that during maturation there exist GAG-COL interactions that provide for a highly compliant (i.e. low tensile modulus) SM. The biological significance of a highly compliant SM is a key property that underlies a mechanism for rapid volumetric expansion during developmental growth. Further, GAG-COL interactions may provide a lower relaxation ratio (i.e. higher peak stress for a given equilibrium stress) and enhanced viscous properties as illustrated in Figure 8. As a creep analysis would reveal lower SM strains during repeated load applications, the biological significance of enhanced viscous properties may be that it is a mechanism that may protect cells against excessive matrix strains, that are know to cause cell death [49], during repeated loading.

A novel outcome of this study was that the experimental tissue VE response with and without GAG depletion treatment was interpreted in terms of a constituent based VE model. Although in general both GAG and COL intrinsic VE parameters may be active during complex *in vivo* loading conditions, in this study intrinsic GAG viscoelasticity was neglected due to the assumption that the UT VE response is dominated by intrinsic COL viscoelasticity. This assumption was invoked in previous studies [34–37] and is supported by experimental time constants and dynamic moduli measured in shear for GAG solutions [26, 50, 51] which were several orders of magnitude less than time constants and dynamic moduli measured for intact AC [26, 52]. Those previous results, combined with the hypothesis that shear resistance is provided by tensile stresses produced in the GAG inflated COL network [53], suggest that intrinsic GAG viscoelasticity can be neglected in UT. Indeed, this assumption was supported by pilot analyses for the current study with four VE coefficients $({g}_{\mathit{1}}^{\mathit{\text{GAG}}},{\tau}_{\mathit{1}}^{\mathit{\text{GAG}}},{g}_{\mathit{1}}^{\mathit{\text{COL}}},{\tau}_{\mathit{1}}^{\mathit{\text{COL}}})$ as unique solutions for these coefficients were not found and the best-fit solutions yielded tensile GAG stresses of ~0.3 MPa which may not be physically relevant as GAGs are thought to provide primarily compressive resistance.^{12}

If the assumption of neglecting intrinsic GAG viscoelasticity is valid, then the results of this study suggest that GAGs interact with the COL network to affect intrinsic COL viscoelasticity as select COL VE parameters changed due to GAG depletion and were correlated with GAG content. Interestingly, the GD-0 and GD-55 specimens exhibited distinct short- and long-term relaxation regimes as evidenced by the ability of the model to find unique solutions using two Prony series terms in the COL relaxation function, whereas this distinction was not possible to detect for GD-85 specimens. The parameter studies suggest that indirect GAG interactions arising from a GAG-COL stress balance as well as direct GAG interactions on COL elastic parameters modulate the peak and equilibrium stresses while direct GAG interactions on COL viscous parameters modulate the relaxation response.

Although the elastic moduli were correlated with GAG content they were not correlated with COL content, in contrast to results in [24]. The protocols, which included two treatment groups that depleted GAGs by ~ 55 and 85%, produced a wide range of GAG content which made it possible to detect correlations with GAG content with a limited number of specimens in each group. However, all groups had comparable COL content making it difficult to detect correlations with COL content; since positive trends existed for *E ^{exp}* (p=0.19),

A previous study [54] suggested that although tissue VE is affected by GAG depletion for aortic valve tissue, the mechanism for such changes is decreasing water content that accompanies GAG depletion. However, additional statistical analyses for the present study (results of which are not presented) revealed no significant correlations between any of the VE parameters and water content, suggesting that GAG-COL interactions were more prominent than water-COL interactions for the experiments performed in this study.

In summary, a constituent based nonlinear VE model was developed that modified a previous model [23] and used to quantify AC tissue remodeling due to altered GAG-COL interactions. The results suggest that the GAGs interact with the COL network to affect tissue VE properties in immature AC tissue. Upon considering results from other studies with more mature tissue, it appears that GAG-COL interactions exist that may be beneficial for rapid volumetric expansion during developmental growth while protecting cells from excessive matrix strains. Furthermore, such interactions appear to diminish as the tissue matures, indicative of a remodeling response during developmental growth.

Funding was received from the National Science Foundation (RS, SK), National Institutes of Health (RS, SK), and Howard Hughes Medical Institute (UCSD for RS).

The elastic stress constitutive law is based on a bimodular polyconvex fiber-reinforced anisotropic strain energy function for AC [40] where an isotropic GAG matrix is reinforced with a COL fiber constituent. The elastic GAG stress ^{e}**S**^{GAG} is based on a two compartmental (i.e. extra- and intra-fibrillar water compartments) isotropic swelling stress model for AC GAGs originally proposed in [38] and used in [5, 39] to derive the polyconvex continuum level model

$${}^{e}{\mathbf{S}}^{\mathit{\text{GAG}}}=\frac{{\alpha}_{\mathit{1}}}{{\mathit{(}{J}^{\mathit{\text{GAG}}}\mathit{)}}^{{\alpha}_{\mathit{2}}-\mathit{1}}}{\mathbf{C}}^{\mathit{\text{GAG}}}$$

(A.1)

where α_{1} and α_{2} are material constants.

The elastic COL stress ^{e}**S**^{COL} is based on [40] and contains both an isotropic component ^{e}**S**^{O} and a anisotropic bimodular component ^{e}**S**^{BIM}:

$${}^{e}{\mathbf{S}}^{\mathit{\text{COL}}}={}^{e}{\mathbf{S}}^{\mathit{0}}+{}^{e}{\mathbf{S}}^{\mathit{\text{BIM}}}.$$

(A.2)

^{e}**S**^{O} is defined as

$${}^{e}{\mathbf{S}}^{O}=\mu \mathit{(}\mathbf{I}-{\mathit{(}{\mathbf{C}}^{\mathit{\text{COL}}}\mathit{)}}^{-\mathit{1}}\mathit{)}$$

(A.3)

where μ is a material constant and −*1* denotes the inverse operator. The fiber reinforcement of the SM appears in ^{e}**S**^{BIM}. In this model, there are three primary fiber directions along mutually orthogonal unit vectors **E**_{1}, **E**_{2}, **E**_{3} which generate structural tensors

$${\mathbf{M}}_{\mathit{1}}={\mathbf{E}}_{\mathit{1}}{\mathbf{E}}_{\mathit{1}},{\mathbf{M}}_{\mathit{2}}={\mathbf{E}}_{\mathit{2}}{\mathbf{E}}_{\mathit{2}},{\mathbf{M}}_{\mathit{3}}={\mathbf{E}}_{\mathit{3}}{\mathbf{E}}_{\mathit{3}}$$

(A.4)

where is the dyadic product. Also, there are two secondary fiber directions in each of the three planes created by the primary unit vectors **E**_{1}, **E**_{2}, **E**_{3} (planes 1–2, 2–3, 1–3) along unit vectors (**E**_{±12}, **E**_{±23}, **E**_{±13}) defined by angles (Ф_{±12}, Ф_{±23}, Ф_{±13}) as

$${\mathbf{E}}_{\mathit{\pm}\mathit{12}}=\mathit{\text{cos}}{\Phi}_{\mathit{\pm}\mathit{12}}{\mathbf{E}}_{\mathit{1}}\mathbf{\pm}\mathit{\text{sin}}{\Phi}_{\mathit{\pm}\mathit{12}}{\mathbf{E}}_{\mathit{2}},{\mathbf{E}}_{\mathit{\pm}\mathit{23}}=\mathit{\text{cos}}{\Phi}_{\mathit{\pm}\mathit{23}}{\mathbf{E}}_{\mathit{2}}\mathbf{\pm}\mathit{\text{sin}}{\Phi}_{\mathit{\pm}\mathit{23}}{\mathbf{E}}_{\mathit{3}},{\mathbf{E}}_{\mathit{\pm}\mathit{13}}=\mathit{\text{cos}}{\Phi}_{\mathit{\pm}\mathit{13}}{\mathbf{E}}_{\mathit{1}}\mathbf{\pm}\mathit{\text{sin}}{\Phi}_{\mathit{\pm}\mathit{13}}{\mathbf{E}}_{\mathit{3}}.$$

(A.5)

These generate structural tensors

$$\begin{array}{l}{\mathbf{M}}_{\pm \mathit{12}}={\mathit{\text{cos}}}^{\mathit{2}}{\varphi}_{\mathit{12}}{\mathbf{E}}_{\mathit{1}}{\mathbf{E}}_{\mathit{1}}+{\mathit{\text{sin}}}^{\mathit{2}}{\varphi}_{\mathit{12}}{\mathbf{E}}_{\mathit{2}}{\mathbf{E}}_{\mathit{2}}\pm \mathit{\text{cos}}{\varphi}_{\mathit{12}}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{sin}}{\varphi}_{\mathit{12}}{\mathit{(}\mathbf{E}}_{\mathit{1}}{\mathbf{E}}_{\mathit{2}}+{\mathbf{E}}_{\mathit{2}}{\mathbf{E}}_{\mathit{1}}\mathit{)}{\mathbf{M}}_{\pm \mathit{23}}={\mathit{\text{cos}}}^{\mathit{2}}{\varphi}_{\mathit{23}}{\mathbf{E}}_{\mathit{2}}{\mathbf{E}}_{\mathit{2}}+{\mathit{\text{sin}}}^{\mathit{2}}{\varphi}_{\mathit{23}}{\mathbf{E}}_{3}{\mathbf{E}}_{\mathit{3}}\pm \mathit{\text{cos}}{\varphi}_{\mathit{12}}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{sin}}{\varphi}_{\mathit{12}}{\mathit{(}\mathbf{E}}_{\mathit{2}}{\mathbf{E}}_{\mathit{3}}+{\mathbf{E}}_{\mathit{3}}{\mathbf{E}}_{\mathit{2}}\mathit{)}{\mathbf{M}}_{\pm \mathit{13}}={\mathit{\text{cos}}}^{\mathit{2}}{\varphi}_{\mathit{13}}{\mathbf{E}}_{\mathit{1}}{\mathbf{E}}_{\mathit{1}}+{\mathit{\text{sin}}}^{\mathit{2}}{\varphi}_{\mathit{13}}{\mathbf{E}}_{\mathit{3}}{\mathbf{E}}_{\mathit{3}}\pm \mathit{\text{cos}}{\varphi}_{\mathit{13}}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{sin}}{\varphi}_{\mathit{13}}{\mathit{(}\mathbf{E}}_{\mathit{1}}{\mathbf{E}}_{\mathit{3}}+{\mathbf{E}}_{\mathit{3}}{\mathbf{E}}_{\mathit{1}}\mathit{)}.\hfill \hfill \hfill \end{array}$$

(A.6)

The strain energy function *W ^{BIM}* is modified from [40] to model a linear stress-strain response typical of newborn bovine tissue. The resulting expression for

$${}^{e}{\mathbf{S}}^{\mathit{\text{BIM}}}=\left({\displaystyle \sum _{i=\mathit{1}}^{\mathit{3}}{\gamma}_{i}\mathit{[}{\lambda}_{i}^{\mathit{2}}\mathit{]}\mathit{(}{\lambda}_{i}^{\mathit{2}}-\mathit{1}\mathit{)}{\mathbf{M}}_{i}}\right)+\delta \mathit{[}{\lambda}_{\pm \mathit{12}}^{\mathit{2}}\mathit{]}\mathit{(}{\lambda}_{\pm \mathit{12}}^{\mathit{2}}-\mathit{1}\mathit{)}{\mathbf{M}}_{\pm \mathit{12}}+\delta \mathit{[}{\lambda}_{\pm \mathit{23}}^{\mathit{2}}\mathit{]}\mathit{(}{\lambda}_{\pm \mathit{23}}^{\mathit{2}}-\mathit{1}\mathit{)}{\mathbf{M}}_{\pm \mathit{23}}+\delta \mathit{[}{\lambda}_{\pm \mathit{13}}^{\mathit{2}}-\mathit{1}\mathit{)}{\mathbf{M}}_{\pm \mathit{13}}$$

(A.7)

where ${\lambda}_{i}^{\mathit{2}}$ terms represent the component of **C**^{COL} in the “*i*-th” direction (i.e. **M**_{1} · **C**^{COL}, etc.) and ${\lambda}_{\pm \mathit{12}}^{\mathit{2}},{\lambda}_{\pm \mathit{23}}^{\mathit{2}},{\lambda}_{\pm \mathit{13}}^{\mathit{2}}$ represent the squared stretch in the secondary fiber directions (i.e. **M**_{+12} · **C**^{COL}, etc.). The material constants *γ _{1}, γ_{2}, γ_{3}*, and δ implement the bimodular feature with the following definitions:

$${\gamma}_{i}\mathit{[}{\lambda}_{i}^{\mathit{2}}\mathit{]=}\{\begin{array}{cc}{\gamma}_{i}& \mathit{\text{for}}\phantom{\rule{thinmathspace}{0ex}}{\lambda}_{i}^{\mathit{2}}\ge \mathit{1}\\ \mathit{0}& \mathit{\text{for}}\phantom{\rule{thinmathspace}{0ex}}{\lambda}_{i}^{\mathit{2}}<\mathit{1}\end{array}\text{,\delta [\lambda +122]={\delta for\lambda +122\u226510for\lambda +1221,etc.}$$

(A.8)

^{13}

Recall (8) for the loading stress *S* for UT with a single constituent:

$$S\mathit{(}t\mathit{)}={}^{e}S\mathit{(}\mathit{0}\mathit{)}+{\displaystyle {\int}_{\mathit{0}}^{t}G\mathit{(}t-\tau \mathit{)}\frac{d{}^{e}S}{d\tau}d\tau .}$$

(B.1)

For a step increase in elastic stress *Δ ^{e}S* at time

$${}^{e}S\mathit{(}t\mathit{)}={}^{e}S\mathit{(}\mathit{0}\mathit{)}+H\mathit{(}t,{t}_{\mathit{\text{step}}}\mathit{)}\Delta {}^{e}S\phantom{\rule{thinmathspace}{0ex}}\text{where}\phantom{\rule{thinmathspace}{0ex}}H\mathit{(}t,{t}_{\mathit{\text{step}}}\mathit{)}=\{\begin{array}{cc}\mathit{0}& \mathit{\text{for}}\phantom{\rule{thinmathspace}{0ex}}t<{t}_{\mathit{\text{step}}}\\ 1& \mathit{\text{for}}\phantom{\rule{thinmathspace}{0ex}}t\ge {t}_{\mathit{\text{step}}}\end{array}$$

(B.2)

so that

$$\frac{d{}^{e}S}{d\tau}\Delta {}^{e}S\delta \mathit{(}\tau ,{t}_{\mathit{\text{step}}}\mathit{)}\phantom{\rule{thinmathspace}{0ex}}\text{where}\phantom{\rule{thinmathspace}{0ex}}\delta \mathit{(}t,{t}_{\mathit{\text{step}}}\mathit{)}=\{\begin{array}{cc}\infty & \mathit{\text{for}}\phantom{\rule{thinmathspace}{0ex}}t={t}_{\mathit{\text{step}}}\\ \mathit{0}& \mathit{\text{for}}\phantom{\rule{thinmathspace}{0ex}}t\ne {t}_{\mathit{\text{step}}}\end{array}$$

(B.3)

where δ(*t,t _{step}*) is the unit impulse function. When

$${\int}_{\mathit{0}}^{t}G\mathit{(}t-\tau \mathit{)}\frac{d{}^{e}S}{d\tau}d\tau}={\displaystyle {\int}_{\mathit{0}}^{{t}_{\mathit{\text{step}}}^{-}}G\mathit{(}t-\tau \mathit{)}\Delta {}^{e}S\delta \mathit{(}\tau ,{t}_{\mathit{\text{step}}}\mathit{)}d\tau}+{\displaystyle {\int}_{{t}_{\mathit{\text{step}}}^{-}}^{t}G\mathit{(}t-\tau \mathit{)}\Delta {}^{e}S\delta \mathit{(}\tau ,{t}_{\mathit{\text{step}}}\mathit{)}d\tau$$

(B.4)

where the first term on the right side of (B.4) is zero because the unit impulse function is zero for *t* ≠ *t _{step}*. Using an expression for

$$\int}_{\mathit{0}}^{t}G\mathit{(}t-\tau \mathit{)}\frac{d{}^{e}S}{d\tau}d\tau}={\displaystyle {\int}_{{t}_{\mathit{\text{step}}}^{-}}^{t}{}^{e}S\delta \mathit{(}\tau ,{t}_{\mathit{\text{step}}}\mathit{)}d\tau}+{\displaystyle {\int}_{{t}_{\mathit{\text{step}}}^{-}}^{t}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}\left(\frac{-\mathit{(}t-\tau \mathit{)}}{{\tau}_{i}}\right)\Delta {}^{e}S\delta \mathit{(}\tau ,{t}_{\mathit{\text{step}}}\mathit{)}d\tau .$$

(B.5)

By noting the integral of the unit impulse function is the Heaviside step function, the first term on the right side of (B.5) becomes:

$${{\displaystyle {\int}_{{t}_{\mathit{\text{step}}}^{-}}^{t}\Delta {}^{e}S\delta \mathit{(}\tau ,{t}_{\mathit{\text{step}}}\mathit{)}d\tau}=\Delta {}^{e}\mathit{\text{SH}}\mathit{(}\tau ,{t}_{\mathit{\text{step}}}\mathit{)}|}_{\tau ={t}_{\mathit{\text{step}}}^{-}}^{\tau =t}=\Delta {}^{e}S.$$

(B.6)

The following equation is needed to integrate the product of a function, *f(t)*, and the impulse function, δ(*t,a*).

$${\int}_{-d}^{d}\delta (t,a)f(t)\mathit{\text{dt}}=f(a)\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{for}}-\infty \le -d<a<d\le \infty .$$

(B.7)

Using (B.7) the second term on the right side of (B.5) becomes

$$\int}_{{t}_{\mathit{\text{step}}}^{-}}^{t}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}\left(\frac{-\mathit{(}t-\tau \mathit{)}}{{\tau}_{i}}\right)\Delta {}^{e}S\delta \mathit{(}\tau ,{t}_{\mathit{\text{step}}}\mathit{)}d\tau =\Delta {}^{e}S{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}}\left(\frac{-\mathit{(}t-{t}_{\mathit{\text{step}}}\mathit{)}}{{\tau}_{i}}\right).$$

(B.8)

Substituting (B.5), (B.6) and (B.8) into (B.1) yields the stress at *t ≥ t _{step}*

$$\mathbf{S}\mathit{(}t\mathit{)}={}^{e}\mathbf{S}\mathit{(}\mathit{0}\mathit{)}+{\Delta}^{e}\mathbf{S}+{\Delta}^{e}\mathbf{S}{\displaystyle \sum _{i=\mathit{1}}^{n}{g}_{i}\phantom{\rule{thinmathspace}{0ex}}\mathit{\text{exp}}}\left(\frac{-\mathit{(}t-{t}_{\mathit{\text{step}}}\mathit{)}}{{\tau}_{i}}\right).$$

(B.9)

GC Thomas and SM Klisch thank Professor Pasquale Vena, Professor Alberto Redaelli, and their colleagues in the Laboratory of Biological Structure Mechanics and the Cellular and Molecular Biomechanics Research Group, respectively, for hosting their research visit at Politecnico di Milano from July 2007 – February 2008.

^{1}These differences are discussed in more detail in Sections 2 and 4.

^{2}These assumptions are discussed further in the Discussion.

^{3}In (4), stresses are VE stresses; a superscript “*e*” will be used to designate elastic stresses.

^{4}Equation numbers preceded by “A” and “B” refer to those in the appendices.

^{5}VE models usually include a term for strain rate (i.e. *d***E**/*dt*) inside the convolution integral to account for the dependence on strain rate. As stress is proportional to strain, the dependence on strain rate can also be achieved with a stress rate term, *d***S**/*dt* [23].

^{6}Note that as *t* → ∞, *S(t)* is the equilibrium (elastic) stress * ^{e}S(0) + Δ^{e}S*.

^{7}Note that relative to a zero stress-zero strain configuration, *R* defined by (19) is strain-independent which is also a characteristic of quasi-linear viscoelasticity.

^{8}The derivation is nearly identical to that in [23] and the reader is referred to that paper for full details.

^{9}This assumption is discussed further in the Discussion.

^{10}Eqn. (25) defines a viscous material constant that is independent of strain, decoupling elastic and viscous properties and is preferable to eqn. (24), which couples elastic and viscous properties, for statistical analysis of COL viscous properties.

^{11}For GD-85 specimens the COL VE parameters $({\tau}_{\mathit{1}}^{\mathit{\text{COL}}},{g}_{\mathit{1}}^{\mathit{\text{COL}}})$ were independent of initial guesses.

^{12}Also, our own pilot simulations using a poroviscoelastic model (results of which are not presented) justified the assumption that fluid flow dependent viscoelasticity can be neglected for out UT protocols.

^{13}It can be seen from (A.8) that when a COL fiber direction is not stretched in tension, it is “turned off” and does not contribute to ^{e}**S**^{BIM}.

1. Linn FC, Sokoloff L. Movement and composition of interstitial fluid of cartilage. Arthritis Rheum. 1965;8:481–494. [PubMed]

2. Buckwalter J, Hunziker E, Rosenberg L, Coutts R, Adams M, Eyre D. Articular cartilage: composition and structure. In: Woo SL-Y, Buckwalter JA, editors. Injury and Repair of the Musculoskeletal Soft Tissues. Park Ridge, IL: American Academy of Orthopaedic Surgeons; 1988. pp. 405–425.

3. Mow VC, Ratcliffe A. Structure and function of articular cartilage and meniscus. In: Mow VC, Hayes WC, editors. Basic Orthopaedic Biomechanics. New York: Raven Press; 1997. pp. 113–178.

4. Klisch SM, Chen SS, Sah RL, Hoger A. A growth mixture theory for cartilage with applications to growth-related experiments on cartilage explants. Journal of Biomechanical Engineering. 2003;125:169–179. [PubMed]

5. Klisch SM, Asanbaeva A, Oungoulian SR, Thonar EJ, Masuda K, Davol A, Sah RL. A cartilage growth mixture model with collagen remodeling: validation protocols. Journal of Biomechanical Engineering. 2008;130:031006:031001–031006:031011. [PMC free article] [PubMed]

6. Schinagl RM, Gurskis D, Chen AC, Sah RL. Depth-dependent confined compression modulus of full-thickness bovine articular cartilage. J Orthop Res. 1997;15:499–506. [PubMed]

7. Huang CY, Stankiewicz A, Ateshian GA, Mow VC. Anisotropy, inhomogeneity, and tension-compression nonlinearity of human glenohumeral cartilage in finite deformation. J Biomech. 2005;38(4):799–809. [PMC free article] [PubMed]

8. Soltz MA, Ateshian GA. A conewise linear elasticity mixture model for the analysis of tension-compression nonlinearity in articular cartilage. J Biomech Eng. 2000;122:576–586. [PMC free article] [PubMed]

9. Laasanen M, Toyras J, Korhonen R, Rieppo J, Saarakkala S, Nieminen M, Hirvonen J, Jurvelin JS. Biomechanical properties of knee articular cartilage. Biorheology. 2003;40:133–140. [PubMed]

10. Donzelli PS, Spilker RL, Ateshian GA, Mow VC. Contact analysis of biphasic transversely isotropic cartilage layers and correlations with tissue failure. J Biomech. 1999;32(10):1037–1047. [PubMed]

11. Krishnan R, Park S, Eckstein F, Ateshian GA. Inhomogeneous cartilage properties enhance superficial insterstitial fluid support and frictional properties, but do not provide a homogeneous state of stress. J Biomech Eng. 2003;125(3):569–577. [PMC free article] [PubMed]

12. Mak AF. The apparent viscoelastic behavior of articular cartilage-the contributions from the intrinsic matrix viscoelasticity and interstitial fluid flows. J Biomech Eng. 1986;108:123–130. [PubMed]

13. Suh J-K, DiSilvestro MR. Biphasic poroviscoelastic behavior of hydrated biological soft tissue. J Appl Mech. 1999;66:528–535.

14. Huang CY, Mow VC, Ateshian GA. The role of flow-independent viscoelasticity in the biphasic tensile and compressive responses of articular cartilage. J Biomech Eng. 2001;123(5):410–417. [PubMed]

15. Korhonen RK, Laasanen MS, Toyras J, Lappalainen R, Helminen HJ, Jurvelin JS. Fibril reinforced poroelastic model predicts specifically mechanical behavior of normal, proteoglycan depleted and collagen degraded articular cartilage. J Biomech. 2003;36(9):1373–1379. [PubMed]

16. Wilson W, van Donkelaar CC, van Rietbergen B, Ito K, Huiskes R. Stresses in the local collagen network of articular cartilage: a poroviscoelastic fibril-reinforced finite element study. J Biomech. 2004;37:357–366. [PubMed]

17. Garcia JJ, Cortes DH. A nonlinear biphasic viscohyperelastic model for articular cartilage. J Biomech. 2006;39(16):2991–2998. [PubMed]

18. Fung YC. Stress-strain history relations of soft tissues in simple elongation. In: Fung YC, Perrone N, Anliker M, editors. Biomechanics: Its Foundations and Objectives. Englewood Cliffs, NJ: Prentice-Hall; 1972. pp. 181–208.

19. Provenzano P, Lakes R, Keenan T, Vanderby R., Jr. Nonlinear ligament viscoelasticity. Ann Biomed Eng. 2001;29(10):908–914. [PubMed]

20. Suh JK, Bai S. Finite element formulation of biphasic poroviscoelastic model for articular cartilage. J Biomech Eng. 1998;120:195–201. [PubMed]

21. DiSilvestro MR, Suh JK. A cross-validation of the biphasic poroviscoelastic model of articular cartilage in unconfined compression, indentation, and confined compression. J Biomech. 2001;34(4):519–525. [PubMed]

22. Park S, Ateshian GA. Dynamic response of immature bovine articular cartilage in tension and compression, and nonlinear viscoelastic modeling of the tensile response. J Biomech Eng. 2006;128(4):623–630. [PMC free article] [PubMed]

23. Vena P, Gastaldi D, Contro R. A constitutent-based model for the non linear viscoelastic behavior of ligaments. Journal of Biomechanical Engineering. 2006;128:449–457. [PubMed]

24. Kempson GE, Muir H, Pollard C, Tuke M. The tensile properties of the cartilage of human femoral condyles related to the content of collagen and glycosaminoglycans. Biochim Biophys Acta. 1973;297:456–472. [PubMed]

25. Kempson GE, Tuke MA, Dingle JT, Barrett AJ, Horsfield PH. The effects of proteolytic enzymes on the mechanical properties of adult human articular cartilage. Biochim Biophys Acta. 1976;428:741–760. [PubMed]

26. Schmidt MB, Mow VC, Chun LE, Eyre DR. Effects of proteoglycan extraction on the tensile behavior of articular cartilage. J Orthop Res. 1990;8:353–363. [PubMed]

27. DiSilvestro MR, Suh JK. Biphasic poroviscoelastic characteristics of proteoglycan-depleted articular cartilage: simulation of degeneration. Ann Biomed Eng. 2002;30(6):792–800. [PubMed]

28. Asanbaeva A, Masuda K, Thonar EJ-MA, Klisch SM, Sah RL. Mechanisms of cartilage growth: modulation of balance between proteoglycan and collagen in vitro using chondroitinase ABC. Arthritis Rheum. 2007;56:188–198. [PubMed]

29. Asanbaeva A, Tam J, Schumacher BL, Klisch SM, Masuda K, Sah RL. Articular cartilage tensile integrity: Modulation by matrix depletion is maturation-dependent. Archives of Biochemistry and Biophysics. 2008;474(1):175–182. [PMC free article] [PubMed]

30. Al Jamal R, Roughley PJ, Ludwig MS. Effect of glycosaminoglycan degradation on lung tissue viscoelasticity. Am J Physiol Lung Cell Mol Physiol. 2001;280(2):L306–L315. [PubMed]

31. Tanaka E, Aoyama J, Tanaka M, Van Eijden T, Sugiyama M, Hanaoka K, Watanabe M, Tanne K. The proteoglycan contents of the temporomandibular joint disc influence its dynamic viscoelastic properties. J Biomed Mater Res A. 2003;65(3):386–392. [PubMed]

32. Elliott DM, Robinson PS, Gimbel JA, Sarver JJ, Abboud JA, Iozzo RV, Soslowsky LJ. Effect of altered matrix proteins on quasilinear viscoelastic properties in transgenic mouse tail tendons. Ann Biomed Eng. 2003;31(5):599–605. [PubMed]

33. Liao J, Vesely I. Relationship between collagen fibrils, glycosaminoglycans, and stress relaxation in mitral valve chordae tendineae. Ann Biomed Eng. 2004;32(7):977–983. [PubMed]

34. Li LP, Herzog W. The role of viscoelasticity of collagen fibers in articular cartilage: theory and numerical formulation. Biorheology. 2004;41(3–4):181–194. [PubMed]

35. Li LP, Herzog W, Korhonen RK, Jurvelin JS. The role of viscoelasticity of collagen fibers in articular cartilage: axial tension versus compression. Med Eng Phys. 2005;27(1):51–57. [PubMed]

36. Wilson W, van Donkelaar CC, van Rietbergen B, Huiskes R. A fibril-reinforced poroviscoelastic swelling model for articular cartilage. J Biomech. 2005;38(6):1195–1204. [PubMed]

37. Garcia JJ, Cortes DH. A biphasic viscohyperelastic fibril-reinforced model for articular cartilage: formulation and comparison with experimental data. J Biomech. 2007;40(8):1737–1744. [PubMed]

38. Basser PJ, Schneiderman R, Bank RA, Wachtel E, Maroudas A. Mechanical properties of the collagen network in human articular cartilage as measured by osmotic stress technique. Archives of Biochemistry and Biophysics. 1998;351:207–219. [PubMed]

39. Oungoulian SR, Chen SS, Davol A, Sah RL, Klisch SM. Extended two-compartmental swelling stress model and isotropic cauchy stress equation for articular cartilage proteoglycans; ASME Summer Bioengineering Conference; Keystone, CO: 2007.

40. Klisch SM. A bimodular polyconvex anisotropic strain energy function for articular cartilage. Journal of Biomechanical Engineering. 2007;129:250–258. [PubMed]

41. Asanbaeva A. Cartilage growth and remodeling: modulation of growth phenotype and tensile integrity. San Diego, La Jolla: PhD, University of California; 2006.

42. Williamson AK, Chen AC, Masuda K, Thonar EJ-MA, Sah RL. Tensile mechanical properties of bovine articular cartilage: variations with growth and relationships to collagen network components. J Orthop Res. 2003;21:872–880. [PubMed]

43. Williamson AK, Masuda K, Thonar EJ-MA, Sah RL. Growth of immature articular cartilage in vitro: correlated variation in tensile biomechanical and collagen network properties. Tissue Eng. 2003;9:625–634. [PubMed]

44. Farndale RW, Buttle DJ, Barrett AJ. Improved quantitation and discrimination of sulphated glycosaminoglycans by use of dimethylmethylene blue. Biochim Biophys Acta. 1986;883:173–177. [PubMed]

45. Woessner JF. The determination of hydroxyproline in tissue and protein samples containing small proportions of this imino acid. Arch Biochem Biophys. 1961;93:440–447. [PubMed]

46. Kim YJ, Sah RLY, Doong JYH, Grodzinsky AJ. Fluorometric assay of DNA in cartilage explants using Hoechst 33258. Anal Biochem. 1988;174:168–176. [PubMed]

47. Sasazaki Y, Shore R, Seedhom BB. Ultrastructure of cartilage under tensile strain. Trans Orthop Res Soc. 2004;29:606.

48. Ficklin TP, Thomas GC, Barthel JC, Asanbaeva A, Thonar EJ, Masuda K, Chen AC, Sah RL, Davol A, Klisch SM. Articular cartilage mechanical and biochemical property relations before and after in vitro growth. Journal of Biomechanics. 2007;40:3607–3614. [PMC free article] [PubMed]

49. Morel V, Quinn TM. Cartilage injury by ramp compression near the gel diffusion rate. J Orthop Res. 2004;22:145–151. [PubMed]

50. Hardingham TE, Muir H, Kwan MK, Lai WM, Mow VC. Viscoelastic properties of proteoglycan solutions with varying proportions present as aggregates. J Orthop Res. 1987;5:36–46. [PubMed]

51. Mow VC, Mak AF, Lai WM, Rosenberg LC, Tang LH. Viscoelastic properties of proteoglycan subunits and aggregates in varying solution concentrations. J Biomech. 1984;17(5):325–338. [PubMed]

52. Spirt AA, Mak AF, Wassell RP. Nonlinear viscoelastic properties of articular cartilage in shear. 1989;7:43–49. [PubMed]

53. Zhu W, Mow VC, Koob TJ, Eyre DR. Viscoelastic shear properties of articular cartilage and the effects of glycosidase treatment. J Orthop Res. 1993;11:771–781. [PubMed]

54. Bhatia A, Vesely I. The effect of glycosaminoglycans and hydration on the viscoelastic properties of aortic valve cusps; Conf Proc IEEE Eng Med Biol Soc; 2005. pp. 2979–2980. [PubMed]

PubMed Central Canada is a service of the Canadian Institutes of Health Research (CIHR) working in partnership with the National Research Council's national science library in cooperation with the National Center for Biotechnology Information at the U.S. National Library of Medicine(NCBI/NLM). It includes content provided to the PubMed Central International archive by participating publishers. |