## Abstract

We investigate the influence that matrix structure, size of diffusing molecules and type and intensity of mechanical loading have on the transport of neutral solutes in articular cartilage. Although this type of investigation has been performed in the past, earlier researchers assumed a constant diffusion coefficient. By contrast, our diffusion coefficient depends on the local deformation of the matrix, and thus varies both in space and in time during an experiment.

We derive a three-dimensional formulation of the problem based on mixture theory and utilize the commercial finite-element code ABAQUS to study it numerically. We also make use of the Cohen–Turnbull–Yasuda model to correlate the decrease of the diffusion coefficient with the increase in tortuosity, owing to the presence of the matrix. Under appropriate circumstances, the equations derived here reduce to the classical convection/diffusion equation and the equations of the biphasic cartilage model. Even though we chose axisymmetric sample geometry for the present calculations, the model can easily be applied to irregular three-dimensional samples.

Our results reinforce and refine previously published studies. The neutral solute's rate of diffusion is reduced under static compression, due to the strain dependence of the diffusion coefficient; an increase in static compression leads to a decrease in the rate of transport of solutes of all sizes. Dynamic loading, on the other hand, augments solute transport due to convection, depending on particle size. The transport of small molecular size solute is moderately enhanced, but only within the surface layer; however, the rate of transport of large molecule solute is greatly increased, even in the deep layer of the cartilage.

## 1. Introduction

Articular cartilage makes the smooth relative motion of the load-bearing surfaces of a joint possible; it also aids the joint in absorbing mechanical shock and distributes joint load more evenly through deformation. In the latter aspect, joint lubrication is analogous to, but more complicated than, soft elastohydrodynamic lubrication. The water content of the cartilage is in the range of 65–80% by net weight. The major structural components of tissue's composite solid matrix are type II collagen fibrils (65% dry tissue weight) interspersed among proteoglycan aggregates (approximately 25% dry weight). Since adult articular cartilage is avascular, the only way chondrocytes can receive nutrients and remove waste product is by diffusive and convective transport across the articular surface. Therefore, an understanding of the mechanism of solute transport under different types and intensities of mechanical loading is essential to our appreciation of the biological activities occurring in the cartilage.

Cartilage is inhomogeneous relative to diffusion. Leddy & Guilak (2003) found the diffusion coefficients for all molecules to be site specific, but also dependent on the size of the molecule. In particular, they observed the diffusion coefficients for 3 and 500 kDa dextrans to vary 1.6 and 2.4 times, respectively, depending on position. Solute transport in cartilage under mechanical load has been studied both experimentally and theoretically. The experimental investigations reveal that static loading leads to reduced diffusion of solutes of all sizes (Burstein *et al*. 1993; Quinn *et al*. 2000, 2001, 2002; Nimer *et al*. 2003). In comparison, dynamic loading enhances the rate of transport of particles, but only the transport of large molecules is affected to a significant degree (Urban *et al*. 1982; Palmoski & Brandt 1984; Katz *et al*. 1986; Ohara *et al*. 1990). Garcia *et al*. (1996, 1998) showed that convection enhanced neutral solute transport by a factor that increased with the size of the solute. Sah *et al*. (1989) and Kim *et al*. (1994) suggested that the mechanical load with frequency higher than 0.001 Hz would increase the solute diffusion in 3 mm diameter discs. Lower‐frequency and static compression would inhibit solute transport.

