## Abstract

The paper presents a three-dimensional microstructural model that was developed to study both equilibrium and time-dependent mechanical behaviour of articular cartilage. The model is based on consolidation theory and makes use of the structure and the mechanical properties of tissue components. Fibrils provide tensile stiffness in this model, whose geometric orientation is described by distribution functions that can be changed freely, making the model applicable to any particular fibril configuration. The model is capable of quantitatively predicting equilibrium properties of the different tests using a single set of model parameters, and explains the anomaly that small compressive modulus and small Poisson's ratio are observed in unconfined compression while tension is governed by large tensile modulus and large Poisson's ratio.

## 1. Introduction

Articular cartilage is a specialized soft tissue that covers the ends of articulating bones, making smooth relative motion of the load-bearing surfaces of the diarthrodial joint possible and, by deforming, distributing joint load more evenly. From a mechanical point of view, cartilage is a multiphase material consisting of a collagen fibril network embedded in proteoglycan skeleton and saturated with interstitial fluids (Minns & Steven 1977; Wilson *et al*. 2004; Li *et al*. 2005). Osteoarthritis, a degenerative disease of the joints, is associated with cartilage deterioration that results from both mechanical and biochemical factors (Setton *et al*. 1999; Nalbant *et al*. 2003; Wilson *et al*. 2004). To understand the onset and subsequent development of osteoarthritis, it is essential to gain full appreciation of both mechanical and biochemical processes as they relate to the functioning of the joint (Mow & Lai 1980; Grodzinsky *et al*. 2000; Quinn *et al*. 2002). A significant component of this understanding is the ability to model cartilage response to various types of mechanical loading, for which we require detailed knowledge of cartilage material properties (Mow *et al*. 1980; Kwan *et al*. 1990; Zhang & Szeri 2005).

The material properties of articular cartilage depend on its composition and microstructure. It is believed that the proteoglycans are mainly responsible for equilibrium stiffness in compression while the collagens govern tensile properties of cartilage (Schinagl *et al*. 1997; Soulhat *et al*. 1999; Elliott *et al*. 2002; Li *et al*. 2005). The transient behaviour is due in part to poroelasticity, i.e. the time dependence of the drag that arises from the relative motion between interstitial fluid and solid matrix (Mow *et al*. 1980; Oloyede & Broom 1991), and in part to the inherent viscoelasticity of the solid matrix (DiSilvestro & Suh 2001; Li *et al*. 2005). The collagen fibrils are inhomogeneously distributed and their orientation varies throughout the tissue giving it an oriented, depth-dependent and layered character (Clark 1991; Mow *et al*. 1992). The deep zone is characterized by fibrils approximately perpendicular to the cartilage–bone interface; the middle zone displays almost random fibril orientation, while the fibrils of the surface zone are arranged in layers parallel to the articular surface (Speer & Dahners 1979; Aspden & Hukins 1981; Mow *et al*. 1992). The concentration of proteoglycans increases with depth from the surface of the cartilage, resulting in depth-dependent increase in the compressive modulus (Schinagl *et al*. 1997; Li *et al*. 2002).

