|Home | About | Journals | Submit | Contact Us | Français|
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 (), and the compressive strain amplitude (ε0). Results show that when Rg > 1, Rd < 1, and > 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 > 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.
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 . 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., [2–11], 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 , and by electroosmotic flow , 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 . Urban and co-workers have attributed the disparate outcome under dynamic loading to the magnitude of the ratio of the 'fluid transport coefficient'  to the solute diffusion coefficient in the tissue [12, 14], thus defining a non-dimensional governing parameter for this problem. Garcia et al.  determined that electroosmotic flow was regulated by the Peclet number, which is a different but related non-dimensional parameter.
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 [18–22]. 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 . 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 . 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 . 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 , 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 . 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 . 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 .” 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 . As observed in the biosynthetic response of cartilage explants in the study by Bonassar et al. , 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 . 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 [31–38] 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 .
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:
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
thus, concentrations and apparent densities can be interchanged using these relations.
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 [31–38]; using the notation of Gu et al. , the equations are given by
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, w, f, are given by
where we assume that the solute volume fraction is negligible, f 1. This assumption is justified for most circumstances as illustrated by the two following examples: For a small solute such as glucose with a typical concentration in media of 1g/L (f = 1/180 mole/L), the volume fraction according to Eq.(1) is given by , 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 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 λs,μs, infinitesimal strain tensor E = (gradu + gradT u)/2, and solid displacement u, similarly to the biphasic analysis of unconfined compression by Armstrong et al. . 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), 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 be the solvent volume accessible to the solute, then represents the volumetric fraction of accessible solvent, or the solubility of the solute in the gel (0 ≤ κf ≤ 1), and is the concentration of solute per volume of accessible solvent. The solute activity is then defined as af = γff/c0 = γfcf/κfc0 where γf is the solute activity coefficient, and the solute standard state  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
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 . Thus, it can be shown that dependent variables in this problem exhibit the dependencies p = p(r,t), ur = ur (r,t), uz/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
where HA = λs + 2μs. The unknowns in these equations are 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 at r = 0 yields
where the diffusion coefficients Dwf, Dsf and permeability k are defined by
The definition of Dwf, the diffusivity of the solute in the solvent, is consistent with the classical treatment of diffusion  as will be shown below. The definition of permeability k is the same as that presented by Lai et al.  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),
where we have substituted and assumed that 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,
To get the relative molar flux per unit mixture area, , multiply this equation by w. Using Eq.(20) in Eq.(14) similarly produces a differential equation which only depends on solute concentration and solid displacement,
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 .
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
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 [46–48], [[σrr]] = 0, [[w]] = 0, and [[f]] = 0, which lead to the relations
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 and ambient pressure p0:
If , we would expect a net transport of solute into the tissue over time, whereas the converse would occur when . Finally, the axial normal traction component and the total axial load also need to be evaluated,
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 (s = 0, w ≈ 1), D reduces to Dwf and and (t) drop out of Eq.(26) to yield Fick's first law of diffusion,
whereas Eq.(25) yields Fick's second law of diffusion (in cylindrical coordinates),
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 (s ≠ 0), then ε(t) = 0 since no solid matrix deformation can take place. In general, 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 , in which case Eq.(26) reduces to
while Eq.(25) yields
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 ((t) = 0), the governing equations are found to differ in general from Fick's classical laws, with Eqs. (25)–(27) reducing to
It is noteworthy that for the special case when D = 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).
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,
where the non-dimensional parameters Rg and Rd are defined as
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/wDwf represents the ratio of the characteristic gel diffusive velocity of the solute (D/r0) to its ideal characteristic gel diffusive velocity (wDwf/r0). The ratio D/Dwf is commonly reported in experimental studies of diffusivity (e.g., [43, 50], see Table 1).
The boundary and initial conditions and the axial normal traction and load, when using non-dimensional variables, remain essentially unchanged,
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.
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,
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
where the non-dimensional loading frequency is given by
Note that is the loading frequency normalized by the characteristic frequency (the inverse of the biphasic theory's gel time constant ) 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 = 0.
Second-order finite difference schemes were employed for spatial discretization of the governing equations. The spatial domain 0 ≤ ≤ 1 was discretized into 200 equally-sized increments Δ. 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ε/d appears in both governing equations of this problem; because the amplitude ε0πRg 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 Δ; 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 ( = 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 . 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 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 is 101–106.
Poisson's ratio, νs = 0, the initial solute concentration in the tissue, , 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, = 10–1000 and ε0 = −0.20 – −0.01 (for transport with dynamic loading), or = 0, ε0 = 0 (for loading-free transport).
Representative results for selected cases are presented first, followed by aggregate results for all tested cases. The solute concentration, ĉf(,), is shown as a function of the radial coordinate at various times 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, = 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 throughout the domain 0 ≤ < 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 = 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,
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 ( = 10) and shows only a slight elevation above the mean at higher frequencies ( = 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 as → ∞, and the time e when the average solute concentration has reached 63.2% of the external solute concentration, , 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 and e = 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.
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 ( = 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 ( = 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 8 shows a representative plot of the radial displacement ûr( = 1,) and the axial load Ŵ() under loading-free conditions, with Rg = 100 and Rd = 0.1. The instantaneous influx of solute at time = 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 , the peak axial load corresponds to 0.0113 N.
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  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/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, 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 ≤ 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 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 described by Urban et al.  and O'Hara et al. , as governing solute transport in a dynamically loaded tissue. In their expression, is a 'fluid transport coefficient' , which these authors described as being dependent on the difference between the swelling and applied pressures and on the hydraulic permeability of the tissue . 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 appearing in Eq.(53) is equal to Rh/ε0. Hence the product Rg appearing in Eq.(52) is equal to Pe/ε0, which indicates that the Peclet number also governs this problem, consistent with the electroosmotic study of Garcia et al. .
It is found that Rd and Rg have very little effect on the loading-free diffusion response ( = 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 [2–6, 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 e = 0.111 over the range of Rd and Rg analyzed in this study; in dimensional form, the characteristic time is and is thus dependent on specimen size and gel diffusivity (a result similar in form to ).
In contrast to loading-free diffusion, under dynamic loading, the values of Rd, Rg, and 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 , 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 e 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, ) in a 5 mm disk (r0 = 2.5 mm) of articular cartilage (HA = 0.5 MPa, k = 10−15m4/N.s, 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 = 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 e does decrease with increasing Rg and (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, ) in cartilage (same parameters as above); then Rd = 0.73, Rg = 1.6 and = 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 (ĉf/) as in classical diffusion, as well as an additional term contributed by the dynamic loading of the gel solid matrix, Rdĉf(ûr/( + (/2)dε/d); the net transport into or out of the tissue depends on these quantities at the circumferential edge of the cylindical sample ( = 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 = 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 = 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 (ûr/(+(/2)dε/d) 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 [61–63]. From the studies of Schneiderman et al. , 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 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 [64–72]. 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, 74–76], 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 . 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., [80–83]), 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 ,
where 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,
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, 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 . 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.
This study was supported with funds from the National Institutes of Health (AR46532, AR46568, AR43628) and a pre-doctoral fellowship from the Whitaker Foundation.