Theoretically, transient transport of solutes in deforming tissue has received little attention in the past. Bowen (1980) investigated the flow of incompressible fluid constituents through a deformable, incompressible solid although the treatment was kept very general. Bowen commented, ‘The precise formulas for **G**_{d} and **Λ**_{dc} (momentum supply terms) are not important to us here’, and, again, ‘where *ξ*_{d} (related to momentum supply) is a function of *s*_{1},…,*s*_{N} (saturations of the *N* components of the mixture).’ Moreover, since the paper does not directly address transport of molecules in the interstitial fluid of the deformable porous material, Bowen deals with a problem that is related to, but not identical to, the subject problem of this paper. All other authors we know of studied the problem that results from keeping the coefficient of diffusion constant and independent of deformation. Mauck *et al*. (2003) examined convective transport under dynamic loading using mixture theory and studied the diffusion phenomena for a range of the governing non-dimensional parameters. Parameter values corresponding to small solute, *M*∼180 Da showed only a very limited response to dynamic loading, but those corresponding to large particles, *M*∼45 kDa, responded vigorously. For the latter particles, the rate of diffusion increased with increasing amplitude of the forcing function. For simplicity of calculation they investigated one-dimensional unconfined compression under the assumption that the volume fractions of matrix and interstitial water, the hydraulic permeability and the solute diffusion coefficient all remain constant during deformation. Based on a finite-element poroelastic model, Ferguson *et al*. (2004) studied the effect of load-induced interstitial flow on the convective mass transport of solutes in intervertebral discs. In their investigation, the diffusion coefficient was kept constant and the permeability varied with strain. Torzilli *et al*. (1987) investigated the transient diffusion of uncharged molecules through unloaded, undeformed cartilage both experimentally and analytically.

The objective of the present study is twofold: to develop a mixture-theory based, three-dimensional formulation of solute transport in cartilage under diverse loading condition, and to investigate how different mechanical loading conditions influence neutral solute transport in the cartilage. It was thought desirable to provide a formulation that can be readily incorporated into commercial finite-element software. The formulation derived in this study combines transport phenomena with deformation in the biological tissue under specified loading, treating this nonlinear problem for the first time. On implementing the model in the finite-element code ABAQUS, we recognize that convection and diffusion of a solute are analogous to changes in temperature in the presence of fluid motion. To simplify the discussion, the present paper focuses on isotropic tissue under infinitesimal deformation but the formulation can describe solute transport in anisotropic, inhomogeneous and hyperelastic soft tissue undergoing large deformation.

## 2. Governing equations

The specimen we wish to model is submerged in a bathing solution of the neutral solute. We envisage the cartilage tissue as consisting of three intrinsically incompressible phases: an inviscid fluid (superscript ‘w’), a hyperelastic solid matrix (superscript ‘s’) and a neutral solute (superscript ‘f’). It is assumed, therefore, that the mechanical behaviour of the tissue can be characterized by the conservation equations of mixture theory (Atkin & Crain 1976; Bowen 1980; Mow *et al*. 1980; Rajagopal & Tao 1995).

### (a) Conservation of linear momentum

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

**T**^{w}, the matrix stress

**T**^{s}and the solute stress

**T**^{f}(2.2)

For the (hyperelastic) solid matrix and the interstitial (inviscid) fluid, the following stress constitutive equations are postulated (parameters defined in table 1)(2.3)and(2.4)

If, on the assumption of negligible solute volume fraction, *ϕ*^{f}≈0, we disregard the solute stress, equations (2.1)–(2.4) yield the momentum equation for the mixture as(2.5)Here, *ρ*^{s} is the density of the solid, * F* is the deformation gradient tensor,

*ψ*is the Helmholtz free energy, is the right Cauchy–Green deformation tensor for the solid,

*ϕ*

^{s}and

*ϕ*

^{w}are the volume fraction of solid matrix and interstitial water, respectively, and

*p*is interstitial fluid pressure.

For infinitesimal deformation, equation (2.5) reduces to(2.6)where * e* is the infinitesimal strain tensor and

*λ*

_{s}and

*μ*

_{s}are Lame's coefficients. Equation (2.6) is the first governing equation of the biphasic theory (Mow

*et al*. 1980).

The momentum equations for the interstitial water and for the neutral solute, respectively, are(2.7)(2.8)Here, *ρ*^{α} represent the density, *μ*^{α} the chemical potential, **v**^{α} the velocity for the *α* component and *f*_{αβ} is the diffusive drag coefficient for momentum exchange between the *α* and the *β* components. Inertia and viscous effects have been neglected in writing equations (2.7) and (2.8).