There is evidence that articular cartilage does not behave as classical materials do, its mechanical response is direction dependent and demonstrates significant tension–compression nonlinearity (Huang *et al*. 2005). The cartilage exhibits large transverse deformation when loaded in tension, indicating an apparent Poisson's ratio in tension greater than 0.5 (Woo *et al*. 1979; Elliott *et al*. 2002; Charlebois *et al*. 2004; Huang *et al*. 2005), not unlike that of other soft tissues (Elliott & Setton 2001; Lynch *et al*. 2003). Woo *et al*. (1979) took samples from all three zones and with three different orientations: parallel to the split line (0°), perpendicular to the split line (90°) and inclined at 45° to the split line. Measurements were obtained at 0.01 s^{−1} constant strain rate loading. They found that the direction 0° specimens were the stiffest in both surface and middle zones, but less so in the deep zone. Unfortunately, these were unsteady state experiments and thus not directly comparable to our work. Elliott *et al*. (2002), however, measured lateral Poisson's ratio under equilibrium conditions, on specimens oriented along the split line and taken from surface and middle zones. They found Poisson's ratio to be independent of strain and to be significantly higher in the surface zone (1.87±1.11) than in the middle zone (0.62±0.23). Huang *et al*. (2005) measured tensile modulus and Poisson's ratio , again under equilibrium conditions, on specimens from the surface and middle zones, oriented along and perpendicular to the split line. For humeral head specimens from the surface zone, tended to be greater (1.21±0.29) than from the middle zone (1.1±0.26), while there were no significant differences between specimens from the two orientations relative to split line. For the glenoid, neither zonal (surface: 0.92±0.35, middle: 0.99±0.29) nor directional differences of significance were found in the lateral Poisson's ratio . However, when loaded in compression, the apparent Poisson's ratio of articular cartilage was found to be small (Jurvelin *et al*. 1997; Wang *et al*. 2002). To quote Elliott *et al*. (2002, p. 227), ‘these properties will find use in material models for cartilage anisotropy, with explicit representations of the collagen fibres and their unique contributions to tensile loading which will be important for describing the tensile and compressive behaviours of cartilage and other collagen reinforced soft tissues’. How to quantify the relationship between tensile Poisson's ratio and fibril organization, however, is still a challenging task.

We suppose here that microstructural models are capable of predicting the complex behaviour of articular cartilage once the properties of the constituents and their spatial distribution are incorporated in the model. In particular, the tension–compression nonlinearity may result from the mechanical response of the fibrillar network. As the fibrils are slender and buckle easily when compressed (Fung 1993), they give rise to anisotropic behaviour, even when their distribution is otherwise homogeneous and isotropic. There exist many micromechanical models that can predict the mechanical properties of polymers and fibre composites once the internal structure and mechanical properties of components are known (Hill 1964; Christensen & Waals 1972; Advani & Tucker 1987; Chou 1992). Several microstructural models have been developed in the past for drained elastic properties of soft tissues (Aspden 1986; Farquhar *et al*. 1990; Ault & Hoffman 1992; Schwartz *et al*. 1994). However, none of these models incorporates flow-dependent behaviour of soft tissues. Flow-dependent behaviour and fibril behaviour are included, however, in some other models (Soulhat *et al*. 1999; Wilson *et al*. 2004; Li *et al*. 2005).

Though biphasic mixture theory (Mow *et al*. 1980) and consolidation theory (Biot 1941; Oloyede & Broom 1991) can both explain the flow-dependent behaviour of articular cartilage, a more detailed material description is required to successfully model its mechanical response under multiaxial deformation (Cohen *et al*. 1998; Soulhat *et al*. 1999; Soltz & Ateshian 2000). Soulhat *et al*. (1999) were the first to develop a biphasic model that incorporated orthogonally distributed fibrils. The model was later extended to include nonlinear fibrillar response and viscoelastic behaviour, and was shown to successfully simulate transient behaviour in multiaxial deformation (Li *et al*. 1999; Korhonen *et al*. 2003; Li *et al*. 2005). Orthogonal fibril distribution, however, may be unable to predict tensile Poisson's ratio in excess of 0.5, although anisotropy can do so (LeRoux & Setton 2002; Ting & Chen 2005). Wilson *et al*. (2004) developed a poroviscoelastic fibril-reinforced finite-element model incorporating large primary fibrils based on the arcade model of Benninghoff and random secondary fibrils, where the relative density of these two kinds of fibrils could be adjusted.

