PMCCPMCCPMCC

Search tips
Search criteria 

Advanced

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

Modeling of Neutral Solute Transport in a Dynamically Loaded Porous Permeable Gel: Implications for Articular Cartilage Biosynthesis and Tissue Engineering

Abstract

A primary mechanism of solute transport in articular cartilage is believed to occur through passive diffusion across the articular surface, but cyclical loading has been shown experimentally to enhance the transport of large solutes. The objective of this study was to examine the effect of dynamic loading within a theoretical context, and to investigate the circumstances under which convective transport induced by dynamic loading might supplement diffusive transport. The theory of incompressible mixtures was used to model the tissue (gel) as a mixture of a gel solid matrix (extracellular matrix/scaffold), and two fluid phases (interstitial fluid solvent and neutral solute), to solve the problem of solute transport through the lateral surface of a cylindrical sample loaded dynamically in unconfined compression with frictionless impermeable platens in a bathing solution containing an excess of solute. The resulting equations are governed by non-dimensional parameters, the most significant of which are the ratio of the diffusive velocity of the interstitial fluid in the gel to the solute diffusivity in the gel (Rg), the ratio of actual to ideal solute diffusive velocities inside the gel (Rd), the ratio of loading frequency to the characteristic frequency of the gel (f), and the compressive strain amplitude (ε0). Results show that when Rg > 1, Rd < 1, and f > 1, dynamic loading can significantly enhance solute transport into the gel, and that this effect is enhanced as ε0 increases. Based on representative material properties of cartilage and agarose gels, and diffusivities of various solutes in these gels, it is found that the ranges Rg > 1, Rd < 1 correspond to large solutes, whereas f > 1 is in the range of physiological loading frequencies. These theoretical predictions are thus in agreement with the limited experimental data available in the literature. The results of this study apply to any porous hydrated tissue or material, and it is therefore plausible to hypothesize that dynamic loading may serve to enhance solute transport in a variety of physiological processes.

Introduction

Solute transport in biological tissues is a fundamental process of life, providing nutrients to cells and carrying away waste products. In avascular adult articular cartilage, solute transport occurs primarily across the articular surface, with synovial fluid mediating exchanges with the synovium lining the joint capsule [1]. A primary mechanism of solute transport is through diffusion, which has led many investigators to measure the diffusivity of various small and large solutes in articular cartilage (e.g., [211], see Table 1). The mechanism of passive diffusion has been shown experimentally to be enhanced by cyclical loading of cartilage, in the case of a large solute such as serum albumin [12], and by electroosmotic flow [13], both of which mechanisms lead to convective flow within the tissue. However, for smaller solutes, dynamic loading does not enhance transport [12, 14, 15], and convective flow enhances it less significantly than with larger solutes [16]. Urban and co-workers have attributed the disparate outcome under dynamic loading to the magnitude of the ratio of the 'fluid transport coefficient' [17] to the solute diffusion coefficient in the tissue [12, 14], thus defining a non-dimensional governing parameter for this problem. Garcia et al. [13] determined that electroosmotic flow was regulated by the Peclet number, which is a different but related non-dimensional parameter.

Table 1
Molecular weights (MW) and diffusion coefficients for various solutes in aqueous solution (D0) or in a tissue/gel environment (D). Rd is calculated with volumetric water fractions ([var phi]w) obtained from the literature. When not explicitly given, water ...

Studies of articular cartilage metabolism have demonstrated that static loading as well as loading below a characteristic frequency of 0.001 Hz leads to biosynthetic inhibition, whereas dynamic loading stimulates tissue synthesis [1822]. However, whether this enhanced biosynthetic response results from an enhanced nutritional supply remains an unresolved question. Kim et al. have suggested that the stimulation of tissue biosynthetic response under dynamic loading is most likely the result of enhanced fluid flow or changes in cell shape, rather than enhanced nutrient transport [19]. Static compression of articular cartilage has been shown to reduce the diffusivity of various solutes within the tissue, and has been implicated in the altered biosynthetic response of the tissue to static loading [7, 10, 23]. Growth factors, which have been shown to regulate the biosynthetic response of articular cartilage, are generally large solutes with molecular weights on the order of tens of kilodaltons. A further complication of growth factor uptake, their binding to specific protein complexes, was analyzed by Schneiderman et al., who demonstrated that IGF-I binding complexes in normal human articular cartilage are largely excluded from the tissue [8]. Bonassar et al. have shown that dynamic compression accelerates the biosynthetic response of cartilage to free IGF-I and increases the rate of transport of free IGF-I into the matrix, suggesting that cyclic compression may improve the access of soluble growth factors [24]. In a similar fashion, convective diffusion likely aids in the removal of metabolic waste products, creating an environment more suitable for cellular metabolism and matrix biosynthesis [25], as well as influencing the distribution and rate of loss of matrix products from tissue constructs [12, 26].

In studies related to cartilage tissue engineering, it has been suggested that cell growth rates are diffusionally limited [27]. In studies analogous to those carried out on cartilage explants, Buschmann et al. examined the response of chondrocyte-seeded agarose gels and suggested that cell matrix interactions and physicochemical effects may be more important than matrix-independent cell deformation and transport limitations [28]. They also observed that “for dynamic compression, fluid flow, streaming potentials, and cell-matrix interactions appeared to be more significant as stimuli than the small increase in fluid pressure, altered molecular transport, and matrix-independent cell deformation [28].” In our own studies of cartilage tissue engineering, we have found that a dynamic loading regimen applied to chondrocyte-seeded agarose gels over a 28-day period considerably enhances extracellular matrix synthesis as measured by mechanical properties and biochemical composition [29]. As observed in the biosynthetic response of cartilage explants in the study by Bonassar et al. [24], we have also found that the stimulatory effect of dynamic loading is synergistically enhanced by the addition of TGF-β1 or IGF-1 to the culture media [30]. Whether the enhanced extracellular matrix synthesis resulted from the increased nutrient transport alone, cell-matrix interactions alone, or a combination thereof remains unresolved, although the literature findings described above suggest that it is likely that both effects play a role.

The purpose of the current study is to examine the effect of dynamic loading on solute transport within a theoretical context, and to investigate the circumstances under which convective transport induced by dynamic loading might supplement diffusive transport in a porous hydrated tissue or tissue engineered construct. From a theoretical perspective, there is a need to identify the non-dimensional parameters governing this problem, and investigate the response as these parameters are varied over their physiological ranges. The motivation for such a study is the need to estimate, a priori, whether the transport of a particular solute (large or small, such as a growth factor or oxygen) might be enhanced by dynamic loading, and the time constant for the transport of that solute within the tissue. Such theoretical estimates can help determine whether or not the results of a particular experiment might be attributed in part to enhanced nutritional transport. Furthermore, these theoretical results might help to optimize the loading regimen of engineered tissue constructs.

To formulate the governing equations for this problem, we employ the theory of incompressible mixtures [3138] to model the tissue as a mixture of a solid matrix (representing the extracellular matrix/scaffold), a fluid phase representing the solvent (water), and a fluid phase representing the solute. Throughout the text, this mixture model is generically referred to as a 'gel', with the understanding that this gel can represent a biological tissue, an engineered tissue construct, a hydrogel, or any other deformable porous media. The mixture model is useful in this case because it allows us to model transport of a solute in a dynamically loaded gel, and subsumes the classical transport relations, such as Fick's laws, and generalizes them to include the role of gel deformation. For simplicity in this first analysis, the effects of charges in the mixture (e.g., the fixed charge density contributed by proteoglycans and the ionic charges in the interstitial fluid) are neglected. The governing equations of mixture theory are used to solve the problem of solute transport through the lateral surface of a cylindrical gel sample loaded dynamically in unconfined compression, using frictionless impermeable platens. This configuration is employed because it has been frequently used in experimental studies of the biosynthetic response of cartilage explants and chondrocyte-seeded scaffolds [18, 19, 24, 28, 29] as well as in the experimental study of diffusion under cyclical loading [12].

Methods

Only one solute is considered in the analysis of this study, and each of the components of the mixture are considered to be of neutral valence. In mixture theory, given that the mass, volume, and number of moles of phase α are represented by mα, Vα, and nα, respectively, and that V = ΣαVα, the following definitions apply:

φα=dVαdVvolume fractionρα=dmαdVapparent densityρTα=dmαdVαtrue densityc˜α=dnαdVmixture volume - based concentrationcα=dnαdVwsolvent volume - based concentrationMα=dmαdnαmolecular weight,.
(1)

where the superscript ‘w’ indicates the solvent phase. These relations are given in differential form because the variables are defined locally in this continuum model. The molecular weight is constant for a given phase α; furthermore, since each phase is assumed to be intrinsically incompressible, its true density is also constant. From these definitions, it follows that

αφα=1,ρα=φαρTα,c˜α=φwcα=ραMα,
(2)

thus, concentrations and apparent densities can be interchanged using these relations.

Governing Equations