Constitutive equations for the chemical potentials *μ*^{w}, *μ*^{f} are borrowed from Mauck *et al*. (2003):(2.9)(2.10)where the , *α*=w, f, are the chemical potentials for water and solute, *R* is the universal gas constant, *θ* is the absolute temperature of the mixture, *Φ* is the osmotic coefficient, *c*^{f} is the molar concentration of the solute, *γ*^{f} is the mean activity coefficient, is the true density for the solute and is the apparent solute density. This last relation follows from the definitions(2.11)Here, *V* is the volume of the mixture, and represent the volume, the mass, the mole number and the number of molecules of the *α* component. As we have no further interest in and for *α*=s, w, we set *c*^{f}=*c*, *γ*^{f}=*γ* and *M*_{f}=*M* in the sequel.

### (b) Conservation of mass

As interconversion of mass among the components is not allowed, mass conservation for the *α* component reads(2.12)

Alternatively, when taking into account the intrinsic incompressibility of the components of the mixture, , we have(2.13)

For a saturated mixture and summation of equation (2.13) over *α* yields(2.14)

Assuming that the solute volume fraction is negligible, *ϕ*^{f}≪1, we find the velocity of the solid phase from the approximation to equation (2.14) as(2.15)

Equation (2.8) can be written in the following form if we express the relative velocity between solid matrix and solute as :(2.16)

Combining equations (2.16) and (2.7), we express the relative velocity between the interstitial water and solute as(2.17)

Introducing this expression into the momentum equation for the water phase, equation (2.7), we obtain(2.18)

By eliminating the relative velocity **v**^{w}−**v**^{s} between equations (2.15) and (2.18), we further obtain(2.19)

Substituting the constitutive relations (2.9) and (2.10) into equation (2.19) yields(2.20)

We now use the abbreviation(2.21)and simplify equation (2.20) to(2.22)

The physical meaning of *k*, as defined in equation (2.21), becomes clear on combining equations (2.9), (2.10) and (2.18):(2.23)Employing equation (2.11) in equation (2.23), we recover Darcy's law(2.24)For *γ*=*Φ*=1 and *f*_{fs}=0,(2.25)This is the second governing equation of the biphasic theory (Almeida & Spilker 1998).

Equation (2.25) clearly shows that the quantity *k* defined in equation (2.21) is the permeability of the matrix. Permeability in the cartilage varies with porosity (Lai & Mow 1980), and indirectly with strain, according to(2.26)where *ϕ*^{w} is related to through equations (2.35*a*,*b*) below.

### (c) The convection and diffusion equation

Equation (2.17) can be solved for the velocity of the solute(2.27)

By introducing equation (2.27) into equation (2.12) for the solute, we obtain(2.28)

Substituting the constitutive equations (2.9) and (2.10) into equation (2.28) and using *ρ*^{f}=*ϕ*^{w}*cM*,(2.29)From the equation of mass conservation for water, equation (2.13), we substitute for into equation (2.29) and obtain the equation for convection and diffusion in a deformable solid,(2.30)Defining the diffusion coefficient by(2.31)we can write equation (2.30) as(2.32)Comparing definitions (2.21) and (2.31), we see that when solute–solid interaction is taken into account, the diffusion coefficient is related to the permeability through(2.33)If, however, we neglect *f*_{fs} in equation (2.32), we obtain(2.34)

Owing to the presence of the matrix, the cross-sectional area available to solute diffusion in the soft tissue is reduced in comparison to its value in unobstructed solvent. The model proposed by Mackie & Meares (1955) that correlates the decrease of the diffusion coefficient with the increase of tortuosity is(2.35a)

This model is based on random walk and is more appropriate for membrane-gels at high degrees of swelling, although Nimer *et al*. (2003) showed that equation (2.35*a*) is suitable for the diffusion of small solute in cartilage.

The effect of loading on the diffusion of large solute should be greater than on diffusion of small solute, a distinction that is not portrayed by equation (2.35*a*). Based on an analysis of the free volume theory of diffusion (Cohen & Turnbull 1959), Yasuda *et al*. (1968) propose(2.35b)The ‘size coefficient’ *K* allows varying the diffusion coefficient according to the relative size of solute molecules, as *K* depends on both the free and the characteristic volumes of the sample (Yasuda *et al*. 1968). Yasuda determines *K*=0.37 for NaCl diffusion (Na∼23 Da, Cl∼35 Da).

