Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Biomech Eng. Author manuscript; available in PMC 2010 June 1.
Published in final edited form as:
PMCID: PMC2842192

Modeling the Matrix of Articular Cartilage Using a Continuous Fiber Angular Distribution Predicts Many Observed Phenomena



Cartilage is a hydrated soft tissue whose solid matrix consists of negatively charged proteoglycans enmeshed within a fibrillar collagen network. Though many aspects of cartilage mechanics are well understood today, most notably in the context of porous media mechanics, there remain a number of responses observed experimentally whose prediction from theory has been challenging.

Method of approach

In this study the solid matrix of cartilage is modeled with a continuous fiber angular distribution, where fibers can only sustain tension, swelled by the osmotic pressure of a proteoglycan ground matrix.


It is shown that this representation of cartilage can predict a number of observed phenomena in relation to the tissue’s equilibrium response to mechanical and osmotic loading, when flow-dependent and flow-independent viscoelastic effects have subsided. In particular, this model can predict the transition of Poisson’s ratio from very low values in compression (~0.02) to very high values in tension (~2.0). Most of these phenomena cannot be explained when using only three orthogonal fiber bundles to describe the tissue matrix, a common modeling assumption used to date.


The main picture emerging from this analysis is that the anisotropy of the fibrillar matrix of articular cartilage is intimately dependent on the mechanism of tensed fiber recruitment, in the manner suggested by our recent theoretical study (G. A. Ateshian. J Biomech Eng, 129(2):240-9, 2007).

1. Introduction

It is well accepted that the tensile properties of articular cartilage arise primarily from its fibrillar type II collagen matrix [1], whereas its compressive properties are dominated by the osmotic contribution of proteoglycans [2]. A number of theoretical frameworks embodying the disparity between tensile and compressive properties of cartilage have been proposed, accounting for the collagen fibers implicitly [35] or explicitly [69]. These models generally propose discrete fiber families to describe the collagen matrix. They are able to capture the most salient features of the cartilage mechanical response, namely, the tension-compression nonlinearity of the stress-strain curve, and the enhanced interstitial fluid pressurization resulting from this disparity [10, 11].

In a recent theoretical study [12], we examined the anisotropy of fibrous tissues in relation to the distribution of tensed and buckled fibers. We found that tissue anisotropy was dependent on the state of strain (even in the range of small strains), as well as the anisotropy of the fiber angular distribution, as a consequence of strain-dependent fiber recruitment. For example, when the strain tensor exhibits a mix of tensile and compressive principal normal values, only a subset of fibers become loaded in tension, and only those fibers contribute to the anisotropy. A significant outcome of this study was the observation that discrete fiber families represent only special cases of material anisotropy; a more general framework that accommodates continuous fiber angular distributions can predict more complex behaviors. Furthermore, we observed that phenomenological constitutive relations that do not model fibers explicitly cannot reproduce the behavior of continuous fiber distribution models. Therefore, the objective of the current study is to examine whether using a continuous fiber distribution, rather than discrete fiber bundles, yields better agreement with experimental results reported in the cartilage mechanics literature.

Many experimental findings in articular cartilage have been successfully predicted with suitable theoretical models. For example, the linear isotropic biphasic theory for soft hydrated tissues [13] can predict the transient response of cartilage deformation and interstitial fluid pressurization in confined compression creep, stress-relaxation, and dynamic loading [1316], and indentation creep [17]. When accounting for the tension-compression nonlinearity arising from the fibrillar solid matrix, the biphasic theory can also predict the response of cartilage deformation and interstitial fluid pressurization in unconfined compression [3, 4, 7, 11, 18, 19]. When accounting for the intrinsic viscoelasticity of the solid matrix, agreement with experimental results in unconfined compression [19, 20] and uniaxial tension [21] also improve.

However, there exist a number of experimental observations that pertain mostly to the tensile response of cartilage, which have yet to be predicted successfully with theoretical models. First, while Poisson’s ratio of articular cartilage can be as low as 0.02 (and always less than 0.5) in unconfined compression [22, 23], its value increases under free swelling of the tissue [24] and can be as high as 4.0 in uniaxial tension [25, 26]. A constitutive model that can predict this large variation from compression to tension has not yet been reported.

Second, the variation in the anisotropic response of cartilage in the transitional range from compression to tension, when measured under different bathing solution osmolarities [24], has not been predicted either.

Finally, another intriguing observation is the seemingly inconsistent response of cartilage to variations in bath osmolarity under uniaxial tension at various applied strains [27]: Under small tensile strains, when the bathing solution is switched from distilled water to normal saline under isometric conditions, the reaction force at the clamped ends of the specimen increases; under larger strains however, the same alteration in bath osmolarity produces the opposite effect.

The objective of this study is to demonstrate that these behaviors can be predicted from theory when the solid matrix is described with a continuous fiber angular distribution, where the fibers can support tensile loading only. For comparison purposes, it is also shown that a model of the solid matrix that incorporates only three discrete, orthogonal fiber bundles is unable to predict most of these observed responses. This model is selected for comparison because it was the first to explicitly address the disparity between the tensile and compressive properties of articular cartilage [7]; it represented a paradigm shift in the modeling of cartilage, as it could predict the response of this tissue to the canonical problem of unconfined compression, a challenge not met previously [28]. Since then, more elaborate models of the fibrillar matrix of cartilage have been proposed, which employ a larger number of discrete fiber bundles [8, 19, 29], though their predictive ability remains to be exhaustively investigated.

2. Methods

This study focuses on the equilibrium experimental responses reported in the literature, where interstitial fluid flow and intrinsic solid matrix viscoelasticity have subsided. Since these experiments were conducted under various bathing solution osmolarities, we use the framework of the triphasic theory for charged hydrated tissues [30] specialized to equilibrium conditions, under finite deformation. The Cauchy stress tensor T representing the total stress in the mixture is given by


where p is the interstitial fluid pressure, I is the identity tensor, and Te is the stress in the solid matrix. From the conservation of momentum for the mixture, under quasi-static conditions, the total stress must satify


2.1. Osmotic Swelling