There are three objectives to the present study: (i) Develop a three-dimensional microstructural poroelastic cartilage model that recognizes both material anisotropy due to tension–compression nonlinearity as well as structural anisotropy due to inhomogeneous fibril distribution and orientation. (ii) Applying the model, analyse apparent equilibrium properties in tension and unconfined compression for various fibril distribution functions and compare with published experimental data. (iii) Implement the model into commercial finite-element software and apply it to a three-layered cartilage construct, estimating model parameters in curve-fitting to published experimental data. The significant finding is that the tensile Poisson's ratio of structurally isotropic matrix can be greater than 0.5. The fibril distribution function of the model is arbitrary, thus the model can be applied to any specialized fibril configuration to quantify the effect of fibril organization. The orthogonal fibril configuration of Soulhat *et al*. (1999) and the arcade model of Benninghoff (Wilson *et al*. 2004) can be viewed as two particular examples of such distribution functions. For analytical convenience, fibril viscoelasticity and nonlinearity are absent in the present work. Nevertheless, the model has potential application in the mechanical modelling of other soft tissues and fibre-reinforced composite materials.

## 2. Governing equations

We envisage the cartilage tissue as a saturated mixture consisting of a fluid phase (identifier *w*), a porous solid matrix (identifier *s*) and a three-dimensional fibril network (identifier *f*) embedded in the solid. It is assumed, therefore, that the mechanical behaviour of the tissue can be characterized by the conservation equations of both poroelastic theory (Biot 1955; Oloyede & Broom 1991) and mixture theory (Atkin & Craine 1976; Bowen 1976; Mow *et al*. 1980; Rajagopal & Tao 1995).

### (a) Field equations

By neglecting inertia and body force, the momentum equation for the mixture as a whole is(2.1)Here, is the mixture Cauchy stress; it consists of the stress in the fluid , the matrix stress and the fibril stress (2.2)

For the (linearly elastic) solid matrix and the interstitial fluid, the following stress constitutive equations are postulated:(2.3)and(2.4)Here, is the excess stress and is the strain in the solid, and are the volume fraction of solid matrix and interstitial water, respectively, *p* is the interstitial fluid pressure, and are Lame's coefficients and is the identity tensor.

For the fibril network we assume(2.5)where is the fibril volume fraction.

The fibril excess stress will be defined next. To calculate the elastic stress arising from the deformation of the fibril network, the following simplifying assumptions are made: (i) the fibril is straight and can support tensile stress, but not bending or compressive stresses; (ii) bonding between solid matrix and fibril is perfect and constrains perfectly against slippage; (iii) interaction between fibrils is neglected.

With these assumptions, the component of strain along the orientation of a generic fibril is , and the stress generated by this strain in the direction is . Here, the homogenized fibril modulus is related to , the true modulus of the collagen, through . This definition refers to average fibril network stress over total tissue area, in accord with the definition of solid matrix and fluid phase stresses.

Relative to the global reference frame , we characterize the fibril orientation ** n** by the Eulerian angles and , therefore(2.6)and the transformation matrix between the local coordinates and the global reference frame (

**), depicted in figure 1, is(2.7)**

*x*We now define the stress tensor for the single fibril directed along ** n**, written relative to the local coordinate system , by(2.8)

The local form transforms to the global form , i.e. the single fibril stress relative to the global reference frame, through(2.9)

Having thus computed the stress due to a single fibril, the properties of the fibril network are found by weighting this single fibril stress by the probability distribution function of fibril orientation (Chou 1992) and summing for all regions of fibrils in tension. The effective elastic stress tensor of the fibril network is given by(2.10)Here, we employ Gauss–Kronrod integration, the lower and upper limits and are obtained from the inequality constraint of tensile direction at a material element. For small deformation, the general condition is that the strain in direction is greater than zero,(2.11)

The fibril network stress is added to the stress generated by the solid matrix and the stress in the fluid, to obtain the mixture stress tensor(2.12)

Additionally, we assume that Darcy's law governs solid–fluid interaction(2.13)Here, and are the displacement of fluid and solid, respectively, and is the permeability tensor.

### (b) Fibril distributions

Any fibril distribution function must satisfy the constraints(2.14)(2.15)the first of which was stated by Advani & Tucker (1987).