The two models depicted in equations (2.35*a*,*b*) were compared for transport of simple electrolyte solutions through membranes by Koter (2002), who used *K*=0.7 and 3.0 for NaCl and HCl solutions. We were unable to find *K* values for large molecules, and arrived at the values in table 2 through comparison with Quinn's data. Comparison of the predictions of equations (2.35*a*) and (2.35*b*) with Quinn's data is shown in figure 1*a* and *b*, respectively; based on these plots, we opt for using equation (2.35*b*).

The diffusion coefficient has the value *D*_{a} in unobstructed aqueous solution, and *D*_{0} is the value it assumes in the interstices of the undeformed matrix of porosity (equation (2.35*b*)). Assuming that an analogous relationship exists between *D*_{a} and the coefficient of diffusion *D* of the deformed matrix of porosity *ϕ*^{w}, and taking into account the relationship(2.36)we obtain the coefficient of diffusion in the deformed matrix:(2.37)

For finite deformation, the convection and diffusion equation for the neutral solute in soft tissue assumes the form(2.38)

For infinitesimal strain and(2.39)

If both *ϕ*^{w} and *D* are constant during deformation and *γ*=1, our convection and diffusion equation (2.32) simplifies to the classical convection/diffusion equation(2.40)

In summary, our model for neutral particle diffusion in a cartilage undergoing large deformation is defined by the following system of coupled, nonlinear equations(2.41)(2.42)(2.43)(2.44)

If interaction between the solute and the solid matrix is neglected, the equations are simplified to(2.45)(2.46)(2.47)(2.48)

## 3. Numerical issues

In the computations, we assume unit value for both the osmotic coefficient *Φ* and the solute activity coefficient *γ*, neglect solute–solid interaction, *f*_{fs}=0 (Huyghe & Janssen 1997; Sun *et al*. 1999), and postulate infinitesimal deformation of the cartilage. Under these conditions, convection/diffusion of the solute is characterized by(3.1)(3.2)(3.3)(3.4)This system of equations is incorporated into the commercial finite-element software ABAQUS, following Ferguson *et al*. (2004). Since convection/diffusion of a solute is analogous to convection/diffusion of heat, numerical solution of equation (3.3) is carried out via the thermal analogy procedure in ABAQUS.

The ABAQUS manual specifies the thermal equilibrium equation(3.5)in a fluid flowing with velocity * v*. Here,

*θ*(

*,*

**x***t*) is the temperature at point

*and time*

**x***t*,

*δθ*(

*,*

**x***t*) is an arbitrary variation of the temperature,

*ρ*(

*θ*),

*c*(

*θ*) and

*k*(

*θ*) are the density, specific heat and conductivity, respectively, of the fluid,

*q*is the heat added per unit volume from external sources,

*q*

_{s}is the heat flowing into the volume across the surface on which temperature is not prescribed (

*S*

_{q}) and

*is the outward normal to*

**n***S*

_{q}. The boundary conditions are over part of the surface,

*S*

_{θ}, and heat flux per unit area entering the domain,

*q*

_{s}(

*), across the rest of the surface.*

**x**The equations are solved as an initial boundary-value problem. The loading protocol is discretized into fixed time-steps. At each time-step, the strain and the interstitial fluid velocity within the tissue are calculated using soil consolidation procedure. Then this new fluid velocity, the new porosity and the new diffusion coefficient are introduced into the convection/diffusion equation.

We intend to validate our model for convection/diffusion of a solute in biological soft tissue by comparing its predictions with the experimental data of Quinn *et al*. (2002). Quinn and co-workers reported on a fluorescence-based method for the observation of solute concentration in statically and dynamically compressed cartilage explants during desorption.