Under equilibrium conditions the fluid pressure has contributions from osmotic effects only. As reviewed by Ehrlich et al. [31] and Chahine et al. [32], this osmotic pressure arises primarily from Donnan and entropic effects. Donnan osmotic pressure results from the presence of negatively charged proteoglycans within cartilage, which produce an imbalance in electrolyte concentrations between the tissue and external bath [33]; this Donnan pressure, which varies with the external bath salt concentration, has been traditionally incorporated in triphasic models of articular cartilage [30]. The entropic pressure arises from the configurational entropy of proteoglycan macromolecules within the cartilage; when measuring the osmotic pressure of chondroitin sulfate solutions, the entropic contribution becomes apparent under high salt concentrations, when Donnan effects have subsided due to charge shielding [31, 32]. This entropic osmotic pressure has been observed to increase quadratically with chondroitin sulfate concentration.

In the current analysis, we assume that the entropic osmotic pressure of proteoglycans trapped within the collagen matrix of cartilage is attenuated relative to that in free solution. This assumption allows us to model the total osmotic pressure as the sum of the Donnan pressure arising from the total charge density cF of proteoglycans in cartilage, and the entropic pressure attenuated by the fraction ϕ of cF, where 0 ≤ ϕ ≤1,


The first term on the right-hand-side represents the Donnan contribution under ideal conditions, in a binary solution of monovalent electrolytes such as NaCl; c* is the external bath osmolarity (c*= 2c* where c* is the salt concentration), p* is the ambient pressure in the external bath, R is the universal gas constant, and θ is the absolute temperature. The remaining polynomial terms represent the entropic contribution, where c1 and c2 are coefficients of the corresponding virial expansion, obtained from experimental measurements [31, 32]. Since the proteoglycans are assumed to be trapped within the solid matrix (at least for the duration of typical experiments), the conservation of mass constrains cF according to


where J= det F is the relative volume of the solid matrix, F is the deformation gradient, and crF and ϕrw are the proteoglycan charge density and water volume fraction of the tissue, respectively, in the strain-free reference configuration (when J= 1). In the absence of an entropic contribution (ϕ= 0), the reference configuration can be achieved under traction-free boundary conditions, in the hypertonic reference state (c*→ ∞ and pp*). However, when entropic effects are present, there is no simple method to achieve the true reference configuration experimentally (p≠ 0, ∀ c*).1

Since the osmotic pressure difference π = pp* depends on the solid matrix deformation via the functional dependence p(J), there is an associated osmotic modulus that describes the resulting tissue stiffness, as presented in our recent study [34],2


Π is a fourth-order spatial elasticity tensor describing an isotropic material, as can be deduced from the general form of Eq.(5).3 This expression may be viewed as the modulus of the ‘ground matrix’ of cartilage, but note that this ground matrix cannot exist on its own, since a tensile load bearing solid matrix must exist to balance the swelling pressure under traction-free conditions. This osmotic modulus contributes to the total elasticity tensor An external file that holds a picture, illustration, etc.
Object name is nihms166611ig1.jpg= Π+An external file that holds a picture, illustration, etc.
Object name is nihms166611ig2.jpg of the tissue, where An external file that holds a picture, illustration, etc.
Object name is nihms166611ig2.jpg is the contribution from the constitutive relation for Te, given below. The spatial elasticity tensor is needed in the finite element implementation of this framework [36].

2.2. Fibrous Matrix

The solid matrix is described using a continuous fiber distribution as proposed by Lanir [37, 38] and Humphrey and Yin [39], constrained such that fibers can only support tensile loading [12]. To formulate this general relation, let the strain be represented by the right Cauchy-Green tensor C = FT · F. For fibers that, in the reference configuration, are directed along the spherical angles (θ, [var phi]) in the global orthonormal basis {ei}, the unit vector along this direction is


where 0 ≤ θ< 2π and 0 ≤ [var phi]< π. Let λn represent the stretch along nr, then


where Nr= nr [multiply sign in circle] nr represents the texture tensor along nr. For convenience, the square of the stretch is denoted by In. The fibers along this direction are in tension when In >1.

The Cauchy stress tensor and spatial elasticity tensor in the solid matrix are given by



where Tne(nr) and Cne(nr) represent the Cauchy stress and spatial elasticity tensor of those fibers initially oriented along nr; H(·) is the unit step function that enforces the tension-only contribution. According to standard relations of finite elasticity, these are related to a fiber strain energy density function Ψ via


2.3. Fiber Constitutive Relation

Under most general conditions, it is possible to model fibers initially oriented along nr as transversely isotropic, with nr representing the axis of symmetry. However, in the current treatment, we seek the simplest possible constitutive relation that may capture the observed experimental response in cartilage. Thus it is proposed to let Ψ= Ψ(nr; In) such that Eqs.(10)(11) reduce to



where N=nn=In1F·Nr·FT, and n=In1/2F·nr represents the current direction of the fibers originally oriented along nr.

In particular, based on earlier experimental results [24], we propose to use a power-law relation


where ξ(nr) (units of force per area) and α(nr) (unitless) are material properties that may vary along the directions nr. (For inhomogeneous materials, it is also possible to vary ξ and α spatially.) The exponent is restricted to α ≥ 2 to prevent singularities in the modulus Cne at the strain origin (J= In = 1). When α > 2, the modulus Cne is zero at the strain origin, yielding a smooth transition for the fiber response Tne from compression to tension, whereas α= 2 produces a jump in modulus at the strain origin (producing a piecewise-linear stress-strain response across the origin).

2.4. Fiber Angular Distribution

For a random fiber distribution, which has an equal fraction of fibers pointing in every direction, the material coefficients ξ and α are the same for all nr. In a spherical coordinate system, this isotropic fiber angular distribution may be represented by a sphere.

To represent more general continuous fiber angular distributions in this study, we propose to use an ellipsoidal distribution. ξ(nr) and α(nr) are formulated in a local orthonormal basis {ai} defined by the local texture relative to the global basis {ei}. In this local basis, ξ(nr) and α(nr) can be expressed in an associated spherical coordinate system (Θ, Φ) as





Here, ξi and αi (i= 1− 3) represent the semi-axes of the respective ellipsoids of ξ and α, along ai. To use Eqs.(15)(16) in the calculation of stresses, a standard coordinate transformation needs to be performed between the local fiber distribution basis {ai} and the global basis {ei}, to relate (Θ, Φ) and (θ, [var phi]) for a given nr [equivalent to equating Eqs.(6) and (17)].