A number of studies has been published on fibril distribution and orientation in articular cartilage. Weiss *et al*. (1968) reported that the middle zone contained randomly oriented collagen fibrils, with little tendency to form bundles. Minns & Steven (1977) employed scanning electron microscopy to find ‘that the superficial zone contains collagen fibrils which are predominantly aligned along one main axis and parallel to the surface’ and ‘fibrils in the deep and middle zones…were generally oriented toward the articular surface’. Clark (1991) observed distinct and well-organized collagen fibres in the radial (deep) zone, while the soft, irregular middle zone ‘resembled cartilage in the initial stages of degeneration’. Fibre organization varied from the centre to the periphery of the tibial plateau cartilage; the deep zone was thick and orderly in the centre, whereas the surface zone was thick and well developed in the periphery. Aspden & Hukins (1981) detected three layers (zones) in patellar cartilage that were characterized by preferred orientation of their collagen fibrils. ‘A surface layer where the preferred orientation is parallel to the articular surface; an intermediate layer in which there is usually no predominant preferred orientation; a deep layer showing a high degree of orientation with the fibrils approximately perpendicular to the articular surface’ (Aspden & Hukins 1981, p. 301).

Our model includes both random and transversely isotropic fibril distributions; we will investigate three kinds of fibril orientations: (i) random, in the middle zone, (ii) axially preferred in the deep zone, and (iii) horizontally preferred in the surface zone. We also compare with the orthogonal fibril network arrangements of Soulhat *et al*. (1999).

For

*random fibril orientation*, the distribution function is constant and its value is computed from(2.16)For

*horizontal fibril alignment*, a Gaussian function is used(2.17)In figure 2, we plot for different values of the variance , which measures the deviation of fibril orientation from the horizontal plane. For , all fibrils are contained in planes parallel to the articulating surface. The distribution function is almost uniform for , meaning that fibril orientation favours all angles, , equally.For

*vertically aligned fibrils*, the distribution function is defined by(2.18)Figure 3 depicts the variation of with , indicating, again, that for we have a uniform distribution while for all fibrils orient normal to the articulating surface.

We note here that both horizontal and vertical distributions are functions of only. The significant difference between vertical and horizontal distributions is that they are centred about different mean values of ( for the horizontal distribution while for the vertical). The variance measures how closely the fibrils are aligned relative to . We may look at these distributions as continuously changing with : when , fibrils are completely in axial directions for vertical fibril distribution, while fibrils are included in transverse planes for the horizontal distribution function.

## 3. Results and discussion

### (a) Equilibrium analysis

The pore pressure vanishes in a state of equilibrium, and the stress–strain relations, equation (2.12), have the form(3.1)when written as a generalized Hooke's law. Here, we used the notation(3.2)

Equations (2.3) and (3.2) use as material characterizers of the solid matrix of the cartilage. It has been customary in cartilage-related work, however, to employ instead the aggregate modulus and Poisson's ratio as the defining elastic parameters. The aggregate modulus is equivalent to the apparent equilibrium modulus in confined compression; this, and Poisson's ratio are related to Lame's constants through(3.3)

In order to compare with experiments from literature, we analyse tension tests on a cubic sample with the loading direction parallel to the cartilage surface (-direction, figure 1), and calculate apparent tensile modulus and Poisson's ratios and , in the lateral and the vertical direction, respectively. In contrast, compression is performed on a cylindrical sample with the loading direction perpendicular to the cartilage surface, resulting in an apparent compressive modulus and one Poisson's ratio. We can compute the apparent equilibrium modulus and Poisson's ratio by solving stress–strain relations to obtain the sequence , by iteration. For example, the iterative method for tension test is as follows.

Impose a small longitudinal strain , and begin by selecting two initial Poisson's ratios , then:

Determine the region of tensile stresses from equation (2.11) and calculate the elasticity tensor from equation (3.2).

Use the stress–strain relations(3.4)to solve for the next value of apparent Poisson's ratios and .

Go to step (1), iterate until convergence is indicated by and

Compute the apparent Young's modulus from equation (3.4)(3.5)