Quinn's experiments were performed on 0.65 mm thick and 2.7 mm diameter cartilage discs, harvested from 18-month-old cows. Prior to the experiments, the discs were soaked in a phosphate-buffered saline (PBS) bath containing either 3 or 40 kDa dextrans, conjugated to tetramethylrhodamine (a fluorescent). After equilibration, the samples were immersed in another PBS bath, large in volume and well agitated to assure zero concentration boundary condition at sample radial edges, and compressed axially between impermeable platens, without radial confinement. Laser confocal measurements of fluorescence intensity were made, after allowing for 1 h of radial solute desorption under the prescribed loading. Unfortunately, the paper does not list the values of some of the parameters pertinent to the experiments; these had to be assumed by us, and are listed in table 2.

To estimate the value of the diffusion coefficient in unstressed cartilage, we curve-fitted the analytical solution for concentration distribution at *D*_{0}=const. (Quinn *et al*. 2002),(3.6)to Quinn's zero load data. Here, *J*_{0} and *J*_{1} are Bessel's functions, *R* is the radius of the sample and *q*_{m} is the *m*th zero of *J*_{0}.

Application of equation (3.6) yielded *D*_{0}=40 μm^{2} s^{−1} for the 40 kDa dextran. However, Quinn's paper does not contain comparable data for the 3 kDa dextran and the value of *D*_{0}=31 μm^{2} s^{−1}, as well as all other parameter values for both solute size, had to be estimated from best fitting our own model to Quinn's non-zero load data.

## 4. Results and discussion

### (a) Validation of the model

We validate our model for convection/diffusion of a solute in biological soft tissue by comparing its predictions with the experimental data of Quinn *et al*. (2002). The input data used in these comparisons are listed in table 2.

Tissue properties vary within different cartilage explants. The tissue properties we estimate for Quinn's two sets of experiments fall within ranges that have been either measured or estimated by other researchers: permeability is based on Mow *et al*. (1980), Poisson's ratio and Young's modulus are derived from Jurvelin *et al*. (1997) and the permeability exponent is from Li *et al*. (2003). Through best fit of our theoretical prediction to the experimental results, we obtain slightly higher diffusion coefficient for 40 kDa dextran than for the 3 kDa dextran. This result might be explained by the site-specificity of the diffusion coefficient that Leddy & Guliak (2003) have discussed: ‘The diffusion coefficients for all molecules varied significantly with depth…’. They showed that the 40 kDa dextran exhibited smaller diffusivity than the 3 kDa dextran when averaged over all zones. However, in the middle and deep zones of the cartilage, the 3 kDa solute may have lower diffusivity than the 40 kDa: Quinn's experiments were carried out in the intermediate and deep zones.

Figure 1*a* plots normalized desorption concentration versus radial distance *R*−*r*, obtained with initial value and boundary value , for all *t*. The symbols in figure 1*a* represent experimental data for 40 kDa dextran diffusion under static compression to *ϵ*=0, 15, and 31% strain. At any radial position *r*, desorption concentration increases as strain increases, that is, , where symbolizes the value of the diffusion coefficient at position *r* under *ϵ*% strain. The continuous curves represent our theoretical prediction at corresponding strain: agreement between theory and experiment is clearly unsatisfactory. We employed the Mackie–Meares tortuosity model, equation (2.35*a*), in these computations. The Mackie–Meares tortuosity model is based on ion diffusion through membranes. Even though Nimer *et al*. (2003) obtained good agreement with this model for small solutes such as Na^{+} (∼23 Da) and leucine (∼131.2 Da), the model seems to be unsatisfactory for large molecules as it makes no allowance for difference in size.

To demonstrate the disparity between tortuosity model equations (2.35*a*) and (2.35*b*), we re-plot the data of figure 1*a* in figure 1*b*, but this time employing the Cohen–Turnbull–Yasuda tortuosity model, equation (2.35*b*), which accounts for the size of the solute; agreement between prediction and data is satisfactory in this case. In all subsequent computations, we shall therefore employ the tortuosity model of Cohen–Turnbull–Yasuda, equation (2.35*b*).