To model discrete fiber bundles instead, consider the case of m bundles oriented along the local basis {ai}. In the local spherical coordinate system, the orientation of each of these bundles may be described by a (Θj, Φj) pair representing the directions njr(j=1m). Thus, m material coefficients ξj and αj need to be specified, and the integrals in Eqs.(8)(9) reduce to the summations



where ΔAj represents that fraction of fibers oriented along njr (a user-specified parameter). For example, for equal fiber fractions in all directions, ΔAj = 4π/m.4

2.5. Implementation

A custom-written three-dimensional finite element analysis that can accommodate tissue swelling under finite deformation [34] was modified to incorporate the osmotic pressure and fibrous tissue constitutive relations summarized above. The double integrals appearing in Eqs.(8)(9) were evaluated numerically using integration over a geodesic dome representation of a unit sphere. This approach was adopted because the more traditional alternative of incrementing the spherical angles θ and [var phi] by uniform increments (or by Gaussian quadrature) over 0 ≤ θ< 2π and 0 ≤ [var phi]< π yields highly variable elemental areas dA = sin [var phi] d[var phi], compromising numerical results (such as producing anisotropic results where isotropy is expected, especially at the poles). Given that the representation of the unit sphere consists of n spherical triangles whose vertices coincide with those of the geodesic dome, the double integration reduces to


where ΔAi is the surface area of the i – th spherical triangle and (θi, [var phi]i) are the spherical coordinates of its centroid. For example, the geodesic dome can consist of an icosahedron with order m regular tesselation, which can be constructed in the symbolic software program Mathematica® (Wolfram Research Inc., version 6) using the command Geodesate[“Icosahedron”,m, {0,0,0},1]. (Additional processing is required to find (θi, [var phi]i) and ΔAi.) A value of m= 4 was used in the current implementation, yielding n= 320 spherical triangles.

In experimental studies, the applied strain is not the same as the true strain in the solid matrix, the latter being represented by C in the above equations, or equivalently, by the Lagrangian strain E= (CI)/2; the applied strain Ea and the true strain E differ by the free-swelling strain E0 in the traction-free initial configuration, Ea= EE0.5 The results presented below employ the applied strain Ea, for consistency with experimental studies.

2.6. Selection of Material Properties

As shown in section 3, the predictions of the model described above are compared to various experimental studies in the literature. The constitutive model presented in this study includes five material coefficients to describe the osmotic response resulting from the proteoglycans, crF,ϕrw for the Donnan effect, and c1, c2, ϕ for the configurational entropy; and six material coefficients to describe the fibrillar matrix, ξi and αi(i= 1−3).

In all cases reported here, the virial coefficients appearing in Eq.(3), which represent the entropic contribution to the osmotic pressure, are taken to be c1= 0.381 × 10−3MPa·m3/Eq and c2= 0.241 × 10−6MPa·m6/Eq2, based on the earlier study of Chahine et al. [32]. The temperature is assumed to be θ= 298K and the universal gas constant is given by R= 8.314J/mol·K. Furthermore, in the absence of definitive experimental evidence, the value of ϕ is taken to be at the midpoint of its range (0 ≤ ϕ ≤ 1).

The water content ϕrw of articular cartilage typically varies from 0.7 to 0.9, depending on animal species, joint, age, and location through the dept of the articular layer. Since variations over that range have relatively negligible influence on the predicted response, the value of ϕrw=0.8 was adopted throughout. The fixed charge density crF of articular cartilage typically ranges from 0 to 200mEq/L [24, 40]; however, since the ideal Donnan behavior adopted in this study overestimates the actual osmotic pressure [31, 32, 41], and since additional contributions from entropic effects are incorporated into the model, representative values of crF are taken to vary over the range 0crF120mEq/L, depending on the experimental study against which the comparison is performed. Using these somewhat lower values of crF produce osmotic predictions more consistent with experimental evidence.

Since the constitutive model for a continuous ellipsoidal fiber angular distribution has not been used previously to model cartilage, there are no existing guidelines for selecting ξi and αi (i= 1− 3). Therefore, representative values of these parameters are proposed for the first time in this study, based on an iterative manual procedure whose objective is to reproduce literature data with reasonable qualitative agreement, as shown in the next section. This iterative procedure is guided by parametric analyses of variations in ξi and αi as outlined in the Appendix. In finite element analyses that specifically explore the anisotropic response of the tissue along two perpendicular directions, different values are adopted for these parameters along two of the three directions. Otherwise, it is assumed that ξ1= ξ2= ξ3 and α1= α2= α3 (representing an isotropic fiber angular distribution), for simplicity.

2.7. Sample Inhomogeneity, Geometry, and Texture Orientation

All finite element analyses reported in section 3 assume that the material properties are homogeneous throughout the cartilage sample. This assumption is based on the fact that comparisons are typically performed against experimental studies where samples were harvested from specific zones of the articular layer; within each zone, it is reasonable to assume that variations in properties are relatively small.

Model predictions are compared to experiments where cartilage samples tested in compression are cylindrical or cubic; and those tested in tension are either dumbbell-shaped or rectangular prismatic bars. In these experimental studies, the gauge section of the specimen (the region over which strain is measured or estimated) is assumed to exhibit a homogeneous strain field. Therefore, to simplify the finite element analyses performed in this study, the finite element model geometry is assumed to be a rectangular prismatic bar, and idealized boundary conditions are applied to neglect end effects (such as clamping of the ends of a rectangular strip). Consequently, the stress and strain fields predicted by the finite element analyses are homogeneous, and the precise dimensions of the tissue sample do not matter. The faces of the rectangular prism align with the global coordinate basis {ei}.

In experimental studies, cartilage samples are usually harvested such that one of the faces of the sample is parallel to the articular surface. Furthermore, in the case of cubic, rectangular, or dumbbell-shaped specimens, the remaining faces are typically oriented parallel and perpendicular to the local split-line direction. Therefore, in the finite element model of the tissue sample, the texture directions {ai} of the fiber angular distribution, described in section 2.4, are taken to be aligned with the global coordinate basis {ei}, where i= 1 denotes the direction parallel to the local split-line; i= 2 denotes the direction perpendicular to the local split-line, in the plane tangent to the articular surface; and i= 3 is the direction normal to the articular surface.

3. Predictions of Experimental Observations

3.1. Stress-Strain Response in Compression and Tension