In order to investigate how the apparent equilibrium modulus and Poisson's ratios change with model parameters, we normalize all moduli with the aggregate modulus .The normalized equilibrium modulus and Poisson's ratios depend only on , and fibril orientation. Note that the collagen modulus of our model is equivalent to the individual fibril modulus of the model of Soulhat *et al*. (1999). The homogenized fibril modulus in our model, whereas in Soulhat's model that has two-thirds of the fibrils distributed along planes perpendicular to the load line, the tensile modulus of fibril network . Thus, the fibril modulus of Soulhat *et al*. is related to our fibril modulus by .

For random fibril distribution, figure 4 shows the apparent equilibrium modulus and Poisson's ratio in tension, changing with and (note that equals here). The tensile Poisson's ratio can be greater than 0.5, according to this figure, provided that the ratio is large enough, even though the tissue is structurally isotropic (no spatial preference of fibril orientation). This raises the possibility of tension–compression induced anisotropy even in randomly oriented fibril network. It may also be observed here that the apparent tensile modulus is significantly smaller than the fibril modulus and decreases further with increasing (figure 4*a*). For example, if we use the middle value for matrix Poisson's ratio, then at and 0.127 at . Thus, the tensile modulus is only about 10% of the fibril modulus in the range. In the middle zone where the tensile modulus of articular cartilage is about one order of magnitude larger than the aggregate modulus (Huang *et al*. 2005), the fibril modulus might be two orders of magnitude larger than that. Accepting this to be the case, in the practical range the tensile Poisson's ratio depends mainly on the ratio, though the solid matrix Poisson's ratio may also influence it slightly (figure 4*b*). Huang *et al*. (2005) measured tensile Poisson's ratio in the middle zone at for humeral head and for glenoid. These values are comparable to our model predictions when is of order 100.

Both the apparent equilibrium modulus and Poisson's ratio change with both and in unconfined compression, as indicated in figure 5. The figure also demonstrates that the compressive modulus approaches the aggregate modulus in numerical value when , in agreement with the experimental observation that the equilibrium modulus measured in unconfined compression differs a little from that in confined compression (Cohen *et al*. 1998; Soltz & Ateshian 2000). Under the condition of large , figure 5 finds the compressive Poisson's ratio to be less than 0.1; this value is slightly below the range measured by Jurvelin *et al*. (1997) on uncalcified bovine humeral cartilage, but is consistent with the measurements ( at 5% compression and at 25% compression) of Wang *et al*. (2002). The compressive Poisson's ratio depends strongly on both and , but, if is fixed, it seems linearly dependent on the solid matrix Poisson's ratio (figure 5*b*).

Tension–compression nonlinearity of the tissue can be explained by comparing figures 4 and 5. Under normal conditions, the tensile stiffness is much larger than the aggregate modulus of the solid matrix and the fibril network carries the tensile loads. In unconfined compression, this fibril network resists lateral expansion of the solid matrix and increases equilibrium stiffness (figure 5), suggesting that fibril-reinforcement is also very important in unconfined compression. An incompressible matrix would be unable to support unconfined compressive load without fibril-reinforcement (i.e. ), though it may support load in confined compression.

We now highlight the influence of tissue microstructure on cartilage tensile properties by comparing different fibril orientations under tension. We chose four distribution functions for this purpose: random, orthogonal, horizontal and vertical . The equilibrium tensile stiffness and Poisson's ratio versus matrix Poisson's ratio for these distribution functions are shown in figure 6, where we maintain . In our examples, the direction of tensile loading is parallel to the articular surface, it is obvious therefore that the vertical fibril distribution yields the smallest tensile modulus. However, contrary to expectations, the horizontal fibril distribution does not result in the largest tensile modulus either. Figure 6*a* shows that the orthogonal fibril distribution has the largest tensile modulus at more than twice the value it has for random fibril distribution. To explain this incongruity, we need to look at the load-carrying efficiency (defined here for tension as ) of the different fibril orientations. When a fibril inclines at angle to the load, its contribution to stiffness parallel to the load drops off as (Farquhar *et al*. 1990). For orthogonal fibril distribution, only one-third of the fibrils are in the direction of the load while two-thirds are in compression and do not contribute to tensile stiffness; the efficiency of the orthogonal fibril network is . For random fibril distribution, most fibrils are off axis and, due to transverse contraction, only a small portion of the fibrils is in tension, resulting in load-sharing efficiency. Similarly, the efficiency is about 0.18 for horizontal and 0.04 for vertical fibril distributions. Compared with random fibril distribution, the orthogonal fibril distribution will obviously overestimate the apparent tensile stiffness for the same fibril modulus.

