|Home | About | Journals | Submit | Contact Us | Français|
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  and the remaining weight is a porous solid matrix (SM) composed of glycosaminoglycans (GAGs), collagens (COLs), and other non-collagenous proteins . 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 .
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  and depth from the articular surface , exhibit anisotropy with respect to directions relative to anatomical split-line direction , 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  R, but not τ, is strain-dependent .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,  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  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 , temporomandibular joint disc tissue , tendon tissue , and mitral valve tissue . While most GAG depletion studies with AC involved the use of mature tissue, a recent study  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 t0 and a configuration κ at time t, such that κ0 and κ represent equilibrium and non-equilibrium states, respectively. The right Cauchy Green deformation tensor C is
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
and is related to Cauchy stress T or first Piola-Kirchhoff stress P using
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 (SGAG) and COL (SCOL) stresses at every time t
Full details of the elastic stress equations are presented in the Appendix. Briefly, the elastic GAG stress eSGAG is based on a two compartmental (i.e. extra- and intra-fibrillar water compartments) isotropic swelling stress model originally proposed in  and used in [5, 39] to derive a continuum level GAG stress equation with material constants α1 and α2. The elastic COL stress eSCOL is based on  and contains both an isotropic component eSO with material constant μ and an anisotropic bimodular fiber-reinforced component eSBIM 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, (Figure 1). Each constituent is then prescribed an initial deformation, 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)=eS(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 and κ0, which is equivalent to and from (A.1) eSGAG(0) = −α1I.4 To balance this swelling stress, the COL network must support a tensile prestress produced by an initial COL deformation gradient tensor, , which yields the initial collagen stress tensor, eSCOL(0)=α1I. 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
In this study τi and gi are assumed strain-independent to represent quasi-linear constituent viscoelasticity; however, the constituent based VE model renders nonlinear tissue viscoelasticity with strain-dependent tissue time constant and relaxation ratios, as derived below for a step increase in strain. To simplify the derivation, only the stress S for uniaxial loading with a single constituent will be considered; thus (6) reduces to
Consider the stress response to a step increase in strain and, consequently, a step increase in elastic stress ΔeS at time tstep>0 with initial equilibrium elastic stress eS(0) (Figure 2). For this special case, the following equation for the stress at t ≥ tstep is obtained in Appendix B:
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(tstep) following a step increase in strain and can be calculated from
where the equilibrium stress is
The peak stress is obtained upon evaluating (9) at t=tstep for each constituent
For a single constituent (15) reduces to
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
For a single constituent (18) reduces to
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  is modified in three ways: (1) the initial equilibrium elastic stresses of the constituents, eSGAG(0) and eSCOL(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, ci(t), and a scalar, cCOL, specific to COL are defined as
The time discretization procedure used to calculate SCOL(t + Δt) after a time increment Δt is
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 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 between groups. Instead, an aggregate COL time constant and an aggregate COL relaxation ratio are defined as functions of 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 it takes to relax by 63.2% can be calculated from (15) which reduces upon neglecting intrinsic GAG VE to
where 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
suggesting a definition for an aggregate COL relaxation ratio as
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 , 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 . 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  and analyzed to quantify content of sulfated GAG  and hydroxyproline . 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 Eexp was defined as equilibrium stress divided by applied strain (0.10). The experimental time constant τexp was defined as the amount of time it takes for the stress to relax by 63.2% from the peak stress that occurs at the end of the ramp phase. The experimental relaxation ratio Rexp was defined as equilibrium stress / peak stress.
Experimental stress and biochemical data were used to calculate elastic GAG-COL parameters in elastic constitutive equations using MATLAB. GAG material constants α1 and α2 (A.1) were calculated using an extended two compartment (i.e. extra- and intra-fibrillar water) swelling stress model for AC PGs developed by  and extended by [5, 39]. For each specimen, biochemical mass measurements were used to calculate GAG swelling stress (normalized to extrafibrillar water area), which is then converted to GAG Cauchy stress (normalized to tissue area). It was assumed that one set of parameters is a valid description for the GAG stress-density relations for each group; therefore, the four specimens for each group were pooled together to curve fit the theoretical GAG Cauchy stress, to yield parameters and α2. The parameter α1 was calculated using , and the continuity equation and equation (3) was used to obtain eSGAG (A.1).
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  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 . 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, , 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 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 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 ; 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 and relaxation ratio 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 (α1,α2), averaged COL elastic parameters (γ1) were adjusted to the GD-55 and GD-85 values while the COL viscous parameters were kept equal to the GD-0 values. Third, the incremental effect of direct GAG interactions on COL viscous properties was isolated: in addition to (α1,α2,γ1), the COL viscous parameters were adjusted resulting in simulations using averaged values for all elastic and viscous parameters for the GD-55 and GD-85 groups.
The effects of GAG depletion treatment on experimental tissue VE parameters (Eexp, τexp, Rexp), aggregate VE parameters s , and COL elastic parameter (γ1) were assessed using analysis of variance (ANOVA) with post-hoc Tukey testing. Since a lesser number of COL viscous parameters were needed to obtain unique solutions for the GD-85 group, they were not directly compared to the GD-0 and GD-55 group values. Consequently, the effects of GAG depletion on COL viscous parameters were compared using paired t-tests for the GD-0 and GD-55 groups only. Correlations between all VE parameters with GAG and COL contents were investigated using linear regression with a t-test analysis of the regression slope. For all statistical analyses, a 0.05 probability of type-1 error was assumed (p = 0.05).
Experimental tissue VE parameters (Eexp, τexp, Rexp) 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 (Eexp, τexp, Rexp) were significantly correlated with GAG content (p<0.05, 0.01, and 0.0001, respectively) but not COL content, although positive trends with COL existed for Eexp (R2=0.17, p=0.19) and Rexp (R2=0.21, p=0.13) (Figure 4).
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 decreased due to increasing GAG depletion treatment (p<0.01) (Table 1; Figure 6), while a significant difference was almost detected for (p=0.06) (Table 1; Figure 6). The aggregate VE relaxation parameter increased due to increasing GAG depletion treatment (p<0.0001) while a significant difference was almost detected for (p=0.08) (Figure 6). The VE parameters 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 γ1 (R2=0.19, p=0.15) and (R2=0.25, p=0.10) (Figure 7).
The parameter studies using averaged GAG elastic parameters (α1,α2) for the GD-55 and GD-85 groups with COL VE parameters 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 (α1,α2,γ1) for the GD-55 and GD-85 groups with COL viscous parameters for the GD-0 group showed that the additional effect of direct GAG interactions on COL elasticity further modulates peak and equilibrium stresses but not relaxation behavior (Figure 8). The simulations with all GAG and COL VE parameters for the GD-55 and GD-85 groups showed that COL VE parameters modulate relaxation behavior in addition to peak and equilibrium stresses.
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 (Eexp, τexp, Rexp) 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, Eexp/τexp values were increased/decreased following GAG depletion and correlated with GAG content in this study; in contrast, Eexp/τexp values were maintained/decreased following GAG depletion in  with more mature AC tissue (18 month old bovine), and Eexp values were decreased following GAG depletion and not correlated with GAG content in [24, 25] with adult human tissue. These differences may be due to the use of a more mature AC tissue in [24–26] than the newborn (~1–3 week old) bovine AC tissue was used here. This hypothesis is supported by recent results that the effects of GAG depletion on tensile properties (i.e. equilibrium modulus, ramp modulus, strength, and failure strain) depend on AC tissue maturational stage .
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 . 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 , 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 , 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 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 . 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 Eexp (p=0.19), Rexp (p=0.13), γ1 (p=0.15) and (p=0.10) it is likely that increasing the number of specimens would allow for the detection of correlations with COL content.
A previous study  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  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  where an isotropic GAG matrix is reinforced with a COL fiber constituent. The elastic GAG stress eSGAG is based on a two compartmental (i.e. extra- and intra-fibrillar water compartments) isotropic swelling stress model for AC GAGs originally proposed in  and used in [5, 39] to derive the polyconvex continuum level model
where α1 and α2 are material constants.
The elastic COL stress eSCOL is based on  and contains both an isotropic component eSO and a anisotropic bimodular component eSBIM:
eSO is defined as
where μ is a material constant and −1 denotes the inverse operator. The fiber reinforcement of the SM appears in eSBIM. In this model, there are three primary fiber directions along mutually orthogonal unit vectors E1, E2, E3 which generate structural tensors
where is the dyadic product. Also, there are two secondary fiber directions in each of the three planes created by the primary unit vectors E1, E2, E3 (planes 1–2, 2–3, 1–3) along unit vectors (E±12, E±23, E±13) defined by angles (Ф±12, Ф±23, Ф±13) as
These generate structural tensors
The strain energy function WBIM is modified from  to model a linear stress-strain response typical of newborn bovine tissue. The resulting expression for eSBIM is
where terms represent the component of CCOL in the “i-th” direction (i.e. M1 · CCOL, etc.) and represent the squared stretch in the secondary fiber directions (i.e. M+12 · CCOL, etc.). The material constants γ1, γ2, γ3, and δ implement the bimodular feature with the following definitions:
Recall (8) for the loading stress S for UT with a single constituent:
For a step increase in elastic stress ΔeS at time tstep>0 with initial equilibrium elastic stress eS(0), the elastic stress response is represented with the Heaviside step function, H(t,tstep)
where δ(t,tstep) is the unit impulse function. When t>tstep, the convolution integral in (B.1) can be evaluated by separating the integration limits at , defined as the moment before the instantaneous step,
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:
The following equation is needed to integrate the product of a function, f(t), and the impulse function, δ(t,a).
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.
1These differences are discussed in more detail in Sections 2 and 4.
2These assumptions are discussed further in the Discussion.
3In (4), stresses are VE stresses; a superscript “e” will be used to designate elastic stresses.
4Equation numbers preceded by “A” and “B” refer to those in the appendices.
5VE models usually include a term for strain rate (i.e. dE/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, dS/dt .
6Note that as t → ∞, S(t) is the equilibrium (elastic) stress eS(0) + ΔeS.
7Note 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.
8The derivation is nearly identical to that in  and the reader is referred to that paper for full details.
9This assumption is discussed further in the Discussion.
10Eqn. (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.
11For GD-85 specimens the COL VE parameters were independent of initial guesses.
12Also, 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.
13It 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 eSBIM.