For a cylindrical specimen placed in a bathing solution which contains a higher (or lower) concentration of the solute phase, and loaded between impermeable frictionless platens, the objective of the analysis is to describe the transient response as the solute transports through the tissue, in the presence of a dynamic axial load or strain (Figure 1). The governing equations for this problem are based on the theory of mixtures [3138]; using the notation of Gu et al. [38], the equations are given by

divσI=0,
(3)

ρwgradμ˜w+fws(vsvw)+fwf(vfvw)=0,
(4)

ρfgradμ˜f+ffs(vsvf)+fwf(vwvf)=0,
(5)

div(φsvs+φwvw+φfvf)=0.
(6)

ραt+div(ραvα)=0,α=w,f.
(7)

Figure 1
Schematic of dynamic unconfined compression of gel construct between frictionless impermeable platens in a bathing solution containing an excess of solute.

The first equation represents the balance of linear momentum for the mixture, neglecting inertial effects (which are significant only in wave propagation problems); the second and third equations represent the balance of linear momentum for the solvent phase (α = w) and the solute phase (α = f), respectively (also neglecting inertial effects, along with fluid viscosity). The fourth equation is the balance of mass for the mixture, taking into account that each phase is intrinsically incompressible; and the fifth equation represents the balance of mass equations for each of the solvent and solute phases. In these expressions vα represents the velocity of phase α and fαβ represents the diffusive drag coefficient for the momentum exchanged between phases α and β as a result of molecular collisions (fαβ = fβα). The constitutive relations for the mixture stress σI and the solvent and solute electrochemical potentials, [mu]w, [mu]f, are given by

σI=pI+λs(trE)I+2μsE,
(8)

μ˜w=μ0w(θ)+1ρTw(pRθΦcf),
(9)

μ˜f=μ0f(θ)+RθMflnaf,
(10)

where we assume that the solute volume fraction is negligible, [var phi]f [double less-than sign] 1. This assumption is justified for most circumstances as illustrated by the two following examples: For a small solute such as glucose (Mf=180g/mole,ρTf=1,500g/L) with a typical concentration in media of 1g/L (cf = 1/180 mole/L), the volume fraction according to Eq.(1) is given by φf=Mfc˜f/ρTf=0.0007, and thus negligible compared to unity. With larger solutes such as growth factors, a representative concentration in supplemented media is on the order of 10µg/L for TGFβ-1 (Mf = 25,000 g/mole) or 300µg/L for IGF-1 (Mf = 7,600 g/mole). If the true density of these growth factors is on the order of 1,000g/L, the volume fraction is [var phi]f ~ 10−8 − 3 × 10−7, again negligible compared to unity.

The gel solid matrix is assumed to behave as a linearly elastic isotropic material, with Lamé constants λss, infinitesimal strain tensor E = (gradu + gradT u)/2, and solid displacement u, similarly to the biphasic analysis of unconfined compression by Armstrong et al. [39]. The electrochemical potentials for the solvent and solute have the same form as in classical physical chemistry [40, 41]. In these expressions, R is the universal gas constant, θ is the absolute temperature, p is the fluid pressure (inclusive of osmotic effects), μ0α(α=w,f) are the chemical potentials in a reference configuration (standard state), and af is the activity of the solute. The model assumes isothermal conditions, which is appropriate for physiological systems. The reference configuration is achieved when af = 1; since the solute is dissolved in the same type of solvent inside and outside the tissue, it is necessary for the reference chemical potential of the solute to be the same in both compartments. When defining the solute reference state, the solute activity af must be formulated to account for the possibility that not all of the solvent inside the tissue is accessible to the solute, due to volume exclusion effect [42, 43]. Let Vfw be the solvent volume accessible to the solute, then κf=dVfw/dVw represents the volumetric fraction of accessible solvent, or the solubility of the solute in the gel (0 ≤ κf ≤ 1), and c¯f=dnf/dVfw=cf/κf is the concentration of solute per volume of accessible solvent. The solute activity is then defined as af = γfcf/c0 = γfcffc0 where γf is the solute activity coefficient, and the solute standard state [40] is then defined in the usual manner, e.g., with c0 = 1 M. For simplicity in this study, it will be assumed that the mixture is ideal, i.e., that the osmotic coefficient Φ and solute activity coefficient γf are equal to unity. Substituting the constitutive relations into the balance of linear momentum equations, Eqs. (3)(5) yields

gradp+(λs+μs)grad(divu)+μs2u=0,
(11)

φwgradp+φwRθgradcf+fws(vsvw)+fwf(vfvw)=0,
(12)

φwRθgradcf+fsf(vsvf)+fwf(vwvf)=0.
(13)

The unconfined compression problem has been shown to be one-dimensional when the loading platens are assumed frictionless, because the axial normal strain is uniform [39]. Thus, it can be shown that dependent variables in this problem exhibit the dependencies p = p(r,t), ur = ur (r,t), [partial differential]uz/[partial differential]z = ε(t), vw = vw(r,t), vf = vf(r,t), cf = cf(r,t), where r and z are the radial and axial coordinates, respectively, and t is time, and that all other dependent variables are equal to zero. The axial strain ε(t) may be prescribed (displacement-control analysis) or may be unknown, to be determined subsequently from the analysis given the axial applied load (load-control analysis). Under these assumptions, the governing equations reduce to

pr+HAr(1rr(rur))=0,
(14)

φwpr+φwRθcfr+fws(vrsvrw)+fwf(vrfvrw)=0,
(15)

φwRθcfr+fsf(vrsvrf)+fwf(vrwvrf)=0,
(16)

1rr[r(φsvrs+φwvrw)]+ε˙(t)=0,
(17)

(φwcf)t+1rr(rφwcfvrf)=0,
(18)

where HA = λs + 2μs. The unknowns in these equations are ur,p,cf,vrw,vrf and the objective is to solve these partial differential equations simultaneously for all five unknowns. Integrating Eq.(17) with respect to r and recognizing that vrs=vrw=vrf=0 at r = 0 yields

φsvrs+φwvrw+r2ε˙(t)=0.
(19)

Now Eqs.(15), (16) and (19) can be solved for p/r,vrw,vrf:

pr=[1k+(1DDwf)RθcfφwDwf](vrs+r2ε˙(t))+(1DDwf)Rθcfr,
(20)

vrw=1φw(φsvrs+r2ε˙(t)),
(21)

vrf=Dcfcfr+(1DφwDwf)vrsDφwDwfε˙(t)r2,
(22)

where the diffusion coefficients Dwf, Dsf and permeability k are defined by

Dwf=φwRθcffwf,Dsf=φwRθcffsf,andk=φw2fws,
(23)

and

1D=1Dsf+1Dwf.
(24)

The definition of Dwf, the diffusivity of the solute in the solvent, is consistent with the classical treatment of diffusion [44] as will be shown below. The definition of permeability k is the same as that presented by Lai et al. [45] and is consistent with Darcy's law for permeation of the solvent through a porous gel solid matrix. Dsf is defined in analogy to Dwf; it is inversely related to the diffusive drag coefficient between the solute and gel solid matrix. D is the gel diffusion coefficient of the mixture, representing the diffusivity of the solute against both the solvent and gel solid matrix of the mixture. The result of Eq.(22) can be used in the balance of mass equation for the solute, Eq.(18),

cft1rr{r[Dcfr(1DφwDwf)cfurt+DφwDwfr2cfε˙(t)]}=0,
(25)

where we have substituted vrs=ur/t and assumed that [var phi]w is constant, to produce a generalized Fick's equation which only depends on solute concentration and solid displacement. Equation (22) can also be rearranged to provide the molar flux (per unit solvent area) of the solute relative to the solid phase,

cf(vrfvrs)=DcfrDφwDwfcf(vrs+r2ε˙(t)).
(26)

To get the relative molar flux per unit mixture area, c˜f(vrfvrs), multiply this equation by [var phi]w. Using Eq.(20) in Eq.(14) similarly produces a differential equation which only depends on solute concentration and solid displacement,

r(1rr(rur))[1HAk+(1DDwf)RθHAcfφwDwf](urt+r2ε˙(t))(1DDwf)RθHAcfr=0.
(27)

Thus, to obtain a final solution for this problem, the coupled system of differential equations in Eqs.(25) and (27) for ur(r,t) and cf(r,t) need to be solved, subject to appropriate boundary conditions. The remaining dependent variables can be obtained by substitution of the solution first into Eq.(20) to get the pressure gradient, which can be integrated using a suitable boundary condition as shown below, then into the velocity expressions in Eqs. (21) and (22). It can be verified that in the absence of a solute, cf = 0, Eq.(27) reduces to the governing equation for unconfined compression of a biphasic tissue [39].

Boundary Conditions

The boundary conditions for this problem need to be determined at r = 0 and r = r0. Because of the axisymmetric conditions at r = 0, all velocities are zero at that location, and according to Eqs. (20)(22), it can be noted that