For orthogonal fibril distribution, the apparent tensile Poisson's ratio equals Poisson's ratio of the solid matrix and, therefore, cannot exceed 0.5 (figure 6*b*). Note that only the lateral Poisson's ratio is shown here, but we should bear in mind that the vertical Poisson's ratio is identical to for random and orthogonal distributions, though not for horizontal and vertical distributions. For random distribution, the tensile Poisson's ratio equals about 0.9 and depends only weakly on within the computed range. The lateral Poisson's ratio for horizontal fibril distribution is larger than that for random distribution and increases as , while the reverse tendency is demonstrated for vertical fibril distribution in figure 6*b*. Through comparison of horizontal and random fibril distributions, we may conclude that both the tensile modulus and the lateral Poisson's ratio are larger in the surface zone than in the middle zone, consistent with the experimental measurements of Elliott *et al*. (2002) and Huang *et al*. (2005). For the vertical Poisson's ratio , our calculations indicate no distinct pattern of zonal differences: it can be smaller or larger in the middle zone than in the surface or the deep zone, depending on the matrix Poisson's ratio. Woo *et al*. (1979) observed large thickness reduction in unsteady tensile tests, but his experiments are not directly comparable with our computations for equilibrium vertical Poisson's ratio . More experiments are needed in this area.

To investigate further the effect of fibril organization, figure 7 plots the tensile modulus and Poisson's ratios as functions of the variance for horizontal fibril distribution. The solid matrix Poisson's ratio has the value of in these calculations. Though similar results have also been obtained for vertical fibril distribution, they are not shown here since we have few exploitable tensile data available on specimens taken from the deep zone. In figure 7*a*, the tensile stiffness is seen to vary little for but it decreases gradually for larger, increasing , converging to a constant value as the distribution approaches random distribution. We observe that lateral Poisson's ratio decreases with the increase of (figure 7*b*), while the vertical Poisson's ratio increases (figure 7*c*), but they both approach constant magnitudes as the fibre distribution tends towards random. It is interesting that can be negative when is small and is large (figure 7*c*).

Negative Poisson's ratios are theoretically permissible (Lakes 1987); for an isotropic material, for example, thermodynamics sets the limits on Poisson's ratio at −1.0 and +0.5, and materials with a negative Poisson's ratio, which requires a transverse expansion on stretching, do exist. These include foams (Lakes 1987), biological membranes, for example lipid bilayers (Lakes 2001; Baughman 2003) and ‘much of the crystalline matter in the Universe’ (Baughman 2003).

In summary, fibril organization has a strong influence on both tensile stiffness and Poisson's ratios, thus, the local morphology of fibril network should be considered in the microstructural modelling of the tissue, supporting the conclusion of Wilson *et al*. (2004).

### (b) Time-dependent behaviour

The time-dependent mechanical response of the microstructural model cannot be evaluated analytically, thus we implement it in the commercial finite-element software abaqus, employing the SOIL consolidation procedure. The input parameters are initial void ratio, permeability, two elastic constants for the isotropic solid matrix, fibril modulus and the fibril distribution function. We employ the 4-node axisymmetric porous element (CAX4P) for the solid matrix and a 4-node axisymmetric solid continuum element (CAX4) for the fibril network, which shares its nodes with the porous element. The material behaviour of fibril element is specified by the user-defined subroutine UMAT. In this subroutine, we update the stresses according to the strain field by equation (2.11) at each increment. We use a continuous integral method to calculate the stresses and the stiffness matrix; this approach differs from the approach of Wilson *et al*. (2004). Starting from initial conditions, we calculate the stresses and solve the governing equations step-by-step.