Although Quinn's paper contains yet another set of data for the 40 kDa solute, taken at 46% static strain, it is not plotted in figure 1 to preserve clarity of presentation. This 46% strain data appear to contradict the previous finding; the diffusion coefficient at this strain is larger than it is at 31%, that is, , whereas our current model yields monotonic decrease of the diffusion coefficient with strain for all strains.

Quinn's 40 kDa dextran data suggest two competing effects that influence the rate of diffusion; one of these, the decrease in diffusivity due to pore closing during deformation, dominates at small strain and is well portrayed by equation (2.35*b*) of the model.

Equation (2.35*b*), which represents only one of a variety of possibilities to account for the size of the diffusing molecules, predicts a diffusion coefficient that decreases monotonically with increasing static compression. However, in unconfined compression, the cartilage fibres will progressively be stretched and aligned with the radial direction as load is increased and motion of the solute will favour this direction of decreasing tortuosity. For small amounts of solute, change in fibre orientation will not affect diffusion too much but fibre alignment will aid the diffusion of large solute. Therefore, beyond a certain degree of compression, diffusion of large solutes will increase with further increase of load rather than decrease. Quinn's experimental results show that 40 kDa solute desorbs from cartilage under 46% compression faster than it does from cartilage under 31% compression. To model this phenomenon, representations of tortuosity other than equation (2.35*b*) must be used.

Figure 2 compares Quinn's experimental data with our theoretical prediction for 3 kDa dextran under static load. Agreement is, again, satisfactory. A comparison of figures 1 and 2 indicates that in unconfined compression the effect of cartilage deformation increases with particle size.

Quinn *et al*. noted slight lateral bulging of their samples during loading, indicating less than ideal frictionless contact between cartilage and platens. In consequence, the diffusion coefficient of their experiment might have varied in the axial direction and the reported results are average values. By contrast, theoretical calculations assume no friction between platens and cartilage, thus experiment and calculation might have been performed under slightly different conditions.

In addition to their static experiments, Quinn *et al*. have performed measurements under dynamic loading. Their data are compared with our prediction in figure 3, which plots against *R*−*r*. Even though our model confirms existence of an inflexion point on the curve, agreement between theory and experiment under dynamic loading is less satisfactory than under static loading. It is possible that ignoring *f*_{fs} in the finite-element calculations is the source of this discrepancy.1 Mauck *et al*. (2003) have shown that although *f*_{fs} has diminishing effect on diffusion under static loading, its influence is considerable under dynamic conditions.

By comparing figures 2 and 3, dynamic loading can be seen to enhance the rate of desorption; during the upward stroke, fluid with zero solute concentration, *c*=0, is drawn into the matrix from the PBS bath, while fluid containing at least some solute, *c*>0, is expelled into the bath during the downward stroke, leading to net loss of solute for the cycle.

### (b) Theoretical prediction

Having investigated the validity of our model, we are ready to consider its implications. We select a confined compression test, shown in figure 4, for this purpose. The relevant physical dimensions and the material parameters chosen for this study are listed in table 3. Solute transport is investigated for two different solutes, molecular weight 400 Da and 400 kDa dextran. The molecular weights were selected based on the work of Leddy & Guliak (2003), who list molecules that are physiologically relevant to diffusion studies in cartilage. These range in molecular weight from 180 Da (glucose) to 2500 kDa (proteoglycan). Leddy and Guliak perform diffusion experiments with 3, 40, 70 and 500 kDa dextran. The 500 kDa dextran is similar in size to some extracellular matrix molecules (e.g. fibrocentin or oligomeric protein). Our chosen molecular weights, 400 Da and 400 kDa, are within the experimental range of Leddy & Guliak. The diffusion coefficients in table 3 were chosen based on published experimental data on solute diffusion in cartilage, ranging from 60 to 4 μm^{2} s^{−1} for molecular weight 400 Da–500 kDa dextran (Maroudas 1970; Torzilli *et al*. 1987; Roberts *et al*. 1996; Quinn *et al*. 2000; Leddy & Guilak 2003).

#### (i) Static loading