ur(0,t)=0,pr|r=0=0,andcfr|r=0=0.
(28)

At r = r0, jump interface (or boundary) conditions must be satisfied to enforce continuity of mass, momentum and energy at the interface with the external solution [4648], [[σrr]] = 0, [[[mu]w]] = 0, and [[[mu]f]] = 0, which lead to the relations

HAurr|r=r0+λs[ur(r0,t)r0+ε(t)]=0,
(29)

cf(r0,t)=κfcf*,
(30)

p=p*Rθ(1κf)cf*,
(31)

where p* is the ambient pressure and cf* is the external bath solute concentration. Since the mixture is assumed to be ideal, the solute activity coefficients are the same inside and outside the tissue, and thus the solubility is equal to the partition factor [2, 49]. In this study, it is assumed that the external bath is well mixed and that cf* remains constant. The initial conditions can be derived in a similar manner, assuming that the tissue is at equilibrium with its initial environment, consisting of an external solute concentration of c0f and ambient pressure p0:

ur(r,0)=0,cf(r,0)=κfc0f,andp(r,0)=p0.
(32)

If c0f<cf*, we would expect a net transport of solute into the tissue over time, whereas the converse would occur when c0f>cf*. Finally, the axial normal traction component and the total axial load also need to be evaluated,

σzz(r,t)=p(r,t)+λs(urr+urr)+HAε(t),σzz*(r,t)=p*(r,t)andW(t)=2π0r0r[σzz(r,t)σzz*]dr.
(33)

Fick's Laws of Diffusion

Classically, diffusion problems have been analyzed using Fick's first and second laws of diffusion. The simplest case of diffusion is represented by a two-phase mixture, consisting of single neutral solute diffusing in a solvent phase. Under these circumstances, since there is no gel solid matrix ([var phi]s = 0, [var phi]w ≈ 1), D reduces to Dwf and vrs and [epsilon with dot above](t) drop out of Eq.(26) to yield Fick's first law of diffusion,

cfvrf=Dwfcfr,
(34)

whereas Eq.(25) yields Fick's second law of diffusion (in cylindrical coordinates),

cft1rr(rDwfcfr)=0.
(35)

This result confirms that mixture theory is consistent with classical laws.

If we next consider diffusion of a solute in a mixture with rigid porous solid matrix ([var phi]s ≠ 0), then ε(t) = 0 since no solid matrix deformation can take place. In general, vrs need not be equal to zero if there is a rigid body motion, however Fick's laws are formulated for problems with no net convective motion, i.e., with vrs=0, in which case Eq.(26) reduces to

cfvrf=Dcfr,
(36)

while Eq.(25) yields

cft1rr(rDcfr)=0.
(37)

In this case, Fick's laws are also recovered, with a diffusivity D instead of Dwf to account for the presence of the porous solid matrix. From Eq.(24), and given that Dwf,Dsf are positive, it is evident that D < Dwf.

Finally, if the gel solid matrix is deformable, but in the absence of dynamic loading ([epsilon with dot above](t) = 0), the governing equations are found to differ in general from Fick's classical laws, with Eqs. (25)(27) reducing to

cf(vrfvrs)=DcfrDφwDwfcfvrs,
(38)

cft1rr(r[Dcfr(1DφwDwf)cfvrs])=0,
(39)

r(1rr(rur))1HAkurt+RθHADsfcf(vrfvrs)=0.
(40)

It is noteworthy that for the special case when D = [var phi]w Dwf, Eqs. (38)(39) become uncoupled from the solid displacement and velocity, and we recover Fick's first and second laws as in Eqs. (36)(37). For this special category of problems (most applicable to small solutes, as discussed below), the classical closed-form solution of the diffusion equation can be obtained for cf(r,t) and substituted into Eqs. (36) and (40) to obtain a solution for ur(r,t).

Dimensionless Parameters

In order to determine the dimensionless parameters governing the general solution of this problem, it is now necessary to non-dimensionalize the governing equations. If the following non-dimensional variables are proposed,

r^=rr0,t^=Dr02t,u^r=urr0,c^f=RθDHADsfcf,p^=pHA,v^rα=r0Dvrα,W^=WHAr02,
(41)

the governing equations in Eqs.(25) and (27) reduce to

c^ft^1r^r^(r^[c^fr^(1Rd)c^fu^rt^+Rdr^2c^fdεdt^])=0
(42)

r^(1r^r^(r^u^r))(1Rg+Rdc^f)(u^rt^+r^2dεdt^)c^fr^=0,
(43)

where the non-dimensional parameters Rg and Rd are defined as

Rg=HAkDandRd=DφwDwf.
(44)

The non-dimensional parameter Rg = HAk/D represents the ratio of the characteristic diffusive velocity of the solvent in the gel (HAk/r0) to the characteristic diffusive velocity of the solute (D/r0) in the gel. The non-dimensional ratio Rd = D/[var phi]wDwf represents the ratio of the characteristic gel diffusive velocity of the solute (D/r0) to its ideal characteristic gel diffusive velocity ([var phi]wDwf/r0). The ratio D/Dwf is commonly reported in experimental studies of diffusivity (e.g., [43, 50], see Table 1).

The pressure gradient, Eq.(20), and molar flux of solute relative to the solvent, Eq.(26), similarly reduce to

p^r^=(1Rg+Rdc^f)(u^rt^+r^2dεdt^)+c^fr^.
(45)

c^f(v^rfv^rs)=c^fr^Rdc^f(u^rt^+r^2dεdt^).
(46)

The boundary and initial conditions and the axial normal traction and load, when using non-dimensional variables, remain essentially unchanged,

u^r(0,t^)=0,c^fr^|r^=0=0,andp^r^|r^=0=0.
(47)

u^rr^|r^=1+λsHA[u^r(1,t^)+ε(t^)]=0,c^f(1,t^)=c^f*,andp^(1,t^)=p^*,
(48)

u^r(r^,0)=0,c^f(r^,0)=c^0f,andp^(r^,0)=p^0.
(49)

σ^zz(r^,t^)=p^(r^,t^)+λsHA(u^rr^+u^rr^)+ε(t),andW^(t^)=2π01r^σ^zz(r^,t^)dr^.
(50)