To show the predictive capability of our microstructural model, we consider the articular cartilage as consisting of three zones, with fibril orientation changing from zone to zone. The assumed relative zonal thicknesses is 15% for the surface, 55% for the middle and 30% for the deep zone, in agreement with the findings of Mow *et al*. (1992) and Elliott *et al*. (2002). Fibril distribution is characterized by , , in the deep zone, the middle zone contains randomly distributed fibrils, while the fibrils of the surface zone are distributed according to , . For simplicity, we assume that both the single fibril modulus and the properties of solid matrix are identical in the three zones, though this assumption is not crucial. Note that the fibril distribution functions used here are only a trial geometric configuration that is intended to mimic fibril organization in the three zones but may not be a true representation of any actual structure. The variance is set at , as this value generates a horizontal distribution function that is not unlike the orientation distribution function of Aspden & Hukins (1981) defined in the surface zone.

The boundary conditions of unconfined compression stress–relaxation problem are specified as(3.6)where are components of the displacement vector of the solid phase.

The model parameters were obtained through fitting to the experimental stress–relaxation data of Soltz & Ateshian (2000). Constant permeability was assumed, and the initial void ratio was set at 4.0. We first fit to unconfined compression stress–relaxation data to obtain , , and , then keep these parameters fixed and curve-fit to confined compression stress–relaxation data to obtain . The optimized parameters are

The pore pressure and stress–relaxation curves in figure 8 show good agreement between experimental data and model prediction. The figure contains curve-fit to experimental stress–relaxation data of Soltz & Ateshian (2000) for six-month-old bovine humeral cartilage; figure 8*a* displays axial stress and pore pressure measured from unconfined compression test, along with axial stress from microstructural model curve-fit and corresponding model prediction of pore pressure, figure 8*b* shows axial stress measured from confined compression test, along with the model curve-fit.

To illustrate further the performance of the model, we calculate apparent equilibrium properties under both tension and unconfined compression in the three zones, using the identical set of parameters. Table 1 lists the results of these computations. On examining the entries in this table, we find that the predicted tensile stiffness and lateral Poisson's ratio fall within the experimental ranges of Elliott *et al*. (2002) and Huang *et al*. (2005).

At , the optimized fibril modulus of our model for the Soltz & Ateshian (2000) data is significantly larger than the value reported by Soulhat *et al*. (effective modulus , fibril modulus ), that was obtained when they curve fitted to their own data (Soulhat *et al*. 1999). Li *et al*. (1999), employing a nonlinear version of the Soulhat model, estimated ; this places Li's average value close to that of the Soulhat model, since tensile strains are very small in Li's simulation. If, however, we use Soulhat's model to curve-fit not Soulhat's data but the data of Soltz & Ateshian (2000), we compute the somewhat larger value . The discrepancy, our 125 MPa opposed to 38.2 MPa of the Soulhat model, is due to the difference in assumed microstructure. The orthogonal fibril distribution has higher load-sharing efficiency than the random, the vertical, or even the horizontal distribution in our model. In accordance with the definition and collagen content (Venn & Maroudas 1977), we estimate the true collagen modulus at from our model, and from the model of Soulhat. Elden (1968) estimated the Young modulus of the collagen to be in the range of , while Sasaki & Odajima (1996) measured . Our predicted collagen modulus is close to these values while the Soulhat model appears to yield an underestimate. For this reason, and because the orthogonal fibril distribution cannot predict tensile Poisson's ratio in excess of 0.5, we strongly advocate taking realistic fibril structure into consideration when modelling cartilage; this has also been explored by Wilson *et al*. (2004).