It is well known that cartilage exhibits a significantly higher modulus in tension [1] than compression [42]. The tensile response has been attributed primarily to the ability of collagen fibers to resist tension [43] whereas the compressive response has been attributed to the negatively charged proteoglycans [2]. The disparity in the tensile and compressive responses was illustrated, for example, in the study of Huang et al. [26], who tested adjacent samples of human glenohumeral joint cartilage in tension and compression (Figure 1a). The tensile moduli were observed to be up to two orders of magnitude higher than the compressive moduli. The tissue was stiffer in tension when loaded parallel to the local split-line direction, compared to the perpendicular direction. It was also stiffer in tension in the superficial zone than in the middle zone, while the opposite was true in compression.

Figure 1Figure 1
(a) Representative stress-strain responses of human glenohumeral joint cartilage in tension and compression, from the superficial zone (SZ) and middle zone (MZ). Tensile responses are reported for cartilage strips harvested parallel and perpendicular ...

These representative results can be simulated with the continuous, ellipsoidal fiber distribution as shown in Figure 1b. In this theoretical prediction, the finite element model of the cartilage sample is loaded with prescribed displacements from compression to tension, in a 0.15 M NaCl bath (c*= 300mol/m3).

The theoretical prediction, which uses different material properties in the superficial zone ( crF=40Eq/m3,ϕrw=0.8, ϕ= 0.5, ξ1= 2.8MPa, ξ2= ξ3= 1.6MPa, α1= α2= α3= 2.5) and middle zone ( crF=80Eq/m3,ϕrw=0.8, ϕ= 0.5, ξ1= 4MPa, ξ2= ξ3= 2MPa, α1= α2= α3= 3.5) is able to capture the rapid rise in the stress-strain response when transitioning from compression to tension. By performing separate loading analyses along direction 1 and direction 2 on each finite element model within each zone (superficial and middle), it is found that the model can also capture the anisotropic response in tension, as well as the inhomogeneous response through the depth. No differences are observed in the compressive response to loading parallel and perpendicular to the split line direction; this outcome is not surprising since the ‘ground matrix’ is isotropic according to Eq.(5), and the anisotropic fibers can only sustain tension.

Qualitatively similar stress-strain responses can be predicted when using three orthogonal fiber bundles with the same material properties (results not shown). Thus, the benefits of using a continuous fiber distribution over a discrete orthogonal distribution is not immediately apparent from the stress-strain response alone.

3.2. Poisson’s Ratio in Compression and Tension

The equilibrium Poisson’s ratio of articular cartilage has been measured in compression using optical techniques [20, 2224, 44], and reported to be approximately in the range 0.01–0.25 in bovine articular cartilage; it increases significantly from the superficial zone to the deep zone [23]. In contrast, studies using optical measurements of Poisson’s ratio of human articular cartilage in tension [25, 26] report much higher values, in the approximate range 1.0–4.0 for human articular cartilage from the superficial zone, and 0.3–1.0 for middle zone cartilage. For human glenohumeral joint cartilage, Huang et al. [26] also reported that there were no differences in the tensile Poisson’s ratio in the plane tangent to the articular surface, when loading along or perpendicular to the split-line direction. Wang et al. [23] and Chahine et al. [24] showed that Poisson’s ratio is higher at smaller compressive strains, and decreases as the applied compressive strain is increased in magnitude. In contrast, in tension, Elliott et al. [25] and Huang et al. [26] have reported that Poisson’s ratio is nearly constant. These experimental results beg the question as to how Poisson’s ratio may vary over such a large range (as low as 0.01 to as high as 4.0) in the transition from compression to tension.

Looking at the theoretical cases examined in the previous section (Figure 1), the predicted incremental Poisson’s ratio6 for the continuous fiber angular distribution is shown in Figure 2a. Remarkably, these predictions are able to capture the salient features reported in experimental studies: Poisson’s ratio is very small in compression (as low as 0.017 in this example), rises in the range of small compressive strains while remaining below 0.5, then rapidly increases to very high values (as high as 3.1 here) in the range of tensile strains. In compression, Poisson’s ratio increases in magnitude from the superficial to the middle zone, whereas its tensile counterpart decreases with depth from the articular surface. Differences in Poisson’s ratio measured parallel versus perpendicular to the split-line direction (ν12 and ν21) are generally small.

Figure 2Figure 2
Theoretical prediction of Poisson’s ratio for in uniaxial tension and compression using (a) a continuous fiber distribution, or (b) three orthogonal fiber bundles. SZ: superficial zone; MZ: middle zone.

When modeling the solid matrix with three orthogonal fiber bundles, it is found that the predictions for Poisson’s ratio (Figure 2b) are unable to reproduce the experimental results; in this case, Poisson’s ratio is found to decrease monotonically from compression to tension, with values ranging from 0.22 down to 0.002.

One feature of these theoretical predictions which may initially appear somewhat inconsistent with experimental results is the near-linear decrease in the tensile Poisson’s ratio in the range of higher tensile strains (as seen for the superficial zone response in Figure 2a), whereas experimental studies have suggested that it should be approximately constant. However, a more careful analysis of experimental data may reveal that this behavior is a subtle feature, potentially overlooked in experimental studies. For example, looking at the representative raw data of Figure 2 in the study of Elliott et al. [25], and replotting it as E11versus E22, a subtle quadratic response emerges which, when differentiated to yield ν12, produces a tensile Poisson’s ratio that linearly decreases from 1.9 to 0.6, as the tensile strain increases by 17% over the tare configuration (Figure 3). Thus, it is plausible that the theoretical predictions for ν12 in the range of large tensile strains (Figure 2a) are in fact consistent with experimental results, as may be further explored in future studies.

Figure 3
Representative experimental tensile response data from Figure 2 of Elliott et al. [25], converted from infinitesimal strain to Lagrangian strain using the formula given by these authors, and re-plotted as lateral versus axial strain. The solid curve is ...

3.3. Effect of Bath Osmolarity on Compressive Response

In the compressive loading experiments of Chahine et al. [24] on immature bovine articular cartilage, the bath NaCl concentration was increased from 0.015M to 2M, producing a softer stress-strain response and lower Poisson’s ratio, as shown in Figure 4a. These effects can be similarly predicted from theory, when using a continuous fiber distribution, as shown in Figure 5. These simulations use the representative material properties crF=80Eq/m3,ϕrw=0.8, ϕ= 0.5, ξ1= ξ2= ξ3= 1.6MPa, α1= α2= α3= 2.5.