The parameters, Rg and Rd, along with the ratio λs/HA = νs/(1 − νs) (where νs is the equilibrium Poisson's ratio of the gel solid matrix), are three of the non-dimensional parameters which govern this solution. The non-dimensionalization of variables proposed in Eq.(41) are valid as long as HA, D and Dsf are non-zero and finite. If special cases need to be examined where one or more of these conditions does not hold, it is simpler to use the dimensional form of the governing equations, or use an alternative non-dimensionalization scheme.

It is evident that the governing equations for this problem, Eqs. (42)(43), are nonlinear so that a closed-form analytical solution is not available in the most general case. A numerical scheme must be employed to achieve a general solution.

Solution

For the general problem considered in this study, a numerical solution can be obtained for the set of nonlinear partial differential equations described above using a finite difference scheme. In this section, results are presented for various choices of the parameters governing the solution. First, a forcing function ε(t) is proposed for the case when dynamic loading is applied under displacement control,

ε(t)=ε02(1cos2πft),
(51)

where ε0 is the peak-to-peak strain amplitude and f is the loading frequency. When solving in the non-dimensional domain, this equation reduces to

ε(t^)=ε02(1cos2πf^Rgt^),
(52)

where the non-dimensional loading frequency is given by

f^=r02HAkf.
(53)

Note that f is the loading frequency normalized by the characteristic frequency (the inverse of the biphasic theory's gel time constant [51]) for the flow of solvent in the porous matrix (e.g., [52, 53]). For the case when no dynamic loading is applied, it suffices to let f = 0.

Second-order finite difference schemes were employed for spatial discretization of the governing equations. The spatial domain 0 ≤ [r with circumflex] ≤ 1 was discretized into 200 equally-sized increments Δ[r with circumflex]. In certain challenging cases, such as those with a high Rg and low Rd, discretization with 2000 points demonstrated that while convergence was not quite attained with 200 increments, the solution slightly underestimated the actual converged solution. Due to limitations in computation time, 200 equally sized increments were used for all cases. The forcing function dε/dt appears in both governing equations of this problem; because the amplitude ε0πfRg of this function can be several orders of magnitude greater than unity under certain choices of the governing parameters (see below), the resulting system of nonlinear differential equations can be very stiff. Consequently, a variety of numerical schemes were explored to obtain a reliable solution, and for the purpose of verification the current study settled on two distinct methods. In the first method, a first-order backward difference scheme was employed for temporal discretization, producing an implicit scheme in the time domain. The resulting system of nonlinear difference equations were solved using the Newton-Raphson scheme with analytical Jacobian; band-diagonal solution algorithms were employed to take advantage of the structure of these equations for increased numerical efficiency. In the second method, the spatially discretized equations were arranged in the form of a system of first-order differential equations in time and solved using the Runge-Kutta-Verner fifth order and sixth order method (subroutine IVPRK of the IMSL library of subroutines, Visual Numerics, San Ramon, CA). In both methods, the time domain was discretized into uniform increments Δt; in the dynamic loading case, up to 2,000 time increments per loading cycle were necessary in some cases to achieve numerical convergence. For free diffusion (f = 0), the system of equations was not stiff and numerical convergence could be achieved with much fewer time iterations.

To achieve a solution in the non-dimensional domain requires the specification of Rg,Rd,νs,c^0f,c^f*,ε0andf^. To help determine characteristic values for these parameters, we look at some typical properties reported in the literature for articular cartilage and agarose gels, as representative examples for the current analysis. According to Table 2, the equilibrium aggregate modulus of articular cartilage in compression is typically in the range HA = 0.1–1.0 MPa, whereas the permeability is on the order of k = 10−16–10−15 m4/N.s, so that HAk varies approximately in the range 10−11–10−9 m2/s in articular cartilage. In 2–3% agarose and various chondrocyte-seeded engineered constructs, these parameters are in the range HA = 10–100 kPa, and k = 10−13–10−12 m4/N.s, thus HAk is in the range 10−10–10−8 m2/s. The diffusivity of various charged and neutral solutes in articular cartilage and agarose is provided in Table 1. In cartilage and agarose gels, the diffusivity D of small solutes, such as Na+, Cl and glucose, is on the order of 10−9 m2/s, whereas for larger solutes such as bovine serum albumin and 40kDa Dextran, the gel diffusivity is on the order of 10−11 m2/s. From these results, it is estimated that Rg is in the range 10−2–103 in cartilage and agarose gels (see Table 2). Typical values of the diffusivity Dwf of solutes in aqueous solutions are also provided in Table 1, along with representative values of the water volumetric fraction [var phi]f, from which it can be estimated that for neutral solutes Rd is typically in the range 0.1–1.0. Physiological loading frequencies for articular cartilage are typically in the range 0.01–1 Hz, and various basic science studies of the response of chondrocyte-seeded agarose gels have employed similar ranges (e.g., [18, 19, 22, 24, 28, 54, 55]). Given the representative values of HAk reported above, and typical radii r0 of 1–3 mm, a representative range for f is 101–106.

Table 2
Aggregate modulus (HA) and hydraulic permeability (k) for human and bovine articular cartilage and tissue engineered constructs. Values of Rg are given for a range of typical small (D = 1 × 10−9 m2/s) and large (D = 1 × 10−11 ...

Poisson's ratio, νs = 0, the initial solute concentration in the tissue, c^0f=0, and the external bath concentration (scaled by the partition factor), κfĉf* = 10−3, were kept constant in this analysis. The parameters which were varied were Rg = 1–100, Rd = 0.1–1, f = 10–1000 and ε0 = −0.20 – −0.01 (for transport with dynamic loading), or f = 0, ε0 = 0 (for loading-free transport).

Results

Representative results for selected cases are presented first, followed by aggregate results for all tested cases. The solute concentration, ĉf([r with circumflex],t), is shown as a function of the radial coordinate [r with circumflex] at various times t in Figure 2a for the case Rg = 100, Rd = 0.1, in the absence of dynamic loading. For this choice of governing parameters, the figure shows Fick's classical diffusion solution in cylindrical coordinates. At the radial edge of the sample, [r with circumflex] = 1, the solute concentration remains equal to ĉf* at all times, according to the boundary condition of Eq.(48)2. Inside the cylindrical sample, the solute concentration is initially equal to κfc^0f throughout the domain 0 ≤ [r with circumflex] < 1, but with increasing time the concentration eventually achieves an equilibrium value, which is the same as κfĉf*. When dynamic loading is applied, e.g., with frequency f = 1000 and strain amplitude ε0 = −0.20, the solute concentration distribution is generally oscillatory but can be plotted in a similar way, as shown in Figure 2b. For this case, the solute concentration distribution is significantly different than in the absence of loading. At early times of loading, in a region near the radial edge of the sample, solute concentration is nearly the same as ĉf*, while inside the tissue solute concentration remains below the external bath concentration ĉf*. As loading continues, solute concentration in the radial edge remains unchanged, while solute accumulates in the inner region of the tissue to levels greater than that of the boundary value. It is also evident, from the comparison of Figure 2a and Figure 2b, that solute concentration increases significantly faster in the dynamically loaded case. Since the solute concentration is generally non-uniform along the radial direction, its average value,

c^avgf(t^)=012πr^c^f(r^,t^)dr^012πr^dr^,
(54)

can be evaluated and plotted as a function of time to indicate the rate at which solute is taken up by the whole tissue, as shown in Figure 3 for select choices of governing parameters. When the Rg and Rd are at unity, it can be noted that although dynamic loading induces an oscillatory variation in average solute concentration, the mean response oscillates about the loading-free diffusion solution at low frequencies (f = 10) and shows only a slight elevation above the mean at higher frequencies (f = 1000). Thus dynamic loading at best slightly enhances solute transport for this particular choice of governing parameters. In contrast, when Rg is greater than unity and Rd is less than unity, increasing the loading frequency produces increasingly faster rises and increasingly higher values of the average solute concentration, as evident in Figure 4. To summarize the salient feature of these results, the steady-state value (averaged over the last loading cycle) of c^avgf(t^)/κfc^f* as t → ∞, and the time te when the average solute concentration has reached 63.2% of the external solute concentration, c^avgf(t^e)/κfc^f*=1e1, are plotted against Rg and loading frequency in Figure 5. For reference, in the absence of loading, regardless of the values of Rg or Rd, the corresponding values are c^avgf()/κfc^f*=1 and te = 0.111. From these results it is evident that for values of Rd below unity, higher Rg and higher loading frequencies may considerably enhance solute transport in this problem. The effect of the strain magnitude is similarly presented in Figure 6, where it is noted that increasing dynamic strain magnitudes also enhance solute transport, as long as Rd < 1, Rg > 1.

Figure 2
Solute concentration, ĉf([r with circumflex], t), at select time points for the case Rg = 100, Rd = 0.1, (a) in the absence of dynamic loading (ε0 = 0), and (b) when ε0 = −0.20 and f = 1000. ...
Figure 3
Average solute concentration normalized by the external bath concentration, c^avgf(t^)/κfc^f*, as a function of time, for various choices of governing parameters (Rg = 1, Rd = 1, ε0 = −0.20 and f = 1000; Rg = 1, ...
Figure 4
Average solute concentration normalized by the external bath concentration, c^avgf(t^)/κfc^f*, as a function of time, for various choices of governing parameters (Rg = 100, Rd = 0.1, ε0 = −0.20, and f = 0; Rg = ...
Figure 5
For Rd = 0.1, 0.5, and 1.0, ε0 = −0.20, (a) steady-state value of c^avgf(t^)/κfc^f* (averaged over a loading cycle) as t→∞, and (b) time te when c^avgf(t^e)/κfc^f*=1e1 ...
Figure 6
For Rd = 0.1 and Rg = 100, and f = 100, (a) transient value of c^avgf(t^)/κfc^f* versus t with increasing strain magnitude (ε0 = 0 to −0.20) and (b) steady-state values of c^avgf(t^)/κ ...

To help understand the solute transport mechanism, the molar flux of solute relative to the solid phase is plotted in Figure 7, for representative loading-free and dynamically loaded cases. In the loading-free case (f = 0), the solute transport is always directed into the tissue as indicated by its negative values, and is greatest in magnitude in the early time response, eventually decreasing to zero at equilibrium. In the dynamically loaded cases (f = 1000), the solute transport is directed into the tissue in the very early time response, with a magnitude much greater than the loading-free case; and the first cycle is seemingly independent of Rd. In subsequent cycles the solute molar flux relative to the solid phase becomes oscillatory, alternating between influx and efflux from the tissue. Interestingly, with Rd = 1, solute molar flux is nearly balanced, with only a slight movement into the tissue while with Rd = 0.1 solute flux is directed into the tissue over nearly the entire loading cycle. In both cases, the magnitude of the solute flux is greatest near the periphery of the cylindrical sample, decreasing to zero at its center.

Figure 7
Solute molar flux relative to solid phase, c^f(v^rfv^rs), at [r with circumflex] = 1 over the first five loading cycles for the case Rg = 100, Rd = 1, in the absence of dynamic loading (ε0 = 0), and when ε0 = −0.20, f ...

Figure 8 shows a representative plot of the radial displacement ûr([r with circumflex] = 1,t) and the axial load Ŵ(t) under loading-free conditions, with Rg = 100 and Rd = 0.1. The instantaneous influx of solute at time t = 0 causes a small, transient, drag induced compressive strain in the radial direction on the order of −0.08% (Figure 8a). This radial compression in turn generates a small transient axial load Ŵ, which reaches a peak of 0.0025 (Figure 8b). For a sample piece of cartilage, with HA = 0.5 MPa, r0 = 0.003 m, where W=W^HAr02, the peak axial load corresponds to 0.0113 N.

Figure 8
(a) Radial displacement ûr ([r with circumflex] = 1,t) and (b) axial load Ŵ(t) over time with Rd = 0.1, Rg = 100, in the absence of dynamic loading (f = 0).

Discussion

The objective of this study was to examine theoretically the effect of dynamic loading of a biological tissue or gel on solute transport, and to identify how the various non-dimensional parameters governing this problem affect its outcome. Using the theory of mixtures has provided a convenient framework to combine transport phenomena with porous media models of biological tissues. Indeed, under the proper limiting conditions, the general governing equations derived in this study reduce to Fick's laws of diffusion [44] or the linear isotropic biphasic theory of Mow and co-workers [34, 39]. A direct outcome of the analysis has been the formulation of familiar dimensional and non-dimensional parameters which govern the response. In particular, the analysis is able to distinguish between the solute diffusivity in the gel (D) and its diffusivity in the solvent (Dwf) by taking into account frictional effects between solute and solvent (fwf) as well as solute and gel solid matrix (fsf). It is interesting that the related non-dimensional parameter, Rd = D/[var phi]w Dwf, is scaled by the solvent volumetric fraction in the gel. If it is assumed that Rd has an upper bound of unity for neutral solutes in a neutral gel, this relation suggests (but does not prove) that even under ideal conditions (Rd = 1), the diffusivity D of a neutral solute in a neutral gel will at most equal the product of the gel porosity times the diffusivity in free solution, [var phi]w Dwf. This relationship implies that under ideal conditions solute diffusivity in the gel is directly proportional to the solvent volumetric fraction, a relationship which remains self-consistent over the entire range 0 ≤ [var phi]w ≤ 1. In general however, steric exclusion limits the volumetric fraction of accessible solvent so that the actual gel diffusivity D is smaller than the ideal gel diffusivity [var phi]w Dwf. As noted in Table 1, the experimental value of Rd decreases with increasing molecular weight, and only small solutes can achieve values that approach unity.

The theoretical analysis also produces a non-dimensional parameter, Rg = HAk/D, for the solute in the gel. This parameter can be identified with the ratio D(H)¯/D described by Urban et al. [14] and O'Hara et al. [12], as governing solute transport in a dynamically loaded tissue. In their expression, D(H)¯ is a 'fluid transport coefficient' [17], which these authors described as being dependent on the difference between the swelling and applied pressures and on the hydraulic permeability of the tissue [14]. The present study formulates this coefficient as the product of the familiar parameters HA (aggregate modulus) and k (hydraulic permeability) of the biphasic theory. The non-dimensional parameter Rg, which is the ratio of the characteristic diffusive velocities of the solvent and solute in the gel, can be related to the non-dimensional Peclet number, which is the ratio of the characteristic convective velocity V of the gel to the diffusive velocity D/r0 of the solute in the gel, Pe = Vr0/D. This relationship is given by Pe = Rg · Rh, where the non-dimensional parameter Rh = Vr0/HAk is the ratio of the convective velocity of the gel to the diffusive velocity of the solvent in the gel. Rh was first described by Lai and Mow [56, 57] when analyzing the problem of a load sliding with a velocity V on the surface of a biphasic articular layer, and has since been used in other related studies [58, 59]. It can be interpreted as a Peclet number for interstitial fluid diffusivity, in contrast to the more common usage of the Peclet number for solute diffusivity as described above, or for heat diffusivity. We now note that the dynamic compression of the gel, as provided in Eq.(51), induces a convective gel velocity of characteristic magnitude V = ε0r0f, so that the non-dimensional parameter f appearing in Eq.(53) is equal to Rh0. Hence the product fRg appearing in Eq.(52) is equal to Pe0, which indicates that the Peclet number also governs this problem, consistent with the electroosmotic study of Garcia et al. [13].

It is found that Rd and Rg have very little effect on the loading-free diffusion response (f = 0, Figure 5). Thus, despite the fact that the governing equations for loading-free diffusion in Eqs. (39)(40) are generally different from Fick's laws (when Rd ≠ 1), the gel deformation resulting from frictional drag between solute and gel solid matrix (Figure 8) does not significantly alter solute transport in the absence of dynamic loading. This result is encouraging because it suggests that in the absence of loading it is entirely consistent to use Fick's laws for the experimental determination of the diffusivity of a solute in a deformable soft tissue or gel [26, 10, 23] despite the fact that these laws do not fully describe transport within a deformable matrix. For loading-free diffusion, the equilibrium response yields a uniform solute concentration across the cylindrical gel (Figure 2a, ,4),4), and ĉavg(∞)/κfĉf* = 1. Because all non-dimensional concentrations are normalized with the same factor (Eq.(41)), this ratio is equal to the ratio of dimensional concentrations, cavg(∞)/κfcf* = 1. The non-dimensional characteristic time required to reach equilibrium in loading-free diffusion remains essentially constant at te = 0.111 over the range of Rd and Rg analyzed in this study; in dimensional form, the characteristic time is te=t^er02/D and is thus dependent on specimen size and gel diffusivity (a result similar in form to [60]).

In contrast to loading-free diffusion, under dynamic loading, the values of Rd, Rg, and f have a great influence on solute transport into the tissue. As Rd decreases from unity and Rg increases from unity, and with increasing non-dimensional loading frequency f, the equilibrium partition factor increases relative to that under loading-free diffusion, i.e., ĉavg(∞)/κfĉf* > 1 (Figure 5). Similarly, the non-dimensional time te correspondingly decreases (Figure 5). A small Rd and large Rg are representative of higher molecular weight solutes which typically have a smaller gel diffusivity D (everything else assumed constant). To illustrate this phenomenon, consider the transport of Gd-ovalbumin (MW= 45 kDa, D = 0.4×10−11m2/s, Dwf = 5.2×10−11m2/s, [11]) in a 5 mm disk (r0 = 2.5 mm) of articular cartilage (HA = 0.5 MPa, k = 10−15m4/N.s, [var phi]w = 0.8) loaded dynamically at a frequency f = 0.1 Hz with a strain amplitude ε0 = −0.2. These parameters produce Rd = 0.096, Rg = 125 and f = 1,250, and according to Figure 5 the solute concentration at steady-state would be cavg (∞)/κfcf* ~ 3 and the characteristic time te would be approximately 2.2 hrs. In contrast, under loading-free diffusion, cavg(∞)/κfcf* = 1 and te ~ 48 hrs (keeping in mind that diffusion is occurring only through the lateral surface of the cylindrical sample). It is important to emphasize that the results of this analysis do not necessarily suggest that the solute concentration inside the tissue will exceed the external concentration, because of the scaling factor introduced by the partition coefficient κf. In general, the larger the molecular weight of the solute, the smaller the partition coefficient will be. For example, if the external concentration of Gd-ovalbumin is cf* = 100 µM and its partition factor under loading-free diffusion is κf = 0.2, the steady-state concentration will be cavg(∞) = 20 µM without loading and cavg(∞)~60 µM with loading. It can be concluded that these theoretical predictions support the experimental finding [12, 14, 24] that dynamic loading enhances the transport of larger molecular weight solutes relative to transport under loading-free diffusion.

Conversely, when 0.5 ≤ Rd ≤ 1, and Rg ≤ 10, as is more typical of smaller solutes, there is less significant enhancement in the equilibrium partition factor relative to loading-free diffusion (Figure 5), though te does decrease with increasing Rg and f (Figure 5). This case can be illustrated with the transport of glucose (MW = 180 Da, D = 0.31×10−9m2/s, Dwf = 0.53×10−9m2/s, [7]) in cartilage (same parameters as above); then Rd = 0.73, Rg = 1.6 and f = 1,250, and under dynamic loading, cavg (∞)/κfcf* ~ 1.1; therefore, there is only a ~10% increase in average equilibrium solute concentration relative to loading-free diffusion in this case. The corresponding characteristic time te is approximately 28 min. with dynamic loading and 37 min. without loading, showing a more significant influence of dynamic loading on the initial transient response of solute transport than on the steady-state response. These predictions are also consistent with experimental findings described in the literature [12, 14].

To understand why dynamic loading can enhance solute transport in a hydrated gel or soft tissue, it is necessary to examine the governing equations along with their numerical solutions. The molar flux of solute relative to the solid phase is given by Eq.(46), which shows solute transport is driven by the concentration gradient ([partial differential]ĉf/[partial differential][r with circumflex]) as in classical diffusion, as well as an additional term contributed by the dynamic loading of the gel solid matrix, Rdĉf([partial differential]ûr/([partial differential]t + ([r with circumflex]/2)dε/dt); the net transport into or out of the tissue depends on these quantities at the circumferential edge of the cylindical sample ([r with circumflex] = 1). In the detailed view of Figure 7, the molar flux of solute relative to the solid is displayed for the first five cycles of loading, showing the loading-free diffusion case, and two comparable dynamically loaded cases which differ only in the value of Rd. In all cases, theory predicts that the solute flux is infinite at t = 0+ because of the initial jump in concentration between inside and outside; because of the finite mesh size used in the numerical analysis, the initial value of the flux is finite but demonstrates an asymptotic trend toward infinity. From these plots, it is immediately apparent that the flux in the presence of dynamic loading achieves greater magnitudes than under loading-free diffusion, over certain portions of the loading cycles. The primary reason for this behavior is that during dynamic loading the solute concentration gradient at [r with circumflex] = 1 oscillates from zero to a peak positive value which exceeds the peak concentration gradient under loading-free diffusion. Over time, this produces a faster rise in average solute concentration in the loaded cases than in the absence of loading. Interestingly, with Rd = 1, the flux is positive (solute efflux from the sample) over a significant portion of a loading cycle, whereas with Rd = 0.1, it remains mostly negative. Though not shown in the figure, the difference between these two loaded cases results from the contribution of the oscillatory term Rdcf ([partial differential]ûr/([partial differential]t+([r with circumflex]/2)dε/dt) to the flux, which is proportionally more significant with higher values of Rd. As a result, cases with small values of Rd produce a larger solute concentration at steady state because of the increased net influx of solute; this increase represents an effective increase in the solubility of the solute in the gel.

Based on the results of this study, it may be possible to estimate whether dynamic loading of tissue engineered constructs produces enhanced transport of growth factors that are present in or added to culture media. In our previous studies [29, 30], it was found that the equilibrium aggregate modulus of chondrocyte-seeded agarose gels increases from HA~5 kPa to HA~100 kPa over 28 days, when the 6.76 mm-diameter constructs are dynamically loaded in unconfined compression with f = 1 Hz and ε0 = −0.1 for one-hour-on, one-hour-off, thrice daily. Though not measured directly in those studies, the permeability of the same type of agarose and chondrocyte-seeded agarose gels has been shown to range from k ~ 10−12 m4/N.s to k ~ 10−13 m4/N.s [6163]. From the studies of Schneiderman et al. [8], the diffusivity of insulin-like growth factor (IGF-I) is estimated at D = 2.2×10−11 m2/s in articular cartilage. Combining these factors together suggests that during tissue elaboration, Rg increases from 250 to 500 and f decreases from 1,250 to 625; if Rd is estimated at <0.5, these results suggest that solute transport is enhanced by dynamic loading throughout the tissue growth process. Consequently, it is plausible that one of the beneficial effects of dynamic loading on cartilage tissue engineering is the enhanced transport of large solutes, such as growth factors, into the tissue. Based on the various examples examined in this study, there is a wide range of conditions which may enhance solute transport in gels and tissues; for our specific cartilage tissue engineering studies [29, 30], it appears that continuous dynamic loading at 1 Hz over a period of approximately 3 hours would be more appropriate to enhance growth factor adsorption into the constructs.

This novel application of mixture theory to describe solute transport in dynamically loaded porous permeable gels may give new insight into the possible mechanisms by which physical stimuli modulate tissue response and enhance tissue development. The limitations of this first theoretical formulation, however, are in considering both the solute and gel in which it diffuses as neutrally charged. Serum growth factors, such as those from the transforming growth factor (TGF) and insulin like growth factor (IGF)/somatomedin family have isoelectric points that range from acidic to basic [6472]. Molecules that are neutral in a physiologic environment will move independent of external charge, while those which carry a charge will be effected by the fixed charges of the matrix and the molecules around them [4, 73]. Thus the incorporation of fixed and movable charges into this mixture theory formulation may broaden its general application. Incorporating charge effects would require modeling the interstitial fluid with at least three species (a neutral solvent phase, an anionic solute phase and a cationic solute phase) [36, 37, 7476], instead of the two employed in this study. Analyses of these more complex models may be easier to interpret given the results of the current study.

In addition to charge effects, the inhomogeneity of material properties in both cartilage and tissue engineered constructs may influence the local distribution of solute transport [77]. Modeling the tension-compression nonlinearity of cartilage, as performed in our earlier studies [78, 79], may also represent an additional refinement of this analysis. Finally, the model proposed herein does not take account of consumption of solute molecules (either by cellular metabolism or binding to the extracellular matrix; e.g., [8083]), which may influence both transient and steady state solute concentration throughout the gel matrix. If solute binding or consumption is to be considered, equations for chemical reactions must be incorporated, and the balance of mass equations must be amended to incorporate mass generation (source or sink) terms. Incorporating the above modifications into this theory will increase its complexity, but also its utility both in understanding natural mechanotransduction mechanisms in native articular cartilage, as well its ability to predict the optimum mechanical environment for the in vitro functional tissue engineering of replacement tissues. This theoretical framework may be further extended to examine the influence of deformational loading on solute transport when engineered constructs are placed in situ, in which case transport across the articular surface may be of greater relevance.

This analysis also assumes that the volume fractions of the gel solid matrix and the solvent, the hydraulic permeability, and the solute diffusion coefficient remain constant under the small deformations being applied to the gel. It is possible to incorporate the higher order effect of deformation on these parameters using the conservation of mass relation for the volume fractions [36],

φwφ0w+(1φ0w)trE,φsφ0s(1trE),
(55)

where φ0wandφ0s are the solvent and solid volume fractions in the absence of deformation. Similarly, constitutive relations can be formulated for the functional dependence of k and D on tr E [10, 45, 84, 85], for example,

k=k0eMktrE,D=D0eMDtrE,
(56)

where k0 and D0 are the hydraulic permeability and solute gel diffusion coefficient in the absence of strain, and Mk,MD are non-dimensional coefficients.

In summary, this study presents an original theoretical analysis on the effect of dynamic loading on neutral solute transport in a neutral gel. It identifies the non-dimensional parameter Rd as an important governing parameter for this problem. It formulates the non-dimensional number Rg = HAk/D in terms of the aggregate modulus, hydraulic permeability, and solute diffusivity in the gel. It establishes that Rg, Rd, f and ε0 together regulate the response of solute transport to dynamic loading in a gel, and that variations in a single parameter are insufficient to predict the full spectrum of the response. It identifies the ranges of these non-dimensional governing parameters over which dynamic loading will enhance solute transport, and the implications for transport in cartilage and cartilage tissue engineering.

As a final remark, the theoretical results of this study are not uniquely relevant to natural and tissue engineered articular cartilage but apply to any porous hydrated tissue or gel. Dynamic loading is inherent to a variety of biological tissues (e.g., muscle, tendon, ligaments, cardiovascular system, etc.). Its effects need not be limited to extracellular matrix but may be equally relevant to transport at the cellular level [86]. It is therefore plausible to hypothesize that dynamic loading may serve to enhance solute transport in a variety of physiological processes, particularly in the case of larger molecular weight solutes. Such effects may be investigated using the theoretical framework provided in this study.

Acknowledgments

This study was supported with funds from the National Institutes of Health (AR46532, AR46568, AR43628) and a pre-doctoral fellowship from the Whitaker Foundation.

References

1. McKibbin B, Maroudas A. Nutrition and metabolism. In: Freeman MAR, editor. Adult articular cartilage. New York: Grune & Stratton; 1974. pp. 461–486.
2. Maroudas A. Physicochemical properties of cartilage in the light of ion exchange theory. Biophys J. 1968;8(5):575–595. [PubMed]
3. Maroudas A. Biophysical chemistry of cartilaginous tissues with special reference to solute and fluid transport. Biorheology. 1975;12(3–4):233–248. [PubMed]
4. Torzilli PA, Arduino JM, Gregory JD, Bansal M. Effect of proteoglycan removal on solute mobility in articular cartilage. J Biomech. 1997;30(9):895–902. [PubMed]
5. Torzilli PA, Adams TC, Mis RJ. Transient solute diffusion in articular cartilage. J Biomech. 1987;20(2):203–214. [PubMed]
6. Torzilli PA. Effects of temperature, concentration and articular surface removal on transient solute diffusion in articular cartilage. Med Biol Eng Comput. 1993;31 Suppl:S93–S98. [PubMed]
7. Burstein D, Gray ML, Hartman AL, Gipe R, Foy BD. Diffusion of small solutes in cartilage as measured by nuclear magnetic resonance (NMR) spectroscopy and imaging. J Orthop Res. 1993;11(4):465–478. [PubMed]
8. Schneiderman R, Snir E, Popper O, Hiss J, Stein H, Maroudas A. Insulin-like growth factor-I and its complexes in normal human articular cartilage: studies of partition and diffusion. Arch Biochem Biophys. 1995;324(1):159–172. [PubMed]
9. Potter K, Spencer RG, McFarland EW. Magnetic resonance microscopy studies of cation diffusion in cartilage. Biochim Biophys Acta. 1997;1334(2–3):129–139. [PubMed]
10. Quinn TM, Kocian P, Meister JJ. Static compression is associated with decreased diffusivity of dextrans in cartilage explants. Arch Biochem Biophys. 2000;384(2):327–334. [PubMed]
11. Foy BD, Blake J. Diffusion of paramagnetically labeled proteins in cartilage: enhancement of the 1-D NMR imaging technique. J Magn Reson. 2001;148(1):126–134. [PubMed]
12. O'Hara BP, Urban JP, Maroudas A. Influence of cyclic loading on the nutrition of articular cartilage. Ann Rheum Dis. 1990;49(7):536–539. [PMC free article] [PubMed]
13. Garcia AM, Frank EH, Grimshaw PE, Grodzinsky AJ. Contributions of fluid convection and electrical migration to transport in cartilage: relevance to loading. Arch Biochem Biophys. 1996;333(2):317–325. [PubMed]
14. Urban JP, Holm S, Maroudas A, Nachemson A. Nutrition of the intervertebral disc: effect of fluid flow on solute transport. Clin Orthop. 1982;170:296–302. [PubMed]
15. Katz MM, Hargens AR, Garfin SR. Intervertebral disc nutrition. Diffusion versus convection. Clin Orthop. 1986;210:243–245. [PubMed]
16. Garcia AM, Lark MW, Trippel SB, Grodzinsky AJ. Transport of tissue inhibitor of metalloproteinases-1 through cartilage: contributions of fluid flow and electrical migration. J Orthop Res. 1998;16(6):734–742. [PubMed]
17. Fatt I, Goldstick TK. Dynamics of water transport in swelling membranes. J. Colloid Sci. 1965;20:962–989. [PubMed]
18. Sah RL, Kim YJ, Doong JY, Grodzinsky AJ, Plaas AH, Sandy JD. Biosynthetic response of cartilage explants to dynamic compression. J Orthop Res. 1989;7(5):619–636. [PubMed]
19. Kim YJ, Sah RL, Grodzinsky AJ, Plaas AH, Sandy JD. Mechanical regulation of cartilage biosynthetic behavior: physical stimuli. Arch Biochem Biophys. 1994;311(1):1–12. [PubMed]
20. Palmoski MJ, Brandt KD. Effects of static and cyclic compressive loading on articular cartilage plugs in vitro. Arthritis Rheum. 1984;27(6):675–681. [PubMed]
21. Guilak F, Sah RL, Setton LA. Physical regulation of cartilage metabolism. In: Mow VC, Hayes WC, editors. Basic orthopaedic biomechanics. Philadelphia: Lippincott-Raven; 1997. pp. 179–207.
22. Guilak F, Meyer BC, Ratcliffe A, Mow VC. The effects of matrix compression on proteoglycan metabolism in articular cartilage explants. Osteoarthritis Cartilage. 1994;2(2):91–101. [PubMed]
23. Quinn TM, Morel V, Meister JJ. Static compression of articular cartilage can reduce solute diffusivity and partitioning: implications for the chondrocyte biological response. J Biomech. 2001;34(11):1463–1469. [PubMed]
24. Bonassar LJ, Grodzinsky AJ, Frank EH, Davila SG, Bhaktav NR, Trippel SB. The effect of dynamic compression on the response of articular cartilage to insulin-like growth factor-I. J Orthop Res. 2001;19(1):11–17. [PubMed]
25. Kuettner KE. Biochemistry of articular cartilage in health and disease. Clin Biochem. 1992;25(3):155–163. [PubMed]
26. DiMicco MA, Sah RL. Dependence of cartilage matrix composition on biosynthesis, diffusion, and reaction. Tranport in Porous Media. 2003;50(1–2):57–73.
27. Freed LE, Vunjak-Novakovic G, Langer R. Cultivation of cell-polymer cartilage implants in bioreactors. J Cell Biochem. 1993;51(3):257–264. [PubMed]
28. Buschmann MD, Gluzband YA, Grodzinsky AJ, Hunziker EB. Mechanical compression modulates matrix biosynthesis in chondrocyte/agarose culture. J Cell Sci. 1995;108(Pt 4):1497–1508. [PubMed]
29. Mauck RL, Soltz MA, Wang CC, Wong DD, Chao PH, Valhmu WB, Hung CT, Ateshian GA. Functional tissue engineering of articular cartilage through dynamic loading of chondrocyte-seeded agarose gels. J Biomech Eng. 2000;122(3):252–260. [PubMed]
30. Mauck RL, Nicoll SB, Seyhan SL, Ateshian GA, Hung CT. Synergistic action of growth factors and dynamic loading for articular cartilage tissue engineering. Tiss Eng. 2003 in press. [PubMed]
31. Mills N. Incompressible mixture of Newtonian fluids. Int J Engng Sci. 1966;4:97–112.
32. Craine RE. Oscillations of a plate in a binary mixture of incompressible newtonian fluids. Int J Engng Sci. 1971;9(12):1177–1192.
33. Bowen RM. Incompressible porous media models by use of the theory of mixtures. Int J Engng Sci. 1980;18(9):1129–1148.
34. Mow VC, Kuei SC, Lai WM, Armstrong CG. Biphasic creep and stress relaxation of articular cartilage in compression? Theory and experiments. J Biomech Eng. 1980;102(1):73–84. [PubMed]
35. Frank EH, Grodzinsky AJ. Cartilage electromechanics--II. A continuum model of cartilage electrokinetics and correlation with experiments. J Biomech. 1987;20(6):629–639. [PubMed]
36. Lai WM, Hou JS, Mow VC. A triphasic theory for the swelling and deformation behaviors of articular cartilage. J Biomech Eng. 1991;113(3):245–258. [PubMed]
37. Huyghe JM, Janssen JD. Quadriphasic mechanics of swelling incompressible porous media. Int J Engng Sci. 1997;35(8):793–802.
38. Gu WY, Lai WM, Mow VC. A mixture theory for charged-hydrated soft tissues containing multi-electrolytes: passive transport and swelling behaviors. J Biomech Eng. 1998;120(2):169–180. [PubMed]
39. Armstrong CG, Lai WM, Mow VC. An analysis of the unconfined compression of articular cartilage. J Biomech Eng. 1984;106(2):165–173. [PubMed]
40. Tinoco I. Physical chemistry: principles and applications in biological sciences. Upper Saddle River, New Jersey: Prentice Hall; 2002.
41. Robinson RA, Stokes RH. Electrolyte solutions; the measurement and interpretation of conductance, chemical potential and diffusion in solutions of simple electrolytes. New York: Academic Press; 1955.
42. Torzilli PA, Askari E, Jenkins JT. Water content and solute diffusion properties in articular cartilage. In: Mow VC, Ratcliffe A, Woo SLY, editors. Biomechanics of diarthrodial joints. New York: Springer-Verlag; 1990. pp. 363–390.
43. Fournier RL. Basic transport phenomena in biomedical engineering. Philadelphia: Taylor & Francis; 1999. Solute diffusion within heterogenous media; pp. 28–32.
44. Katchalsky A, Curran PF. Nonequilibrium thermodynamics in biophysics. Cambridge: Harvard University Press; 1975. Isothermal diffusion and sedimentation; pp. 98–102.
45. Lai WM, Mow VC. Drag-induced compression of articular cartilage during a permeation experiment. Biorheology. 1980;17(1–2):111–123. [PubMed]
46. Hou JS, Holmes MH, Lai WM, Mow VC. Boundary conditions at the cartilage-synovial fluid interface for joint lubrication and theoretical verifications. J Biomech Eng. 1989;111(1):78–87. [PubMed]
47. Sun DN, Gu WY, Guo XE, Lai WM, Mow VC. A mixed finite element formulation of triphasic mechano-electrochemical theory for charged, hydrated soft tissues. Int J Num Meth Eng. 1999;45(10):1375–1402.
48. Ateshian GA, Soltz MA, Mauck RL, Basalo IM, Hung CT, Lai WM. The role of osmotic pressure and tension-compression nonlinearity in the frictional response of articular cartilage. Transport in Porous Media. 2003;50:5–33.
49. Van Holde KE, Johnson WC, Ho PS. Principles of physical biochemistry. Saddle River: Prentice-Hall; 1998. Thermodynamics of transport processes; p. 574.
50. Deen WM. Hindered transport of large molecules in liquid-filled pores. AIChE J. 1987;33:1409–1425.
51. Holmes MH, Lai WM, Mow VC. Singular perturbation analysis of the nonlinear, flow-dependent compressive stress relaxation behavior of articular cartilage. J Biomech Eng. 1985;107(3):206–218. [PubMed]
52. Lee RC, Frank EH, Grodzinsky AJ, Roylance DK. Oscillatory compressional behavior of articular cartilage and its associated electromechanical properties. J Biomech Eng. 1981;103(4):280–292. [PubMed]
53. Soltz MA, Ateshian GA. Interstitial fluid pressurization during confined compression cyclical loading of articular cartilage. Ann Biomed Eng. 2000;28(2):150–159. [PubMed]
54. Grodzinsky AJ, Levenston ME, Jin M, Frank EH. Cartilage tissue remodeling in response to mechanical forces. Annu Rev Biomed Eng. 2000;2:691–713. [PubMed]
55. Lee DA, Bader DL. Compressive strains at physiological frequencies influence the metabolism of chondrocytes seeded in agarose. J Orthop Res. 1997;15(2):181–188. [PubMed]
56. Mow VC, Lai WM. Recent developments in synovial joint biomechanics. SIAM Review. 1980;22:275–317.
57. Lai WM, Mow VC. Flow fields in a single layer model of articular cartilage created by a sliding load. ASME Adv Bioengng. 1979:101–104.
58. Ateshian GA, Wang H. A theoretical solution for the frictionless rolling contact of cylindrical biphasic articular cartilage layers. J Biomech. 1995;28(11):1341–1355. [PubMed]
59. Ateshian GA, Wang X. Sliding tractions on a deformable porous layer. Journal of Tribology. 1998;120(1):89–96.
60. Bonassar LJ, Grodzinsky AJ, Srinivasan A, Davila SG, Trippel SB. Mechanical and physicochemical regulation of the action of insulin-like growth factor-I on articular cartilage. Arch Biochem Biophys. 2000;379(1):57–63. [PubMed]
61. Buschmann MD, Gluzband YA, Grodzinsky AJ, Kimura JH, Hunziker EB. Chondrocytes in agarose culture synthesize a mechanically functional extracellular matrix. J Orthop Res. 1992;10(6):745–758. [PubMed]
62. Soltz MA, Stankiewicz A, Mauck RL, Hung CT, Ateshian GA. Direct hydraulic permeability measurements of agarose hydrogels used as cell scaffolds. ASME Adv Bioeng. 1999;43:229–230.
63. Andarawis NA, Seyhan SL, Mauck RL, Soltz MA, Ateshian GA, Hung CT. A novel permeation device for hydrogels and soft tissues; Proceedings ASME IMECE, paper no. 23149; 2001.
64. Perdue JF. Chemistry, structure, and function of insulin-like growth factors and their receptors: a review. Can J Biochem Cell Biol. 1984;62(11):1237–1245. [PubMed]
65. Enberg G, Carlquist M, Jornvall H, Hall K. The characterization of somatomedin A, isolated by microcomputer-controlled chromatography, reveals an apparent identity to insulin-like growth factor 1. Eur J Biochem. 1984;143(1):117–124. [PubMed]
66. Kuffer AD, Herington AC. Proteolytic conversion of insulin-like growth factors to an acidic form(s) Biochem J. 1984;223(1):97–103. [PubMed]
67. Herington AC, Kuffer AD. Insulin-like growth factor characteristics of an acidic non-suppressible insulin-like activity. Biochem J. 1984;223(1):89–96. [PubMed]
68. Herington AC, Cornell HJ, Kuffer AD. Recent advances in the biochemistry and physiology of the insulin-like growth factor/somatomedin family. Int J Biochem. 1983;15(10):1201–1210. [PubMed]
69. Kim MK, Warren TC, Kimball ES. Purification and characterization of a low molecular weight transforming growth factor from the urine of melanoma patients. J Biol Chem. 1985;260(16):9237–9243. [PubMed]
70. Yamaoka K, Hirai R, Tsugita A, Mitsui H. The purification of an acid- and heat-labile transforming growth factor from an avian sarcoma virus-transformed rat cell line. J Cell Physiol. 1984;119(3):307–314. [PubMed]
71. Wang CC-B, Soltz MA, Mauck RL, Valhmu WB, Ateshian GA, Hung CT. Comparison of the equilibrium axial strain distribution in articular cartilage explants and cell-seeded alginate disks under unconfined compression. Trans Orthop Res Soc. 2000;25:131.
72. Maroudas A. Physicochemical properties of articular cartilage. In: Freeman MAR, editor. Adult Articular Cartilage. Kent: Pitman Medical; 1979. pp. 215–290.
73. Johnson EM, Berk DA, Jain RK, Deen WM. Diffusion and partitioning of proteins in charged agarose gels. Biophys J. 1995;68(4):1561–1568. [PubMed]
74. Lai WM, Mow VC, Sun DD, Ateshian GA. On the electric potentials inside a charged soft hydrated biological tissue: streaming potential versus diffusion potential. J Biomech Eng. 2000;122(4):336–346. [PubMed]
75. Lai WM, Sun DD, Ateshian GA, Guo XE, Mow VC. Electrical signals for chondrocytes in cartilage. Biorheology. 2002;39(1–2):39–45. [PubMed]
76. Mow VC, Ateshian GA, Lai WM, Gu WY. Effects of fixed charges on the stress-relaxation behavior of hydrated soft tissues in a confined compression problem. Int J Solids Structures. 1998;35(34–35):4945–4962.
77. Wang CC, Hung CT, Mow VC. An analysis of the effects of depth-dependent aggregate modulus on articular cartilage stress-relaxation behavior in compression. J Biomech. 2001;34(1):75–84. [PubMed]
78. Soltz MA, Ateshian GA. A Conewise Linear Elasticity mixture model for the analysis of tension-compression nonlinearity in articular cartilage. J Biomech Eng. 2000;122(6):576–586. [PMC free article] [PubMed]
79. Huang CY, Mow VC, Ateshian GA. The role of flow-independent viscoelasticity in the biphasic tensile and compressive responses of articular cartilage. J Biomech Eng. 2001;123(5):410–417. [PubMed]
80. Horner HA, Urban JP. 2001 Volvo Award Winner in Basic Science Studies: Effect of nutrient supply on the viability of cells from the nucleus pulposus of the intervertebral disc. Spine. 2001;26(23):2543–2549. [PubMed]
81. Bhakta NR, Garcia AM, Frank EH, Grodzinsky AJ, Morales TI. The insulin-like growth factors (IGFs) I and II bind to articular cartilage via the IGF-binding proteins. J Biol Chem. 2000;275(8):5860–5866. [PubMed]
82. Pedrozo HA, Schwartz Z, Gomez R, Ornoy A, Xin-Sheng W, Dallas SL, Bonewald LF, Dean DD, Boyan BD. Growth plate chondrocytes store latent transforming growth factor (TGF)-beta 1 in their matrix through latent TGF-beta 1 binding protein-1. J Cell Physiol. 1998;177(2):343–354. [PubMed]
83. Chintala SK, Miller RR, McDevitt CA. Basic fibroblast growth factor binds to heparan sulfate in the extracellular matrix of rat growth plate chondrocytes. Arch Biochem Biophys. 1994;310(1):180–186. [PubMed]
84. Lai WM, Mow VC, Roth V. Effects of nonlinear strain-dependent permeability and rate of compression on the stress behavior of articular cartilage. J Biomech Eng. 1981;103(2):61–66. [PubMed]
85. Gu WY, Yao H, Huang CY, Cheung HS. New insight into deformation-dependent hydraulic permeability of gels and cartilage, and dynamic behavior of agarose gels in confined compression. J Biomech. 2003;36(4):593–598. [PubMed]
86. Luby-Phelps K, Castle PE, Taylor DL, Lanni F. Hindered diffusion of inert tracer particles in the cytoplasm of mouse 3T3 cells. Proc Natl Acad Sci U S A. 1987;84:4910–4913. [PubMed]
87. Allhands RV, Torzilli PA, Kallfelz FA. Measurement of diffusion of uncharged molecules in articular cartilage. Cornell Vet. 1984;74(2):111–123. [PubMed]
88. Roger P, Mattisson C, Axelsson A, Zacchi G. Use of holographic laser interferometry to study the diffusion of polymers in gels. Biotechnol Bioeng. 2000;69(6):654–663. [PubMed]
89. Johnson EM, Berk DA, Jain RK, Deen WM. Hindered diffusion in agarose gels: test of effective medium model. Biophys J. 1996;70(2):1017–1023. [PubMed]
90. Mow VC, Hou JS, Owens JM, Ratcliffe A. Biphasic and quasilinear viscoelastic theories for hydrated soft tissues. In: Mow VC, Ratcliffe A, Woo SLY, editors. Biomechanics of diarthrodial joints. New York: Springer-Verlag; 1990. pp. 215–260.
91. Williamson AK, Chen AC, Sah RL. Compressive properties and function-composition relationships of developing bovine articular cartilage. J Orthop Res. 2001;19(6):1113–1121. [PubMed]
92. Freed LE, Langer R, Martin I, Pellis NR, Vunjak-Novakovic G. Tissue engineering of cartilage in space. Proc Natl Acad Sci U S A. 1997;94(25):13885–13890. [PubMed]
93. Chang SC, Rowley JA, Tobias G, Genes NG, Roy AK, Mooney DJ, Vacanti CA, Bonassar LJ. Injection molding of chondrocyte/alginate constructs in the shape of facial implants. J Biomed Mater Res. 2001;55(4):503–511. [PubMed]