Under static loading, solute transport occurs by diffusion alone. The solute diffuses into the tissue across its upper surface into the surface layer and then penetrates into the deep layer of the tissue. Figure 5*a*,*b* plots normalized concentration distribution against axial position for various magnitudes of strain, with initial condition and boundary condition . As indicated, static compression leads to decreasing diffusion for both small and large solutes. However, the effect of loading is stronger on the large solute than on the small, as demonstrated in figure 6; this figure plots the normalized reduction in concentration suffered when strain is increased from 0 to 20%, that is . When compressed, the tortuosity of the matrix increases, relative to its stress-free state, equally for solutes of all sizes. However, the reduction of cross-section available to the diffusing molecule is proportionately larger for the large molecule than for the small. As compression is increased, the diffusion coefficient decreases accordingly. At the bottom plate, this reduction for the large solute is almost eight-fold of its value for the small solute.

#### (ii) Dynamic loading

Dynamic loading enhances solute transport by cyclic convection. In confined compression, the top of the sample is covered with the solute-saturated PBS bath. During the upward motion of the platen, solute-saturated fluid, , enters the cartilage through the pores of the platen, in accordance with the boundary condition. As some of the solute diffuses from the surface into the deeper layers while the platen is still in its upward motion, the fluid that is squeezed out through the pores of the platen during compression is already solute-depleted, . The consequence is a net gain of solute by the cartilage for the cycle.

The process of solute enrichment by convection is more effective for the large solute; the large solute tends to lag behind when the fluid is being squeezed out during the compression part of the cycle, thus less large than small solute is carried out by convection. This results in convection having a larger effect on the large solute than on the small one.

The Peclet number, , in the deep layers of the tissue is smaller than in the surface layers, going to zero as *Z*→1.0, due to the diminishing velocity. Additionally, the Peclet number for small solute is smaller than for large solute in the same position in the tissue because of equation (2.35*b*). For higher Peclet numbers, the convective term in equation (3.3) contributes more to solute transport compared with lower Peclet numbers. Therefore, the influence of dynamic loading on the transport of large molecules is stronger than it is on small molecules, and it is stronger in the surface layers than it is in the deep layers. This is demonstrated in figures 7 and 8. Here, we show concentration change with time under three different loading conditions: (i) no load, *ϵ*=0%; (ii) static load only, *ϵ*_{stat}=20%, *ϵ*_{dyn}=0%; and (iii) dynamic load only, *ϵ*_{dyn}=20%, *ϵ*_{stat}=0%. For the small solute, the effect of cyclic convection on solute transport is negligible in the bottom layers, figure 7*a*, becoming somewhat more noticeable within the top layers, figure 7*b*. For the 40 kDa solute, the transport-enhancing effect of convection is significant even in the deep layers; figure 8*a* shows conditions in the surface zone and figure 8*b* shows conditions in the deep layer.

## 5. Conclusions

We modelled the convection and diffusion transport of neutral molecules in articular cartilage. The model is based on mixture theory and accounts for the effects of deformation during loading. Under the appropriate conditions, the model reduces to the well-known bi-phasic theory. The model is implemented in the finite-element code ABAQUS.

The results obtained with our model support and, at the same time, extend past findings. In particular, we show the following.

The diffusion rate of neutral solutes is reduced when the sample is compressed, due to the strain dependence of the diffusion coefficient; an increase in static compression leads to a decrease in the rate of transport of solutes of all sizes.

Reduction in the rate of diffusion under static loading is larger for large solute than it is for small solute.

When the surrounding bath is maintained at saturation level, dynamic loading enhances solute transport due to convection.

Dynamic loading enhances the transport of large particles more than small particles.

The beneficial effect of dynamic loading on solute transport is stronger in the surface layer than in the deep layer.

Even though an axisymmetric element was chosen for the present calculations, the model can be applied to irregular, three-dimensional cartilage samples undergoing large deformation.

## Acknowledgments

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. We also wish to acknowledge Dr J. Novotny for making us aware of the problem of diffusion in cartilage.

## Footnotes

↵We are indebted to one of the reviewers for this observation.

- Received July 28, 2004.
- Accepted February 1, 2005.

- © 2005 The Royal Society