Figure 4Figure 4
Average experimental response from the study of Chahine et al. [24] on immature bovine cartilage at various NaCl bathing solution concentration, for compressive loading along the split-line direction, and lateral contraction perpendicular to it. (a) Stress ...
Figure 5Figure 5
Theoretical predictions for testing conditions similar to the study of Chahine et al. [24], using a continuous fiber angular distribution.

When modeling the solid matrix with three orthogonal fiber bundles, using the same material properties, the variation of the stress-strain response with increasing salt concentration exhibits the same qualitative variation as experimental results. However, as reported above, Poisson’s ratio does not increase in the range of small compressive strains (Figure 2b). Furthermore, the responses at 0.015M and 0.15M are nearly identical (results not shown), contrary to experimental studies.

3.4. Effect of Proteoglycan Digestion on Tensile Response

Schmidt et al. [45] reported on the tensile response of bovine knee cartilage (18 months old) before and after proteoglycan extraction. In a result that may seem counter-intuitive, they observed that the stress-strain response of PG-depleted samples (loaded at 0.05cm/min ) exhibited a stiffer response for the digested tissue, in the range of higher tensile strains (Figure 6a). They attributed this observation to “…the collapse of the inflated collagen network in specimens depleted of proteoglycan glycosaminoglycans, resulting in a collagen network that is more aligned along the axis of the specimens prior to tensile testing.”

Figure 6Figure 6
(a) Stress-strain response of normal and proteoglycan-depleted bovine articular cartilage in uniaxial tension, adapted from Figure 6 of Schmidt et al. [45], with permission. (b) Theoretical simulation of similar testing conditions using a continuous fiber ...

To examine the effect of PG depletion from theory, a uniaxial tensile loading analysis was simulated for specimens with a normal ( crF=120mEq/l) and hypothetical depleted level ( crF=12mEq/l) of PG. The remaining material properties were set to ϕrw=0.8, ϕ= 0.5, ξ= 0.2MPa and α= 2.5 (isotropic fiber angular distribution, for simplicity). The equilibrium stress-strain results for these two cases (Figure 6b) show remarkable qualitative agreement with the representative experimental results of Schmidt et al. [45] (Figure 6a), confirming that the PG-depleted tissue first exhibits a softer response, then a stiffer response as strain increases. In contrast, an orthogonal fiber distribution, with the same material properties, shows virtually no change in the tensile stress-strain response before and after PG-depletion (results not shown).

3.5. Effect of Bath Osmolarity on Tensile Response

Grodzinsky et al. [46] tested the response of bovine femoropatellar articular cartilage (1.5–2 yrs-old) under isometric tensile strain (15–20 percent), in baths of varying NaCl concentration (from 0 to 1M). As the salt concentration increased, they observed that the equilibrium clamping stress decreased in magnitude (Figure 7a). In a subsequent, related study, Akizuki et al. [27] showed that the variation in the equilibrium stress with increasing salt concentration (from 0 to 0.15 M NaCl) depended on the magnitude of the isometric clamping strain: At low strains, the clamping stress increased in magnitude with increasing salt concentration; however, as the clamping strain was raised, the clamping stress decreased in magnitude with increasing salt concentration, consistent with the earlier study.

Figure 7Figure 7
(a) Equilibrium uniaxial tensile stress response of bovine articular cartilage in bathing solutions of various NaCl concentrations, under 20% isometric strain; adapted from Figure 8 of Grodzinsky et al. [46], with permission. (b) Theoretical predictions ...

A theoretical analysis of these experimental configurations, using a random continuous fiber angular distribution, reveals a complex interaction of swelling and fiber responses. Using crF=80mEq/l,ϕrw=0.8, ϕ= 0.5, ξ= 2MPa, and α= 2.5, and varying the external NaCl bath concentration from 10−4M to 1M, it is observed from theory that the variation in clamping stress with salt concentration depends very significantly on the clamping strain (Figure 7b). The theoretical prediction at a clamping strain of 7% is qualitatively consistent with the results of Grodzinsky et al. [46], showing a similar sigmoidal rise in stress with decreasing NaCl concentration. Focusing on the clamping stresses in the isotonic state at 0.15 M, relative to the hypotonic state, it is noted that the clamping stress rises in the isotonic bath when the clamping strain is less than 2.1%, and decreases in the isotonic bath when the strain is between 2.1% and 11%, qualitatively consistent with the experimental results of Akizuki et al. [27]. Interestingly, according to theory, the trend reverses itself once again at higher clamping strains.

When using a discrete fiber distribution with three orthogonal bundles, the clamping stress increases with increasing NaCl concentration irrespective of the magnitude of the clamping strain (results not shown). Thus, a discrete fiber bundle representation is unable to capture the complex response observed experimentally.

4. Discussion

The concept of modeling continuous fiber angular distributions explicitly can be traced to the work of Lanir [37, 38]. This modeling approach has been used for lung tissue [47], pleura and mycoardium [39], the aortic valve cusp [48], and pericardium [49]. In articular cartilage however, constitutive models that have incorporated fibers explicitly have used a discrete number of bundles, most typically three orthogonal bundles [7, 5052], though the studies of Wilson et al. [8, 19, 29] have used up to four primary fiber bundles and 13 randomly oriented secondary bundles. Other studies used the framework of Conewise Linear Elasticity theory [35] to represent the tension-compression nonlinearity of cartilage without explicit modeling of its fibers [4, 5, 53]; more recently, we showed that this framework is substantially equivalent to modeling three orthogonal fiber bundles as well [12]. Therefore the current study, which adopts the continuous fiber angular distribution model for cartilage, provides a novel opportunity to examine those experimental observations that had yet to be predicted from theory.

Results demonstrate that modeling the solid matrix of articular cartilage with a continuous fiber angular distribution, where the fibers can support tensile loading only, predicts many experimental responses previously reported in the literature. Most of these predictions could not be obtained from constitutive models which use either a few discrete fiber bundles, as illustrated with orthogonal fiber bundles in this study, or phenomenological constitutive relations which do not model fibers explicitly. It is possible that some of the salient features observed with the continuous fiber distribution model adopted here may also be captured, to variable extents, by discrete fiber models that use more than three orthogonal bundles [8, 19, 29], though these investigations have not yet been reported.