The micromechanical model establishes the connection between bulk behaviour of the tissue and structure and mechanical properties of its components. The bimodular material model of Soltz & Ateshian (2000) is based on macroscopic properties, but our model shows that the bimodular response is a result of fibril behaviour. The elasticity tensor from equation (3.2) is dependent on strain field, thus the macroscopic behaviour depends on loading configurations. The constitutive models that are based on bulk behaviour can only approximate mechanical behaviour of cartilage under some but not all loading conditions. For example, the isotropic model (Mow *et al*. 1980) is a good approximation for confined compression, and the transversely isotropic model (Cohen *et al*. 1998) performs well in unconfined compression. Although DiSilvestro & Suh (2001) have demonstrated that the biphasic poroviscoelastic model could fit unconfined, confined and indentation responses through one set of material properties, it may have difficulty in explaining tension–compression nonlinearity. Our model also shows that the apparent anisotropy can be caused by asymmetric deformation of the tissue when the fibril distribution is isotropic, and that the anisotropy may be enhanced by anisotropic fibril distribution. To account for orthotropic symmetry, the bimodular material model of Soltz & Ateshian (2000) needs 12 elastic constants, yet it still cannot interpret all observed facts (Wang *et al*. 2003). With modern techniques, fibril organization can be determined (Aspden & Hukins 1981) and thus it is possible to construct a realistic fibril distribution function for the tissue. With this fibril distribution, combined with the intrinsic material properties of the tissue components, we will be able to explain the apparent tension–compression nonlinearity and anisotropy of articular cartilage.

The principal limitations of the present study are assumptions of linearity of fibril response, homogeneity of solid matrix and neglect of swelling of the tissue. Wilson *et al*. (2004) and Li *et al*. (2005) have demonstrated that both viscoelasticity and nonlinearity of fibrils play important roles in the mechanical behaviour of the tissue. Wilson *et al*. (2005) have incorporated inhomogeneous structure, nonlinearity and viscoelasticity of fibrils, inhomogeneity and swelling effects of solid matrix into a cartilage model, which is so far almost complete. In the absence of nonlinearity of fibrils, our model will not be able to curve-fit multi-step stress–relaxation in tension or unconfined compression experiments. Therefore, in order to predict all phenomena in cartilage, the limitations of this study should be eliminated in future studies.

## 4. Conclusions

We developed a three-dimensional microstructural poroelastic model to predict the mechanical behaviour of articular cartilage. This model is based on consolidation theory, and is characterized by a fluid phase, a permeable isotropic solid matrix and a three-dimensional fibril network. The solid matrix and fibril network forms a solid composite, whose mechanical properties are determined by the geometric configuration and properties of components. The fibrils are assumed to support tensile load only and are responsible for tension–compression nonlinearity. This model is general and employs arbitrary fibril distribution function, thus it can be applied to any particular fibril configuration.

Three different fibril distribution functions are proposed for the three zones of the articular cartilage. With these fibril configurations, the model has been employed to predict equilibrium properties of the tissue in tension and unconfined compression. Parametric studies have shown that fibril organization strongly influences the tensile stiffness and Poisson's ratios. By incorporating appropriate fibril distribution, the model is capable of quantitatively predicting equilibrium properties of different tests using one set of model parameters, and it explains the anomaly that small modulus and Poisson's ratio are observed in unconfined compression while large modulus and lateral Poisson's ratio govern in tension. The orthogonal fibril distribution, however, cannot explain the tensile Poisson's ratio in excess of 0.5.

Assuming a particular three-layered fibril configuration, the set of model parameters is obtained by fitting to stress–relaxation data in unconfined and confined compression tests. We predict a much higher, and possibly more correct, fibril modulus, than that provided by the model of Soulhat *et al*. (1999). This, again, emphasizes the influence of fibril organization. We, therefore, advocate that realistic fibril structure should be taken into account in mechanical modelling of articular cartilage, supporting the conclusion of Wilson *et al*. (2004). Moreover, we provide a freedom of choice of fibril distribution functions, thus the model has potential application to other fibrous soft tissues and fibre-reinforced composite materials.

## Acknowledgements

During the course of this research, we received financial support from the National Institutes of Health, under the COBRE grant P20RR16458, Dr T. S. Buchanan, PI. This support is gratefully acknowledged.

## Footnotes

- Received March 2, 2006.
- Accepted April 6, 2006.

- © 2006 The Royal Society