|Home | About | Journals | Submit | Contact Us | Français|
A biphasic mixture model is developed which can account for the observed tension-compression nonlinearity of cartilage by employing the continuum-based Conewise Linear Elasticity (CLE) model of Curnier et al. (J Elasticity 37:1–38, 1995) to describe the solid phase of the mixture. In this first investigation, the orthotropic octantwise linear elasticity model was reduced to the more specialized case of cubic symmetry, to reduce the number of elastic constants from twelve to four. Confined and unconfined compression stress-relaxation, and torsional shear testing were performed on each of nine bovine humeral head articular cartilage cylindrical plugs from 6 month old calves. Using the CLE model with cubic symmetry, the aggregate modulus in compression and axial permeability were obtained from confined compression (H−A =0.64±0.22 MPa, kz = 3.62 ± .97 × 10−16 m4/N.s, r2 =0.95±0.03), the tensile modulus, compressive Poisson ratio and radial permeability were obtained from unconfined compression (E+Y = 12.75 ± 1.56 MPa, ν− =0.03±0.01, kr =6.06±2.10×10−16 m4/N.s, r2 =0.99±0.00), and the shear modulus was obtained from torsional shear (µ=0.17±0.06 MPa). The model was also employed to successfully predict the interstitial fluid pressure at the center of the cartilage plug in unconfined compression (r2 =0.98±0.01). The results of this study demonstrate that the integration of the CLE model with the biphasic mixture theory can provide a model of cartilage which can successfully curvefit three distinct testing configurations while producing material parameters consistent with previous reports in the literature.
The primary mechanical function of articular cartilage is to transmit loads across diarthrodial joints while maintaining low friction and wear. The remarkable ability of this tissue to sustain itself over the human lifetime of seven to eight decades has been investigated both experimentally and theoretically. Mechanical tests of articular cartilage under various configurations have established that the material is anisotropic and inhomogeneous in tension [1,2,3,4], as well as in compression [5,6], that it exhibits a viscoelastic response [7,8,9,10,11,12], and that its response differs in tension and compression [13,14,15,16,17,18]. The latter behavior may be referred to as tension-compression nonlinearity or bimodular response. The tissue has also been shown to exhibit a nonlinear response under finite deformation [2,19,20,21]. Over the years, several theories have been proposed to describe one or more of these characteristic behaviors; some of these theories have been quite successful in predicting experimental responses under one set of testing conditions, but no single model has yet been able to describe the full complexity of the material.
It is now generally accepted that porous media models, such as the biphasic mixture theory of Mow et al.  or the electromechanical model of Frank and Grodzinsky  can appropriately explain the time-dependent viscoelastic response of cartilage as resulting primarily from the dissipative drag of interstitial fluid flowing through the porous solid matrix. These models have been employed most frequently under the assumption of a linear isotropic solid matrix undergoing infinitesimal deformation. Good curvefitting and predictive ability have been found between the linear isotropic biphasic theory and experimental results in confined compression creep and stress-relaxation [11,23,24] and dynamic loading [22,25]. The theory has been somewhat less successful in curvefitting the early time response of in situ cartilage indentation tests with a flat porous indenter [26,27], and unsuccessful in predicting the response in unconfined compression between impermeable platens [28,29,30,31]. However, Cohen et al. [14,32] demonstrated that the linear transversely isotropic biphasic theory could provide successful curvefits of indentation, unconfined compression and confined compression tests, suggesting that difficulties encountered with the isotropic model could be overcome by modeling tissue anisotropy. Similarly to the earlier work of Lanir , these authors proposed that the use of a transversely isotropic model was a first approximation to the modeling of tension-compression nonlinearity, with the axial Young's modulus representing the compressive modulus, and the transverse Young's modulus representing the tensile modulus. Alternatively, Suh and Bai  were able to show successful curvefits of confined and unconfined compression and indentation experiments using a linear isotropic poroviscoelastic biphasic model first proposed by Mak  and also used by Setton et al. .
In a recent study, Bursac et al.  suggested that the linear transversely isotropic biphasic model, while successful in curvefitting confined and unconfined compression data, could not properly predict the radial compressive stress measured on the side wall of the testing chamber during confined compression, using the methodology of Khalsa and Eisenberg . At the same time, Soulhat et al.  proposed a fibril-network reinforced model for the solid phase of biphasic cartilage which incorporates tension-compression nonlinearity, where the collagen fibrils can only sustain tensile stresses. They demonstrated that this model could successfully curvefit the response of cartilage in unconfined compression.
In the present study, we propose to similarly formulate a biphasic model which can account for the observed tension-compression nonlinearity of cartilage. This formulation employs the continuum-based Conewise Linear Elasticity (CLE) model of Curnier et al.  to describe the solid phase of a biphasic mixture, rather than the transversely isotropic linear biphasic theory proposed by Cohen et al.  or the composite theory approach proposed by Soulhat et al. . It is the hypothesis of this study that such a model can successfully curvefit and predict experimental measurements in confined and unconfined compression. The long term objective of this study is to develop a single constitutive model of cartilage which can describe its mechanical response in multiaxial, tensile, compressive and shear testing configurations under small strains. The significance of this long term aim is to provide a better understanding of the state of stress within articular cartilage under normal physiologic conditions and the mechanism by which changes in the stresses may be implicated in the etiopathogenesis of osteoarthritis; furthermore, a constitutive model which more accurately reflects the mechanics of cartilage under various loading conditions can provide a better assessment of the role of interstitial fluid pressurization in reducing the friction and wear of articular cartilage .
In the biphasic theory of Mow et al. , articular cartilage is modeled as a binary mixture of an intrinsically incompressible solid phase, representing primarily the collagen fibers, proteoglycans and chondrocytes, and an intrinsically incompressible fluid phase representing the interstitial water (up to 85% by weight in healthy cartilage). The total stress tensor σ inside the tissue is the sum of the interstitial fluid pressure p and the elastic or effective stress σe resulting from deformation of the solid matrix,
where I is the identity tensor. The momentum equation for the mixture, under quasi-static conditions and in the absence of body forces is given by
where · denotes the divergence operator and is the gradient operator. The continuity equation for the mixture is given by
where w = f (vf − vs) is the flux of fluid relative to the solid, f is the fluid volume fraction (tissue porosity) and vs, vf are the solid and fluid phase velocities, respectively.
These governing equations need to be complemented by constitutive relations to provide a complete formulation of the theory. As in the original biphasic theory (Mow et al., 1980), the relation between fluid flux and pressure gradient is taken to be linear, as per Darcy’s law:
where k is the permeability tensor.
Since cartilage has been shown to have different properties in tension and compression, it is proposed to use the Conewise Linear Elasticity theory developed by Curnier et al.  to account for this “bimodular” response. Furthermore, because it has been shown experimentally that the tensile properties of cartilage differ along directions parallel and perpendicular to the split lines [1,2,3,4], it is proposed that a most general infinitesimal strain formulation for cartilage should accommodate orthotropic symmetry, since two preferred directions (or planes of symmetry) which are orthogonal automatically define a third direction (plane) orthogonal to the first two. This is achieved by employing the “orthotropic octantwise linear elasticity” model of Curnier et al.  for the elastic stress of the biphasic mixture, which reduces to:
with summations implicit over a and b (a,b = 1,3; b ≠ a), and λab = λba ; tr(.) is the trace operator which yields the first invariant of its tensorial argument, and . In this expression, E is the infinitesimal strain tensor related to the solid displacement u through E = (u + uT)/2. Aa is a texture tensor corresponding to each of the three preferred material directions defined by the unit vector aa (aa · aa = 1, no sum over a, • denoting the dot product of vectors), with Aa = aa aa ( denoting the dyadic product of vectors); from its definition it can be recognized that Aa is symmetric, i.e., . For orthotropic material symmetry, aa · ab = 0 when b ≠ a, and the three planes of material symmetry are defined by their unit normal vectors aa. The three preferred directions are generally taken to be: a1 parallel to the split line direction, a2 perpendicular to the split line direction, and a3 normal to the articular cartilage surface (Fig. 1). The term Aa : E represents the component of normal strain along the preferred direction aa. Note that the tension-compression nonlinearity is embodied in the fact that
This signifies that the material properties λaa differ whether the normal strain component along the direction aa is compressive or tensile. Therefore, the material properties of this model are λ−aa, λ+aa, λab, and µa (a,b = 1,3; b ≠ a, no sum), or more explicitly, the twelve constants λ−11, λ+11,λ−22,λ+22,λ−33,λ+33, λ23,λ13,λ12,µ1,µ2,µ3. There are eight fourth-order elasticity tensors for this orthotropic octantwise linear elasticity model, since Aa : E could be either negative or positive along any of the three possible directions, a=1, 2 or 3. The convexity of the energy potential is ensured by the positive definiteness of each of these eight elasticity tensors ; the general relations for satisfying these constraints can be found in the standard treatment of orthotropic materials (e.g., ).
In smooth linear elasticity, there are nine material constants in an orthotropic model; thus the CLE theory requires three additional constants. A major advantage of this theory is that it is formulated in invariant form. It is based on the formulation of a piecewise quadratic energy potential which yields the above piecewise linear stress-strain law with piecewise constant elasticity modulus .
A complete set of governing equations for this biphasic mixture model employing the CLE theory for the solid phase is obtained by substituting the constitutive relation of Eq.(4) into the continuity equation, Eq.(3), and the relation of Eq.(5) into the momentum equation, Eq.(2), recognizing that vs = Dsu/ Dt in Eq.(3) (where Ds/ Dt is the material derivative with respect to the solid phase). The primary unknowns of these two differential equations are therefore the solid displacement u and fluid pressure p, which can be obtained for any given problem by solving these equations with the proper boundary conditions. This model is hypothesized to be the most general formulation for predicting the experimental behavior of cartilage under infinitesimal strains with a piecewise linear stress-strain law, which remains to be verified from experimental studies. It can potentially model tension-compression nonlinearity as well as tissue anisotropy. As in all constitutive relations, spatial tissue inhomogeneity can be accounted for by varying the material properties as a function of position, e.g., through the depth of the cartilage layer; since the direction of split lines is also likely to vary from site to site on an articular surface, the orientation of the preferred material directions may also be inhomogeneous.
Experimental validation of the biphasic CLE model with its twelve elastic constants cannot be achieved realistically within a single study. The model described above is likely to require several studies to fully investigate its validity and predictive ability. Typically, it would require testing of cartilage samples in tension, compression and shear with respect to all three planes of symmetry; however, due to specimen size limitations, there remain many challenges to testing articular cartilage in tension along the direction perpendicular to the articular surface (a3), in compression parallel and perpendicular to the split lines (a1, a2), or in torsional shear about a1 and a2. As an intermediate step in this process, it is proposed to raise the level of material symmetry of the octantwise orthotropic CLE model to that of cubic symmetry, which is still defined by the same three texture tensors and constitutive relation of Eq.(5), but with the following identities:
In this cubic symmetry model, the number of elastic material constants reduces from twelve to four. The primary limitation of this model, in relation to the known experimental behavior of articular cartilage, is that it cannot distinguish between the tensile responses of cartilage along directions parallel and perpendicular to the split lines; the tensile properties along those two directions would be averaged into a single value. Its major advantage however is that it reduces the number of material constants considerably, while preserving the tension-compression nonlinearity response, which may be useful as a first step in the process of validating this theoretical framework for cartilage. Therefore, the experimental validation performed in this study is achieved using the testing configurations of confined compression along a3, unconfined compression along a3, and torsional shear about a3, which unlike tensile tests tend to be less sensitive to the averaging of tensile properties along a1 and a2.
For cubic symmetry, it is straightforward to derive constraint equations for the material constants to satisfy positive definiteness of each of the elasticity tensors which arise for the various combinations of Aa : E positive or negative; if it is also required that uniaxial loading produce an axial normal strain of the same sign as the corresponding normal stress, the constraints can be combined into the following expressions:
where H−A = λ−1 + 2µ and H+A = λ+1 + 2µ. It will become apparent below that H−A is the confined compression modulus and that the expressions in Eq.(8 c) are the unconfined compressive (E−Y) and tensile (E+Y) Young's moduli, respectively. For cartilage, the experimental results presented below will demonstrate that H−A ≤ H+A, so that in practice Eq.(8 b)1 will be more constraining than Eq.(8 b)2.
In a bimodular material, the stress-strain relation is strictly non-linear (even if piecewise linear). Therefore, since the principle of linear superposition does not apply, the choice of an appropriate reference configuration for a particular analysis may be of significant importance. It is assumed in this model that the articular cartilage reference configuration of zero-elastic stress and zero-strain is that which prevails under hypertonic conditions, in analogy to the experimental work of Eisenberg and Grodzinsky  and the theoretical framework of the triphasic theory for soft-hydrated tissues  which accounts for swelling (pre-stress) effects. Therefore, under non-hypertonic initial conditions (e.g., in isotonic conditions when cartilage is immersed in 0.15 M NaCl solution) there exists a tensile swelling strain which can be represented by the tensor Es. For the cubic symmetry model, it can be shown that this state of strain is isotropic, i.e., Es = Es I, where Es is the scalar measure of the swelling strain. In theory, every analysis should account for this initial state of strain since choosing between λ+1 and λ−1 is dependent on the sign of the strain relative to the reference configuration, not the strain relative to the initial configuration for a particular problem. In practice however, the magnitude of the swelling strain is relatively small [42,43] in comparison to tissue strains under physiologic loading conditions; this tensile swelling strain may also be offset by the presence of a tare compressive load, which is typically applied under experimental testing conditions. Therefore, it can be shown that in most cases (and for the testing configurations of this study in particular) the choice of tensile versus compressive moduli according to Eq.(6) is only slightly affected whether the initial or reference configuration is used.
All three testing configurations of this study employ the same cylindrical plug of articular cartilage, of thickness h and radius r0. In this study it is assumed for simplicity that the material is homogeneous, with no variations of the properties through the depth of cartilage. Furthermore, all analyses assume that axisymmetric conditions prevail; implicit in this assumption is that the split lines of cartilage are oriented radially relative to the specimen axis of symmetry. These simplifying assumptions are conceptually similar to those of Cohen et al.  and Soulhat et al. . Because cartilage cylindrical plugs are easily harvested with their longitudinal axis along a3, the cylindrical coordinate axes for the analysis coincide with the preferred material directions, i.e., [a1 ] = [1 0 0]T, [a2 ] = [0 1 0]T and [a3 ] = [0 0 1]T. From these vectors the texture tensors can be easily deduced and it is found that A1 : E = Err, A2 : E = Eθθ and A3 : E = Ezz, which are the normal strain components along the radial, circumferential and axial directions, respectively. The permeabilty tensor is also assumed to have its principal directions aligned with these preferred material directions, so that its matrix form in this coordinate system is given by
where kr, kθ, kz are the radial, circumferential and axial permeability constants, respectively.
In the confined compression configuration of this study, the cartilage plug is compressed by a free-draining porous filter inside a cylindrical chamber with impermeable side wall and base (Fig. 2b). It can be shown that for this configuration, the strain tensor, when expressed in matrix form in the current coordinate system, reduces to
where uz is the axial deformation (along a3). Since the only non-zero normal strain is compressive, it follows that A3 : E < 0 and from Eq.(5) and Eq.(7) the matrix of the elastic stress tensor in the current coordinate system reduces to
Substituting these relations into the governing equations yields the following differential equation for the axial displacement:
where uz = uz (z,t). The boundary conditions are formulated for a confined compression chamber with impermeable bottom surface (z = 0) where the displacement and fluid flux are zero, and where the deformation is applied on the opposite surface (z = h) via a free-draining porous indenter. Thus, for stress-relaxation problems,
where ua (t) is the prescribed surface ramp displacement given by
and V0 is the ramp velocity. The general form of Eq.(12) subject to the boundary conditions of Eq.(13)–Eq.(14) is the same as that for an isotropic, linear biphasic model and the solution has been provided previously (e.g., ). Once the solution for uz (z,t) has been obtained, the total reaction force is calculated from
and the interstitial fluid pressure is given by
In the present theory, the confined compression response is governed by the material constants H−A and kz.
In unconfined compression, the cartilage plug is loaded between impermeable platens (Fig. 2c) and it is assumed that the friction between the cartilage and platens is negligible . For this configuration, the matrix of the strain tensor reduces to
where ur is the radial displacement. The radial and circumferential normal strains (Err = ur/r and Eθθ = ur/r) are tensile while the axial strain (Ezz = uz/ z) is compressive; thus, A1 : E > 0, A2 : E > 0 and A3 : E < 0. The matrix of the stress tensor is then given by
For this type of problem, it has been shown  that the axial strain is not a function of the spatial coordinates, so that uz/ z ε(t). The governing differential equation for the radial displacement is given by
where ur = ur (r, t). The boundary conditions for this problem are
with the second boundary condition representing, , while also recognizing that the fluid pressure is ambient at the radial edge, p(r0,t) = 0. In a stress-relaxation experiment, the axial strain is given by
where ua (t) has the same form as in Eq.(14). Once the solution has been obtained, the fluid pressure is given by
and the total reaction force at the platens is
As can be observed from these equations, the unconfined compression response is governed by the material constants H+A, H−A, λ2 and kr. The mathematical solution for this equation is identical to the linear transversely isotropic biphasic analysis provided by Cohen et al. , though the physical interpretation of the material constants is somewhat different.
At equilibrium, when the fluid pressure has subsided, the radial displacement is a linear function of the radial coordinate and thus Err = Eθθ; from this relation as well as Eq.(18), Eq.(19) & Eq.(20), the equilibrium Young's modulus and Poisson's ratio ν−=−Err/ Ezz in compression can be derived,
The expression for E−Y is the same as that of Eq.(8 c)1; re-arranging the above expressions produces E−Y = H−A −2λ2 ν−. When using the constraint of Eq.(8 b)2, then −1 < ν− <1/2, however, for cartilage, the more restrictive constraint on Poisson's ratio in compression derives from Eq.(8 b)1, i.e.,
thus, for positive values of ν−, the greater the ratio of tensile to compressive modulus H+A/ H−A, the smaller the upper bound on Poisson's ratio will be . (Note that conversely, for uniaxial tension, Poisson's ratio would be ν+ = λ2/ (H−A + λ2) and the range of ν+ would effectively be restricted to −1 <ν+ <1/2.)
An additional useful analysis is that of the response of the cylindrical plug under rapid compressive loading, which becomes nearly identical to the response of an incompressible elastic solid [44,45]. Under such isochoric conditions, uz/ z = ε and ur/ r = ur/ r = −ε/ 2 at t = 0+, from which it follows that
Since these stress components are all spatially homogeneous, the interstitial fluid load support at t = 0+ is given by the ratio of the pressure to the total axial stress, i.e.,
Thus, for a linear biphasic material (H+A = H−A), the instantaneous fluid load support is equal to 1/3  while for a bimodular material with H+A > H−A, it can approach unity either as λ2 approaches H−A or when H+A H−A.
In torsional shear about the a3 direction, the specimen is clamped between two porous platens which interdigitate with the tissue to prevent slippage (Fig. 3). The matrix of the strain tensor for this problem reduces to
where uθ is the circumferential displacement. From the CLE constitutive relation,
Substituting these stress components into the momentum equation indicates that there is no fluid pressurization and no transient response for this type of problem; the equation yields the classical torsion solution,
where α is the twist angle per unit cylinder height. From this solution, the non-zero shear stress component is and the torque is obtained by integration,
The angular displacement imparted by the platens to the specimen is θ0 = αh. The torsional shear solution is governed by the material constant µ. The above analysis demonstrates that such measurements can yield the shear modulus from a simple testing configuration; it also demonstrates that torsion about one of the preferred material directions involves neither tensile nor compressive moduli even though the principal normal strains assume both tensile and compressive values.
Three glenohumeral joints of 6 month old bovine were obtained from a local abattoir. After dissection, cylindrical plugs of cartilage (n=9, diam. = 4.78 mm) were punched out from the humeral head and microtomed (Leica Instruments GmbH, Nussloch, Germany, Model SM2400) by 20–100µm at the deep zone to produce a parallel surface to the intact articular side, resulting in a mean±std.dev. thickness of h=0.85±0.14 mm. Specimens were stored at −25°C until ready for use.
Confined and unconfined compression tests were performed on the same apparatus, fitted with different testing chambers and loading platens depending on the type of test (Fig. 2). This apparatus consists of a computer controlled stepper micrometer (Model 18500, Oriel Instruments, Stratford, CT, step size 0.7–1.2 µm) connected in series with a linear variable differential transformer (LVDT) for measuring displacements (HR 100, Schaevitz Sensors, Hampton, VA), and a multiaxial load cell (JR3, 20E12A-M25B, 6-axis load cell, Woodland, CA) mounted on its base. The LVDT provides measurements of the specimen axial deformation (uz) while the load cell measures the axial reaction force (Fc, Fu). The stepper micrometer is driven by a personal computer equipped with a National Instruments data acquisition card (PCI-MIO-16XE, Austin, TX) and running the LabView data acquisition and control software. The output from the LVDT and load cell are also monitored with this card and software. The stepper micrometer was operated under load control, via a feedback control loop using the load transducer, for the application of the specimen tare load; displacement control was employed for the stress-relaxation tests in confined and unconfined compression. The order of testing for confined and unconfined compression was randomized and all tests were conducted in physiological buffered saline (PBS) solution.
For unconfined compression experiments, the tissue was compressed between two impermeable flat platens (Fig. 2c). The fixed, bottom aluminum platen had a 1 mm hole connecting to a piezoresistive microchip pressure transducer (NPC 1210-100G-3N; Lucas NovaSensor, Fremont, CA; range 0–690 kPa gauge pressure), similar to the one used in previous studies [24,25]. The top platen was covered with a 1 mm thick glass slide. The specimen was mounted with its articular surface facing the fluid-filled hole leading to the pressure transducer. A tare load of 0.89 N (50 kPa) was first applied; following equilibrium, as indicated by fluid pressure reducing to zero, a stress relaxation experiment was initiated where the displacement ua (t) of the top stainless steel platen was ramped at constant velocity (V0 =0.25 µm/s) until 10% tissue strain was reached (V0t0 = 0.1h) at which time the displacement was maintained constant. The unconfined compression reaction force Fu(t) and fluid pressure at the center of the articular surface, p(r = 0,t), were recorded as a function of time. At the completion of the stress relaxation test the tissue was unloaded and allowed to re-equilibrate for 1 hour.
For confined compression tests, the specimen was loaded into a confining chamber, 4.78 mm in diameter (Fig. 2b). Loading was applied on the specimen via a porous indenter, 4.76 mm in diameter. The same loading protocol was used as for unconfined compression, starting with the tare load and followed by the axial tissue displacement ua (t). The confined compression reaction force Fc (t) at the articular surface was recorded as a function of time. At the completion of both stress relaxation tests the tissue was frozen at −25°C for shear testing on a subsequent day.
A second apparatus was used for shear tests (Fig. 3), which employed a load cell (Sensotec 0–45 N, Columbus OH) and torque transducer (RTS-10, 0–0.035 Nm, Transducer Techniques, Temecula CA) placed inline and connected to a stepper micrometer (Oriel Instruments, Stratford, CT). For these tests, performed within 15 days of the confined and unconfined compression stress-relaxation tests, the specimen was thawed, allowed to equilibrate at room temperature in PBS for 1 hour, and then mounted between two porous loading platens. As for stress-relaxation tests, the specimen was compressed between the platens to a tare load of 0.89N. Following tare equilibrium, a step rotation of θ0 =0.016 rad was applied about the axis of the cylindrical specimen using a stepper motor (Slo-Syn Model 440, Warner Electric, Ann Arbor MI). The resultant equilibrium reaction torque, T, was measured with the torque transducer; for all specimens, torque equilibrium was achieved within 60 seconds.
For confined and unconfined compression, a least-squares curvefitting algorithm was used to find the best-fit material parameters by matching the experimentally measured transient reaction force with the corresponding theoretical response (Eq.(15) or Eq.(23)). For computational efficiency, the theoretical responses were determined using a numerical finite difference scheme for solving the governing differential equations (Eq.(12) or Eq.(19)); central difference was used for spatial differentiation and backward difference in time, yielding a fully implicit scheme. The spatial solution domain (0 ≤ z ≤ h for confined compression and 0 ≤ r ≤ ro for unconfined compression) was divided into 40 equal intervals and the finite difference solution was obtained at each time step using a linear solver for band-diagonal systems of equations. Numerical integration was employed for calculating the reaction force in unconfined compression, Eq.(23), using the trapezoidal rule. Curvefitting was performed using the BFGS quasi-Newton optimization method for minimizing a function with simple bounds using a finite-difference gradient (IMSL Fortran Numerical Libraries, Visual Numerics, Houston, TX), with the objective function given by the root-mean-square of the residual error between the experimental and theoretical load response.
The parameters H−A and kz were obtained by curvefitting the total load response Fc(t) to the confined compression solution (Fig. 4). Using this value of H−A in Eq.(23), H+A, λ2 and kr were obtained by curvefitting the total load response Fu(t) from unconfined compression (Fig. 5). The shear modulus µ was obtained from Eq.(31) using the equilibrium torque response T, given the specimen radius r0 and thickness h.
Once the material parameters had been determined from curvefitting, the fluid pressure at the articular surface of the specimen was subsequently predicted from theory for unconfined compression using Eq.(22) at r = 0, and compared with the corresponding experimental data to test the predictive ability of the model. Curvefits and predictive comparisons between experimental and theoretical results were assessed with a nonlinear coefficient of determination r2 [24,46].
The mean±std.dev. values of the material properties (n=9) were found to be: H−A=0.64±0.22 MPa and kz=3.62×10−16 ±0.97×10−16 m4/N.s from curvefitting the transient confined compression response (r2=0.95 ±0.03, Fig. 4), H+A=13.2±1.7 MPa, λ2= 0.48±0.23 MPa, and kr= 6.06×10−16 ±2.10×10−16 m4/Ns from curvefitting the transient unconfined compression data (r2 = 0.99±0.002, Fig. 5). A paired two-tailed t-test performed on the radial and axial permeability yielded p=0.001. The model predicted the transient fluid pressure response with an r2 =0.98±0.01 in unconfined compression (Fig. 5). The shear tests produced µ = 0.17 ± 0.06 MPa. Table 1 summarizes the material properties obtained for each specimen by curvefitting the experimental responses, as well as r2 values for the curvefitting of the confined compression load response (c.c.f.) and unconfined compression load response (u.c.f.), and for the prediction of the fluid pressure in unconfined compression (u.c.p.). The elastic constants provided in Table 1 (H−A, H+A, λ2, µ) represent a complete set of properties for a material with cubic symmetry; these properties can be combined to produce alternative material constants as summarized in Table 2 (E−Y, E+Y, λ−1, λ+1, ν−) calculated from the results of Table 1. Using a paired Student t-test, λ−1 (Table 2) and λ2 (Table 1) were found to be statistically different (p<0.001), which technically precludes the possibility that the material is isotropic in compression.
The average material properties and tissue dimensions described above were substituted into the governing equations for unconfined compression to calculate the load supported by interstitial fluid pressure as a function of time during the stress-relaxation test,
where p(r,t) is given in Eq.(22). The magnitudes of Pu (t) and Fu (t) and the fraction of the total load supported by fluid pressure, Pu/ Fu, are plotted as a function of time in Fig. 6 for the average material properties listed in Table 1.
Good curvefits of experimental data were achieved with the CLE-biphasic model described in this study, as may be assessed from the results of Table 1. Importantly, no a priori assumption was made about Poisson's ratio in the process of curvefitting the transient response of articular cartilage in unconfined compression [14,16]. In addition to the successful curvefits, it was demonstrated that the interstitial fluid pressure response in unconfined compression could be predicted from the theory with good accuracy (Fig. 5), as determined from the near-unity correlation between theory and experiment (Table 1, r2 u.c.p.), providing further support toward the validation of this model. This prediction serves as an important further test of the theory since it employs the material properties obtained from other measurements to calculate the fluid pressure.
The long-term hypothesis of this study is that a biphasic mixture model which accounts for the tension-compression nonlinearity and anisotropy of the tissue solid phase can predict the experimental response of articular cartilage under a wide variety of testing configurations. The CLE theory of Curnier et al.  was employed to model a piecewise-linear, or bimodular response of an orthotropic solid phase under infinitesimal strains. Due to the large number of material constants arising in this model, it was proposed to first adopt a simpler model with cubic symmetry to reduce the number of elastic properties from twelve to four. In order to characterize these material properties, along with the tissue permeabilities which govern the transient response, three different testing configurations (confined and unconfined compression and torsional shear) were employed on the same cartilage specimens, each of which yielded a subset of the tissue material properties. Confined compression subjects the material to compressive strains only whereas unconfined compression produces both compressive and tensile strains, thereby providing the ability to estimate both compressive and tensile moduli from these two tests. It was demonstrated that the theory could successfully curvefit the transient stress-relaxation responses of confined and unconfined compression, as assessed by the nonlinear correlation coefficient for each of these tests. The material properties thus obtained demonstrated a very significant disparity between the tensile and compressive moduli of cartilage (H+A H−A or E+Y E−Y), which differ by a factor greater than twenty, confirming the necessity of modeling this bimodular response. The value of H+A, when substituted into the expression for Young's modulus in tension, Eq.(8 c)2, produces a value of E+Y =12.75 ± 1.56 MPa (Table 2). Tensile properties of bovine humeral head cartilage have been reported previously by Woo et al. , under a strain rate of 1%/s, whereas the current study reports the tensile modulus under equilibrium conditions which precludes a direct comparison. For human humeral head cartilage, a recent study has shown that the tensile stress-strain response may be somewhat nonlinear, with the equilibrium tensile modulus averaging 4.2 MPa at 0% strain for specimens harvested parallel and perpendicular to the local split line direction, from the surface and middle zones of cartilage ; however, at 16% tensile strain, this equilibrium modulus increased to 19.6 MPa. In this context, the tensile modulus reported in the current study may properly represent an average modulus for an assumed linear tensile response, supporting the hypothesis that the CLE constitutive model produces material constants that are physically meaningful; nevertheless, in a future development of this theory, it may become necessary to generalize the constitutive relation to accommodate a nonlinear response in tension. Importantly, considering that three independent tests were employed to gather the four elastic constants of the theory, the stability constraints of Eq. (8) were satisfied by the material constants of every tested cartilage sample (Table 1 & Table 2).
Lanir  and Cohen et al. [14,32] were the first to suggest that tension-compression nonlinearity of cartilage could explain the observed experimental response of cartilage in unconfined compression, which until then had remained an elusive goal [28,29,30,31]. They demonstrated this point by simulating the cartilage tension-compression nonlinearity with a transversely isotropic biphasic model using an axial elastic Young's modulus much smaller than the modulus in the transverse plane. This approach was deemed reasonable in unconfined compression, where the normal strains are compressive in the axial direction and tensile in the transverse plane. Bursac et al.  adopted this model in their analysis of cartilage unconfined compression and produced similarly good curvefits of their experimental responses; however, when they tested the same cartilage samples in confined compression, measuring the radial normal stress on the confining chamber’s side wall , they found stress magnitudes two to three times smaller than predicted by the transversely isotropic model. They suggested that their study thus pointed out the limitation of simulating tension-compression nonlinearity with transverse isotropy. Concurrently, Soulhat et al.  proposed to overcome the limitations of this model by using a true bimodular response for the solid phase of the biphasic mixture. This was achieved by employing a composite model for the solid phase where the high-modulus collagen fibrils, which can sustain tensile stresses only, are embedded in an isotropic ground matrix having a lower modulus. These authors were also able to achieve good curvefits of their experimental measurements in unconfined compression, though it remains to demonstrate whether their model could explain the experimental findings of Bursac et al. . The number of elastic material constants in the transversely isotropic model is five, whereas the model of Soulhat et al.  requires three elastic constants and an assumed fiber volume fraction distribution. For the specific configuration of unconfined compression, the governing equations of the models employed by Cohen et al.  and Soulhat et al.  are identical in form to those presented in the current study [Eq.(19)–Eq.(22)], though the constitutive relations and material constants are not the same. As a result, all three models would be expected to equally describe the response of cartilage in unconfined compression; furthermore, a term-by-term comparison of these equations can serve to relate some of the material constants of these three models.
The CLE model adopted in this study offers similar advantages to the model of Soulhat et al. , using a somewhat different approach. Though the bimodular response of cartilage described by the current model is a manifestation of the fibrillar nature of the collagen matrix, there is no explicit attempt to model the microstructure of the material in the CLE model [48,49,50]. This approach is similar to the analysis of composite materials where, for example, the stacking of multiple laminates with different fiber orientations may still achieve a material response with orthotropic symmetry, i.e., with only three preferred material directions. Thus, in the present analysis, the preferred directions aa (a=1,3) do not represent three fibrillar bundle directions per se, but rather the planes of symmetry resulting from fibrillar bundles oriented in multiple directions whose compounded effect is to impart an orthotropic response to cartilage. The choice of orthotropic symmetry is predicated by the presence of split lines in cartilage, as described above in the motivation for the choice of this model.
According to Eq.(11), the radial normal stress on the side wall of the confined compression chamber is given by , therefore a measurement of the equilibrium value of the radial stress under a known axial strain can produce a direct measurement of λ2. Such direct measurements have been performed by Khalsa and Eisenberg  who found that on average H−A =0.40 MPa and λ2 =0.13 MPa (using the current notation for the corresponding material parameters), for 1–2 year old bovine femoral knee cartilage, in the range of 5% to 26% compressive strain. Using the same methodology for radial stress measurement in confined compression, Bursac et al.  reported strain-dependent results of material constants from both confined and unconfined compression of 3–4 week old bovine patellar groove cartilage, with H−A =1.72−0.79 MPa, E−Y =1.50−0.68 MPa, and λ2 =0.40−0.39 MPa for compressive strains ranging from 5.6% to 13.9%; these authors noted that the observed strain-softening response of their study was unusual. In the higher range of compressive strain, their results are relatively close to those of the current study (Table 1 & Table 2; H−A =0.64 MPa, E−Y =0.60 MPa, and λ2 =0.48 MPa), considering differences in the source of material and their removal of ~0.5mm from the superficial zone of their cartilage samples; substitution of their experimentally derived material constants into Eq.(24) above would produce ν− ≈ 0.055 and H+A ≈ 6.7 MPa.
An interesting consequence of using the CLE-biphasic model over the linear isotropic biphasic theory is that the former predicts a greater magnitude of interstitial fluid pressurization than the latter in unconfined compression. Equation (27) and Fig. 6 demonstrate that for the average properties derived from the experiments of this study, the peak magnitude of the fluid load support reaches 97.5% of the total reaction force; this fluid load support remains above 90% for a duration of 210 seconds, confirming that fluid pressurization remains significant for physiological loading durations. In contrast, the linear isotropic biphasic theory predicts a peak fluid load support of only 33% in unconfined compression , even while it predicts a peak of 100% in confined compression  and up to 96% in the contact of cylindrical layers [51,52]. Clearly, as per Eq.(27) and the experimental results in Table 1, the large difference between H+A and H−A, and the closeness of λ2 and H−A both contribute to increase this fluid load support. Effectively, the tension-compression nonlinearity of cartilage promotes greater interstitial fluid pressurization in the tissue than a linear material response would otherwise provide because the higher tensile stiffness causes a greater resistance to radial tissue expansion under axial compressive loading; the implication of this result is that fibrillation of the surface layer of articular cartilage, which is often a pre-cursor to full-blown osteoarthritis, produces a lower tensile stiffness (e.g., ), thus potentially a smaller disparity between the tensile and compressive moduli, that would reduce the cartilage fluid load support according to the present findings. This loss of fluid load support has been previously hypothesized to be a significant mechanism in the pathogenesis and progression of osteoarthritis. Indeed, several recent studies have proposed that cartilage interstitial fluid pressurization is of major importance in relation to mechanical function, such as load support and stress-shielding of the solid matrix [24,25,45,51,52,54] as well as cartilage friction and wear [38,55,56].
A material which exhibits a much greater modulus in tension than compression will necessarily exhibit a small Poisson ratio in compression because the tendency for the tissue to expand radially is strongly resisted by the tensile stiffness. In the CLE model, this is manifested by the restriction of Eq. (25) which, for the average material constants of this study, restricts the compressive Poisson ratio to an upper bound of 0.045. A similar restriction was observed in the fibril-network reinforced model of Soulhat et al. ; these authors expressed concern about this apparent limitation of their study, in light of earlier findings by Jurvelin et al.  who reported optical measurements of Poisson's ratio in compression for bovine articular cartilage averaging 0.186 ± 0.065. In the current study, using Eq.(24)2, we predict from theory that ν− =0.034 ± 0.014 (Table 2). However, in recent preliminary experimental studies , we have performed direct measurements of immature bovine articular cartilage in unconfined compression using chondrocytes as fiducial markers, observed under epifluorescence microscopy as described by Schinagl et al. , and found that Poisson's ratio averaged 0.058±0.023, smaller than previously reported by Jurvelin et al. and in better agreement with the methodology of the present study. It is important to note that because very few studies have reported on direct measurements of the compressive Poisson's ratio of articular cartilage to date, this issue requires further experimental investigations, supported by theoretical predictions which account for the recognized bimodular response of cartilage. In tension, Poisson's ratio has been reported from several studies to exceed unity [2,17,59,60]; whereas the assumption of cubic symmetry in the current study's implementation of the octantwise orthotropic CLE model imposes an upper bound of 0.5 on ν+, the more general orthotropic CLE model can potentially account for such measurements and experiments are currently in progress which analyze the tensile response of cartilage using this theoretical framework.
From the result of this study, the permeability was found to be statistically different in the axial and radial directions, suggesting that the permeability tensor k may not be isotropic; the larger permeability observed along the radial direction may be a reflection of the ultrastructural organization of the collagen fibrils, particularly in the superficial zone where the fibrils are oriented tangentially to the surface, potentially facilitating radial flow while restricting axial flow. It has also been suggested that interdigitation of a porous loading platen with the cartilage may introduce artifacts in the measurement of permeability from confined compression , which may serve as an alternative explanation to the observed differences. However, a more definitive verification of this observation would require direct measurements of the permeability.
A tare load was applied to the cartilage samples which produced an average compressive tare strain of approximately 8%, primarily to ensure complete confinement in the confined compression test. Based on our prior experience [24,25], proper confinement is essential to produce good agreement between theory and experiment. The equilibrium strain of 10% prescribed for the confined and unconfined compression experiments was selected as a compromise between the constraint of using infinitesimal strain theory and the desire to apply more physiological loads. As can be deduced by comparing the general expressions for the Lagrangian and infinitesimal strain tensors, the error introduced when using infinitesimal strain theory for strains of 10% is on the order of 1% strain (or ten percent relative error), which is generally considered acceptable in soft tissue mechanics. If the infinitesimal strain formulation of this model is properly validated from experiments, the natural subsequent step will be to formulate a corresponding finite deformation model to properly account for physiological loading magnitudes.
The results of this study suggest that the CLE model, when integrated within a biphasic mixture, has the potential to model the response of cartilage in a greater number of experimental testing configurations than possible with previous models. This assessment is based on the ability of the model to curvefit the stress-relaxation response of the same cartilage samples in confined and unconfined compression and in shear, as well as its ability to predict the interstitial fluid response in unconfined compression, while accounting for the known tension-compression nonlinearity of the tissue. All four elastic constants of the CLE model with cubic symmetry were obtained from this study, in addition to two out of possibly three distinct permeability constants. These encouraging results motivate further studies whose aim may be to characterize the response of cartilage under other testing configurations, such as confined and unconfined compression and torsional shear along the directions a1 and a2 (parallel and perpendicular to the split lines), as well as indentation tests, and uniaxial or biaxial tensile tests. These series of tests can help determine whether the set of material constants which appear in the more general orthotropic CLE-biphasic model (i.e., the twelve elastic constants and three permeability constants) may reduce to a smaller set through general identities such as, e.g., µ1 ≈ µ2 ≈ µ3 (as appears to be the case for meniscal tissue ) or E−Y1 ≈ E−Y2 ≈ E−Y3 (for bovine cartilage, ). Furthermore, since the tensile properties (e.g., ) and compressive properties (e.g., ) of cartilage are known to vary significantly through the depth of the articular layer, it will be necessary to investigate the depth-variation of the CLE-biphasic material constants as well. Such variations may be important when interpreting the measurements of Poisson’s ratio from unconfined compression  and confined compression , within the context of the CLE-biphasic model.
This study was supported by the National Institutes of Health (1R29-AR43628). We would like to thank Professors P. Zysset, A. Curnier and Q-C He of the École Polytechnique Fédérale de Lausanne, Switzerland, for their helpful comments and for reviewing our derivation of the CLE model with cubic symmetry from the more general formulation of their theory.