As shown in our earlier study [12], a purely fibrillar matrix where fibers can only support tensile loading is inherently unstable, unable to support uniaxial compression or tension without collapsing on itself. Stability of such structures is provided by the ground matrix, represented here by the Donnan and entropic osmotic swelling of the proteoglycans. In this context, the continuous fiber angular distribution, coupled with the osmotic swelling response of proteoglycans, can predict the transition of the incremental Poisson’s ratio from very small values in unconfined compression to very high values in tension. The response in unconfined compression can be explained by the large resistance to lateral expansion imparted by the fibers loaded in tension, relative to the small resistance to axial compression imparted by the osmotic modulus of the proteoglycan ground matrix, thus producing a small Poisson’s ratio. In uniaxial tension, the natural tendency of fibers to collapse in the lateral direction is only resisted by the much softer ground matrix, whereas the axial extension is resisted by the much stiffer fibers themselves, thus producing a very high Poisson’s ratio. However, as more and more fibers realign themselves along the loading direction, fiber recruitment diminishes, thus leading to a reduction in Poisson’s ratio as tensile strains further increase (Figure 2a). This tendency for fibers to collapse along the uniaxial tensile loading direction is not manifested by orthogonal fiber bundles when loaded along their principal axes; consequently, these discrete bundle models are unable to predict a large Poisson’s ratio in tension (Figure 2b).

A particularly intriguing outcome of this study is the ability of the continuous fiber angular distribution model to predict the complex uniaxial tensile response of cartilage in varying bathing solution salt concentrations (Figure 7). The experimental results of Akizuki et al. [27] have been challenging to model using existing frameworks [30, 54]. The most comprehensive effort appears to be the studies of Loret an Simoes [5, 53, 55], which appeal to partitioning the interstitial water of cartilage into intra-fibrillar and extra-fibrillar compartments [56, 57]. The results of the current study suggest that the behavior observed by Akizuki et al. [27] need not be attributed necessarily to elaborate concepts of proteoglycan physical chemistry in compartmentalized domains of the collagen-proteoglycan matrix. A continuous fiber angular distribution is able to predict their experiments by accounting for the evolution of matrix anisotropy with increasing applied strain, as a result of evolving fiber recruitment. It should be noted that a further challenge in the modeling of the experimental results of Akizuki et al. [27] rests on their use of distilled water. In the absence of sufficient ions in solution, the carboxyl and sulfate groups of proteoglycans may not ionize and the effective fixed charge density may be significantly lower than assumed in the model adopted here [58]. Furthermore, the behavior of charged molecular species in distilled water may deviate substantially from the equations of ideal Donnan equilibrium, due to the process of counterion condensation which may further alter the molecular charge density [59].

The aim of this study was to present a constitutive formulation for articular cartilage which captures several salient phenomena observed experimentally. However, no direct attempt was made to curve-fit experimental responses exactly. It is anticipated that exact agreement with experimental measurements may require further refinements in the various constitutive models formulated above. For example, the Donnan osmotic pressure contribution given in Eq.(3) may need to account for non-ideal physico-chemical behavior, by the introduction of experimentally-validated osmotic and activity coefficients. Similarly, the entropic contribution of proteoglycans enmeshed within cartilage may require more elaborate validations, particularly in the likely event that this contribution depends on their composition and molecular weight distribution. The simple power law proposed for the response of fiber bundles in Eq.(14), or the ellipsoidal fiber angular distributions proposed in Eqs.(15)(16) may possibly need to be reformulated as well, to reflect experimental responses more accurately. It is also important to recognize that experimental results reported in the literature may not necessarily be represented accurately using the ideal conditions employed in the theoretical predictions. For example, when examining the ‘equilibrium’ tensile response of articular cartilage, true equilibrium may not necessarily have been achieved experimentally, as cartilage has been known to require upward of 14 hours to reach stress-relaxation equilibrium in tension [21].

Nevertheless, taken together, the results of this study strongly suggest that modeling the solid matrix of articular cartilage with a continuous fiber distribution provides a more accurate representation of the response of this tissue to various loading conditions. This representation is able to predict observed experimental phenomena that are otherwise difficult to predict using other existing frameworks. While continuous fiber distributions impose an additional computational overhead in numerical implementations, this additional burden appears to be justified by the resulting increase in predictive accuracy.


This study was supported by the National Institute of Arthritis and Musculoskeletal and Skin Diseases of the National Institutes of Health, U.S.A. (AR46532).


The sensitivity of the stress-strain and Poisson’s ratio-strain response to changes in the material coefficients ξi and αi (i= 1− 3) of the fibrous matrix is examined via a parametric analysis. The remaining material coefficients are set to representative values of crF=120Eq/m3,ϕrw=0.8, and ϕ= 0.5. An isotropic distribution is assumed for the fibers, such that ξ1= ξ2= ξ3 and α1= α2= α3. By varying ξi from 2MPa to 10MPa by increments of 2MPa, while keeping αi= 2.5, it is found that the slope of the stress-strain response increases more rapidly in tension, exhibiting a stiffer behavior (Figure 8a). Similarly, it is observed that Poisson’s ratio achieves increasingly higher peak values with increasing ξi, in the range of tensile applied strains (Figure 8b). %

Figure 8Figure 8
Parametric variation of the material parameter ξi from 2MPa to 10MPa, by increments of 2MPa. (a) Sensitivity of the stress-strain response. (b) Sensitivity of the Poisson’s ratio-strain response.

By varying αi from 2.0 to 3.0 by increments of 0.25, while keeping ξi= 6MPa, it is observed that increasing values of αi produce a more gradual rise in stress with strain, effectively extending the “toe region” of the stress-strain response (Figure 9a). Similarly, the peak value of Poisson’s ratio becomes depressed, and its onset is delayed to higher strains, with increasing αi(Figure 9b).

Figure 9Figure 9
Parametric variation of the material parameter αi from 2.0 to 3.0, by increments of 0.25. (a) Sensitivity of the stress-strain response. (b) Sensitivity of the Poisson’s ratio-strain response.

These sensitivity analyses indicate that ξi and αi have distinctive influences on the stress-strain and Poisson’s ratio-strain responses. For example, if experimental responses show a relatively longer toe region in the stress-strain curve, but a relatively elevated Poisson’s ratio in tension, this behavior may be captured with a relatively larger αi (to extend the toe region) and relatively larger ξi (to increase Poisson’s ratio).


1A complete enzymatic digestion of the proteoglycans is an option, albeit imperfect.

2The tensor dyadic products [multiply sign in circle], [multiply sign in circle], and [otimesmac] are described by Curnier et al. [35].

3For comparison purposes, the elasticity tensor of a linear isotropic Hookean elastic solid is λI [multiply sign in circle] I + 2μI [otimesmac] I, where λ and μ are the Lame constants. Thus, in the range of small strains, λ [equivalent] −(π+ J[partial differential]π/[partial differential]J) and μ[equivalent]π.

4When modeling anisotropic fiber bundle properties, it is arguably simpler to assume equal fiber fractions in all directions and allow the m material coefficients ξj and αj to vary with orientation.

5The traction-free initial configuration, denoted with 0, generally differs from the traction-free reference configuration, denoted with r, because of the osmotic swelling pressure (p0≠ 0, pr = 0).

6For example, ν12 [equivalent][partial differential]E22/[partial differential]E11


1. Kempson GE, Freeman MA, Swanson SA. Tensile Properties of Articular Cartilage. Nature. 1968;220(5172):1127–8. [PubMed]
2. Maroudas A. Physicochemical Properties of Cartilage in the Light of Ion Exchange Theory. Biophys J. 1968;8(5):575–95. [PubMed]
3. Cohen B, Lai WM, Mow VC. A Transversely Isotropic Biphasic Model for Unconfined Compression of Growth Plate and Chondroepiphysis. J Biomech Eng. 1998;120(4):491–6. [PubMed]
4. 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(6):576–86. [PMC free article] [PubMed]
5. Loret B, Simoes FMF. Mechanical Effects of Ionic Replacements in Articular Cartilage. Part I: The Constitutive Model. Biomech Model Mechanobiol. 2005;4(2–3):63–80. [PubMed]
6. Lanir Y. Biorheology and Fluid Flux in Swelling Tissues. I. Bicomponent Theory for Small Deformations, Including Concentration Effects. Biorheology. 1987;24(2):173–87. [PubMed]
7. Soulhat J, Buschmann MD, Shirazi-Adl A. A Fibril-Network-Reinforced Biphasic Model of Cartilage in Unconfined Compression. J Biomech Eng. 1999;121(3):340–7. [PubMed]
8. 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–204. [PubMed]
9. Li LP, Herzog W. The Role of Viscoelasticity of Collagen Fibers in Articular Cartilage: Theory and Numerical Formulation. Biorheology. 2004;41(3–4):181–94. [PubMed]
10. 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–7. [PubMed]
11. Park S, Krishnan R, Nicoll SB, Ateshian GA. Cartilage Interstitial Fluid Load Support in Unconfined Compression. J Biomech. 2003;36(12):1785–96. [PMC free article] [PubMed]
12. Ateshian GA. Anisotropy of Fibrous Tissues in Relation to the Distribution of Tensed and Buckled Fibers. J Biomech Eng. 2007;129(2):240–9. [PMC free article] [PubMed]
13. Mow VC, Kuei SC, Lai WM, Armstrong CG. Biphasic Creep and Stress Relaxation of Articular Cartilage in Compression: Theory and Experiments. J Biomech Eng. 1980;102(1):73–84. [PubMed]
14. Lee RC, Frank EH, Grodzinsky AJ, Roylance DK. Oscillatory Compressional Behavior of Articular Cartilage and Its Associated Electromechanical Properties. J Biomech Eng. 1981;103(4):280–92. [PubMed]
15. Soltz MA, Ateshian GA. Experimental Verification and Theoretical Prediction of Cartilage Interstitial Fluid Pressurization at an Impermeable Contact Interface in Confined Compression. J Biomech. 1998;31(10):927–34. [PubMed]
16. Soltz MA, Ateshian GA. Interstitial Fluid Pressurization During Confined Compression Cyclical Loading of Articular Cartilage. Ann Biomed Eng. 2000;28(2):150–9. [PubMed]
17. Mow VC, Gibbs MC, Lai WM, Zhu WB, Athanasiou KA. Biphasic Indentation of Articular Cartilage--Ii. A Numerical Algorithm and an Experimental Study. J Biomech. 1989;22(8–9):853–61. [PubMed]
18. Lanir Y. Biorheology and Fluid Flux in Swelling Tissues, Ii. Analysis of Unconfined Compressive Response of Transversely Isotropic Cartilage Disc. Biorheology. 1987;24(2):189–205. [PubMed]
19. 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(3):357–66. [PubMed]
20. Disilvestro MR, Zhu Q, Wong M, Jurvelin JS, Suh JK. Biphasic Poroviscoelastic Simulation of the Unconfined Compression of Articular Cartilage: I--Simultaneous Prediction of Reaction Force and Lateral Displacement. J Biomech Eng. 2001;123(2):191–7. [PubMed]
21. 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–30. [PMC free article] [PubMed]
22. Wong M, Ponticiello M, Kovanen V, Jurvelin JS. Volumetric Changes of Articular Cartilage During Stress Relaxation in Unconfined Compression. J Biomech. 2000;33(9):1049–54. [PubMed]
23. Wang CC, Deng JM, Ateshian GA, Hung CT. An Automated Approach for Direct Measurement of Two-Dimensional Strain Distributions within Articular Cartilage under Unconfined Compression. J Biomech Eng. 2002;124(5):557–67. [PubMed]
24. Chahine NO, Wang CC, Hung CT, Ateshian GA. Anisotropic Strain-Dependent Material Properties of Bovine Articular Cartilage in the Transitional Range from Tension to Compression. J Biomech. 2004;37(8):1251–61. [PMC free article] [PubMed]
25. Elliott DM, Narmoneva DA, Setton LA. Direct Measurement of the Poisson’s Ratio of Human Patella Cartilage in Tension. J Biomech Eng. 2002;124(2):223–8. [PubMed]
26. 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]
27. Akizuki S, Mow VC, Muller F, Pita JC, Howell DS, Manicourt DH. Tensile Properties of Human Knee Joint Cartilage: I. Influence of Ionic Conditions, Weight Bearing, and Fibrillation on the Tensile Modulus. J Orthop Res. 1986;4(4):379–92. [PubMed]
28. Brown TD, Singerman RJ. Experimental Determination of the Linear Biphasic Constitutive Coefficients of Human Fetal Proximal Femoral Chondroepiphysis. J Biomech. 1986;19(8):597–605. [PubMed]
29. Wilson W, Huyghe JM, Van Donkelaar CC. Depth-Dependent Compressive Equilibrium Properties of Articular Cartilage Explained by Its Composition. Biomechanics and Modeling in Mechanobiology. 2007;6(1–2):43–53. [PubMed]
30. Lai WM, Hou JS, Mow VC. A Triphasic Theory for the Swelling and Deformation Behaviors of Articular Cartilage. J Biomech Eng. 1991;113(3):245–58. [PubMed]
31. Ehrlich S, Wolff N, Schneiderman R, Maroudas A, Parker KH, Winlove CP. The Osmotic Pressure of Chondroitin Sulphate Solutions: Experimental Measurements and Theoretical Analysis. Biorheology. 1998;35(6):383–97. [PubMed]
32. Chahine NO, Chen FH, Hung CT, Ateshian GA. Direct Measurement of Osmotic Pressure of Glycosaminoglycan Solutions by Membrane Osmometry at Room Temperature. Biophys J. 2005;89(3):1543–50. [PubMed]
33. Overbeek JT. The Donnan Equilibrium. Prog Biophys Biophys Chem. 1956;6:57–84. [PubMed]
34. Azeloglu EU, Albro MB, Thimmappa VA, Ateshian GA, Costa KD. Heterogeneous Transmural Proteoglycan Distribution Provides a Mechanism for Regulating Residual Stresses in the Aorta. Am J Physiol Heart Circ Physiol. 2008;294(3):H1197–205. [PubMed]
35. Curnier A, Qi-Chang H, Zysset P. Conewise Linear Elastic Materials. J Elasticity. 1994;37(1):1–38.
36. Bonet J, Wood RD. Nonlinear Continuum Mechanics for Finite Element Analysis. Cambridge University Press; Cambridge; New York, NY: 1997.
37. Lanir Y. The Rheological Behavior of the Skin: Experimental Results and a Structural Model. Biorheology. 1979;16(3):191–202. [PubMed]
38. Lanir Y. Constitutive Equations for Fibrous Connective Tissues. J Biomech. 1983;16(1):1–12. [PubMed]
39. Humphrey JD, Yin FC. A New Constitutive Formulation for Characterizing the Mechanical Behavior of Soft Tissues. Biophys J. 1987;52(4):563–570. [PubMed]
40. Narmoneva DA, Wang JY, Setton LA. A Noncontacting Method for Material Property Determination for Articular Cartilage from Osmotic Loading. Biophys J. 2001;81(6):3066–76. [PubMed]
41. Buschmann MD, Grodzinsky AJ. A Molecular Model of Proteoglycan-Associated Electrostatic Forces in Cartilage Mechanics. J Biomech Eng. 1995;117(2):179–92. [PubMed]
42. Armstrong CG, Mow VC. Variations in the Intrinsic Mechanical Properties of Human Articular Cartilage with Age, Degeneration, and Water Content. J Bone Joint Surg Am. 1982;64(1):88–94. [PubMed]
43. 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(2):456–72. [PubMed]
44. Jurvelin JS, Buschmann MD, Hunziker EB. Optical and Mechanical Determination of Poisson’s Ratio of Adult Bovine Humeral Articular Cartilage. J Biomech. 1997;30(3):235–41. [PubMed]
45. Schmidt MB, Mow VC, Chun LE, Eyre DR. Effects of Proteoglycan Extraction on the Tensile Behavior of Articular Cartilage. J Orthop Res. 1990;8(3):353–63. [PubMed]
46. Grodzinsky AJ, Roth V, Myers E, Grossman WD, Mow VC. The Significance of Electromechanical and Osmotic Forces in the Nonequilibrium Swelling Behavior of Articular Cartilage in Tension. J Biomech Eng. 1981;103(4):221–31. [PubMed]
47. Lanir Y. Constitutive Equations for the Lung Tissue. J Biomech Eng. 1983;105(4):374–80. [PubMed]
48. Billiar KL, Sacks MS. Biaxial Mechanical Properties of the Native and Glutaraldehyde-Treated Aortic Valve Cusp: Part Ii--a Structural Constitutive Model. J Biomech Eng. 2000;122(4):327–335. [PubMed]
49. Sacks MS. Incorporation of Experimentally-Derived Fiber Orientation into a Structural Constitutive Model for Planar Collagenous Tissues. J Biomech Eng. 2003;125(2):280–287. [PubMed]
50. Li LP, Buschmann MD, Shirazi-Adl A. A Fibril Reinforced Nonhomogeneous Poroelastic Model for Articular Cartilage: Inhomogeneous Response in Unconfined Compression. J Biomech. 2000;33(12):1533–41. [PubMed]
51. 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–9. [PubMed]
52. Ateshian GA, Ellis BJ, Weiss JA. Equivalence between Short-Time Biphasic and Incompressible Elastic Material Responses. J Biomech Eng. 2007;129(3):405. [PMC free article] [PubMed]
53. Loret B, Simoes FMF. Mechanical Effects of Ionic Replacements in Articular Cartilage. Part Ii: Simulations of Successive Substitutions of Nacl and Cacl(2) Biomech Model Mechanobiol. 2005;4(2–3):81–99. [PubMed]
54. Myers ER, Lai WM, Mow VC. A Continuum Theory and an Experiment for the Ion-Induced Swelling Behavior of Articular Cartilage. J Biomech Eng. 1984;106(2):151–8. [PubMed]
55. Loret B, Simoes FMF. Articular Cartilage with Intra- and Extrafibrillar Waters: A Chemo-Mechanical Model. Mech Mater. 2004;36(5–6):515–541.
56. Wachtel E, Maroudas A. The Effects of Ph and Ionic Strength on Intrafibrillar Hydration in Articular Cartilage. Biochim Biophys Acta. 1998;1381(1):37–48. [PubMed]
57. Huyghe JM. Intra-Extrafibrillar Mixture Formulation of Soft Charged Hydrated Tissues. J Theor Appl Mech. 1999;37:519–536.
58. Frank EH, Grodzinsky AJ, Phillips SL, Grimshaw PE. Physicochemical and Bioelectrical Determinants of Cartilage Material Properties. Springer-Verlag; New York: 1990.
59. Waigh TA. Applied Biophysics: A Molecular Approach for Physical Scientists. J. Wiley \& Sons; Chichester, England; Hoboken, NJ: 2007.