INTRODUCTION

Top Abstract Introduction Results & Discussion References

THE EXCHANGE OF WATER AND macromolecules from circulating plasma to the interstitial space is a basic phenomenon that regulates important biological processes (3). Several mathematical models have been developed in the past to describe the hydraulic permeability of the capillary wall, with the aim of identifying the physical forces responsible for this function. The majority of models proposed for the hydraulic permeability of extracellular space within the capillary wall, known as fiber matrix models, are based on the assumption that the extracellular space is filled with fibrous material.

Viscous flow through fibrous media (and through porous media in general) is usually described by Darcy's law, according to which the average fluid velocity vector U is proportional to the average pressure gradient P

(1)

where µ is the fluid viscosity, and K_D is the Darcy permeability of the medium. In the case of extracellular matrix, it is impossible to determine experimentally the dependence of K_D from the dimensions and the arrangement of the fibers, because of technical difficulties. In fact, local water flow rates and pressures, on a microscopic scale, cannot be directly measured and, in addition, ultrastructural techniques do not allow obtainment of reliable quantitative informations on the actual dimensions of extracellular protein fibers. Several mathematical expressions have been proposed in the literature to calculate K_D from geometrical parameters of a fiber matrix, such as the radius of the fibers and their length per unit volume. However, there are large differences in the values of K_D calculated using the available expressions. Furthermore, none of these expressions is based on direct solution of the local hydrodynamic flow inside a three-dimensional fiber matrix.

The aim of our study was then to investigate the hydraulic permeability properties of a specific fiber matrix model of basement membrane by numerical simulation of water flow across the fibers and to compare the results with the most popular analytical expressions reported in literature. Extracellular matrix has been usually modeled as a random array of fibers (3). However, recent ultrastructural studies of the basement membrane (GBM) of the glomerular capillary wall show a rather regular fibrillar organization with junctions connecting three or four short fibers (strands) to form a three-dimensional meshwork (15, 16, 21). We then considered an ordered array of cylindrical strands connected as shown in Fig. 1. We determined K_D by numerical simulation of water flow in a periodic element of the array, under a given pressure gradient. We calculated K_D assuming several values of fiber radius and length of the geometrical model and compared the results of the numerical analysis with several analytical expressions proposed in the literature. We also obtained a simple interpolating formula that relates K_D to the geometrical parameters of the fiber array. We discuss the applications of our results to the problem of water filtration through the GBM. Since recent studies on glomerular membrane permeability indicate that changes in K_D of the GBM are associated with development of renal diseases (7, 11-13, 29), our model could be used to derive more insight on possible changes in GBM ultrastructure that are responsible for loss of hydraulic permeability. Furthermore, our geometrical model can be applied to other fibrous media, and our results can be used to describe viscous flow in such media.

View larger version (37K): [in this window] [in a new window]   Fig. 1.   Schematic representation of the tetrahedral periodic array of cylinders. A: isometric view. B: top view. S1, S2, and S3 are symmetry planes.

	PROBLEM FORMULATION

Available studies based on high resolution electron microscopy (17, 18, 24, 30) suggest that the GBM is composed of a three-dimensional meshwork of fibrils with a polygonal structure with pores ranging from 4 to 6 nm in diameter. There are biochemical evidences (1) that collagen type IV macromolecules are flexible rods of about 400 nm in length and that they form the mechanical support of the GBM meshwork. Assembly of collagen IV molecules is due to formation of dimers (with the association of the their COOH-terminals) or tetramers (by association of the NH₂-terminal) (1). There are also evidences that collagen type IV molecules align laterally and twist around each other (33) to form a scaffold into which other components of the GBM are incorporated (such as laminin, heparan sulfate proteoglycans, and others) (31). To mimic the spatial arrangement of fibrils within the GBM, as shown by Moriya et al. (24), we reasoned that the most symmetric junction of four filaments is the tetrahedral structure in which any two adjacent filaments form the same angle (109.5°). We thus considered a tetrahedral periodic arrangement of cylinders of uniform radius (r) and length (L), as shown schematically in Fig. 1, where the periodic unit element consists in the hexagonal prism shown in Fig. 1A. We calculated the volume occupied by the solid () as

(2)

Assuming that the structure is oriented parallel to the endothelial surface, we considered that water flows through the meshwork along the z-axis (see Fig. 1A). With this assumption, planes S₁, S₂, and S₃ (see Fig. 1B) are symmetry planes. Therefore, flow analysis can be restricted to the volume delimited by these three planes. Within this volume the flow field is periodic along the z-axis; thus we considered the triangular prism represented in Fig. 2, composed of three identical elements in the z direction, with periodic boundaries in planes ₁ and ₄. Planes ₂ and ₃ are also periodic with ₁ but with a rotation of 120° around the z-axis. For this reason, we then furthered restricted the hydrodynamical domain to the part of the triangular prism delimited by the planes ₁ and ₂.

View larger version (22K): [in this window] [in a new window]   Fig. 2.   Computational domain for the 3-dimensional array. Planes 1, 2, 3, and 4 are periodic boundaries. Side surfaces of the triangular prism are symmetry boundaries.

The approach used to calculate K_D was to simulate numerically the flow field inside the described domain under a given average pressure gradient P directed along the z-axis. From the results of the numerical analysis, the average velocity U was calculated, and K_D was determined using Eq. 1. Numerical simulations for several values of , ranging from a very dilute ( = 0.005) to a very concentrated array ( = 0.7), were performed. For convenience, results are expressed in terms of the dimensionless permeability

(3)

The advantage of using K * instead of K_D is that the latter depends on two quantities (r and ), whereas K * depends only on the ratio r/L or, equivalently, on .

The local Reynolds number (Re) of water moving within the assumed meshwork of fibers was calculated as Re = Ur/µ, where is the water density, µ is the water viscosity, and U = |U|. For water filtration through the GBM, with r of the order of a few nanometers (17, 18, 24), Re is less than 10⁷. In this flow regimen the proportionality between U and P holds strictly; therefore K_D is a well-defined parameter (independent from Re). We then assumed Re = 0 and considered the stationary Stokes and continuity equations for an incompressible fluid

(4)

(5)

where u is the local velocity field, and p is the local pressure field. Input values of P were imposed by means of a boundary condition on the surface stress . The following equality holds

(6)

where the integral on the left-hand side extends to the surface of the simulation cell, with exclusion of the solid parts. In Eq. 9b V_p = 8/(9)L³ is the volume of the simulation cell, which is one-third of the entire periodic cell (the triangular prism delimited by ₁ and ₂ in Fig. 2). To impose a given P, it is sufficient to impose  in such a way that Eq. 6 is satisfied. Therefore the following boundary conditions were imposed

(7)

(8a)

(8b)

(9)

(10a)

(10b)

Equation 7 expresses the no-slip condition at the cylinder's surface _S. Equation 8a was used to impose that the velocity on the symmetry boundaries S₁, S₂, and S₃ has a zero normal component (n_i being the unit vector normal to the plane S_i). Equation 8b imposes that tangential stress is equal to zero on symmetry boundaries. Equation 9 expresses the periodicity of the velocity profile between outflow (₂) and inflow (₁) surfaces that are identical but rotated by 120° (R is the 120° rotation matrix around the z-axis).¹ Equations 10a and 10b represent the condition for normal stress in inflow and outflow boundaries. For any given value of the average pressure gradient P the problem of Eqs. 7-10 admits a unique solution, Eq. 14. The average velocity was then calculated as

(11)

where the integral of the vertical component of the velocity u_z is extended over the volume  available to water.

	METHOD OF SOLUTION

A Galerkin finite element solution of the boundary value problem represented by Eqs. 4 and 5 and 7-10b was obtained using the computational fluid dynamics package FIDAP (version 7.5; Fluid Dynamics International, Evanston, IL), running on an HP 9000/750 workstation. The typical finite element mesh used is shown in Fig. 3. The mesh was built in such a way to make the ₁ and ₂ surfaces identical, allowing imposition of the periodic boundary condition (Eq. 9) on corresponding nodes of ₁ and ₂. The number of nodes ranged from ~5,000 to ~7,000 according to different values of . We used 8-node isoparametric brick elements with linear basis functions for the velocity and with a discontinuous constant pressure approximation. A consistent penalty formulation was chosen, replacing the continuity Eq. 5 by

(12)

where we use the dimensionless variables

(13)

and V is the estimated average velocity. Typical dimensionless element sizes ranged from 2 × 10² to 9 × 10². We employed a value of the penalty parameter  = 10¹⁰ for all simulations. This choice was dictated by the requirement of a small compressibility error. To verify this condition, we performed a series of simulations in which several distributions of velocity at the inflow were imposed, leaving the velocity at the outflow free and relaxing the periodic boundary condition (Eq. 8). The flow rates across the inflow and the outflow were then numerically determined, and their relative difference taken as the compressibility error. We found that with a value of = 10¹⁰ the compressibility error was less than 0.5%. The discretized problem was solved using a successive substitution method. Once the solution was obtained, the integral in Eq. 11 was numerically evaluated using a built-in function of the software package.

View larger version (82K): [in this window] [in a new window]   Fig. 3.   Three-dimensional finite element mesh ( = 0.15).

To verify that the calculated K * did not depend appreciably on element type and on pressure approximation, we performed simulations with 27-node isoparametric brick elements, with quadratic basis functions for the velocity, and with three different pressure approximations (discontinuous trilinear pressure, discontinuous linear pressure on a local basis, and discontinuous linear pressure on a global basis). The differences calculated among the results of these simulations were less than 0.05%. We also tested the effect of varying the penalty parameter and noticed that changing  from 10¹² to 10⁹ did not affect the calculated K* appreciably.

	VALIDATION OF THE METHOD

To check the reliability of the described computational method, we performed a series of numerical simulations of the incompressible viscous flow perpendicular to the two-dimensional equilateral triangular array of parallel cylinders shown in Fig. 4 and compared the results with theoretical analysis available in the literature. This problem has been studied analytically and numerically by several authors with good agreement among them and with experimental results (see Ref. 21 for review). One of the most updated analysis is that by Sangani and Acrivos (27), who found a general series solution for the Stokes flow and determined the coefficients of the series by a numerical method. For small values of the fractional solid volume , their results for the dimensionless permeability (K* = K_D/r², where r is the cylinder radius) agree very well with an analytical expression found previously by the same authors (28) and subsequently improved by Drummond and Tahir (6). Edwards et al. (9) solved the Navier-Stokes problem in the same geometry by a Galerkin finite element method for several values of Re confirming, in the case Re = 0, the results of Sangani and Acrivos (27).

View larger version (26K): [in this window] [in a new window]   Fig. 4.   Two-dimensional equilateral triangular array of cylinders with the computational domain.

In the present work, we assumed the flow to be directed along the y-axis in the unit element delimited in Fig. 4 by the bold line. Since lines S₁ and S₂ are symmetry boundaries for velocity and pressure, and lines ₁, ₂, and ₃ are periodic boundaries for the velocity, we restricted the computational domain to one-fourth of the unit element, as shown in Fig. 4 (hatched region). We considered the boundary value problem represented by Eqs. 4, 5, and 7 and the following boundary conditions

(14a)

(14b)

(15)

(16a)

(16b)

which correspond, respectively, to the boundary conditions expressed by Eqs. 8-10b for the three-dimensional problem. As for the three-dimensional problem, the two conditions expressed by Eq. 15 are not independent.¹

A consistent penalty formulation was used with a penalty parameter  = 10⁸ (chosen to obtain a negligible compressibility error). We used four-node isoparametric quadrilateral elements with linear basis functions for the velocity and with a discontinuous constant pressure approximation. Typical dimensionless element size ranged from 2 × 10² to 9 × 10². The dimensionless permeability was obtained calculating the average velocity U under a given pressure gradient P. Variations of  from 10¹⁰ to 10⁶ produced changes in K * of less than 0.3% for all values of considered.

The results of our calculations are summarized in Table 1, where the values of K * derived from Ref. 27 are also reported for comparison. The latter values of K * were obtained from the dimensionless drag f = F/(µU ) (F being the drag per unit length on a cylinder) by means of the relation K * = /(f ), which follows from force balance on the unit element. As shown in Table 1, there is an excellent agreement between values of K * obtained by the present numerical analysis and those derived from Ref. 27. We also tested the dependence of the calculated K * values on the direction of the driving pressure gradient. To this purpose, we performed two simulations (for  = 0.4 and  = 0.1) of the flow directed along the x-axis of Fig. 4, using the same procedure. Calculated K* values differ from those determined for flow along the y-axis by less than 0.1%, indicating that, for Re = 0, K* is not importantly affected by the direction of the driving pressure gradient in respect to the cylinder arrangement.

	RESULTS AND DISCUSSION

Top Abstract Introduction Results & Discussion References

Numerical results for the dimensionless permeability K * of the three-dimensional array, for several values of , are reported in Table 2 and in Figs. 5 and 6. We obtained a satisfactory interpolation of these data using the following formula, with three freely adjustable parameters

(17)

Optimal values of the parameters, calculated by nonlinear fitting of the numerical data, were a = 0.05539, b = 2.382, and c = 0.7275. This fit is accurate within 3% throughout the interval 0.005    0.7 (see Fig. 5).

View larger version (11K): [in this window] [in a new window]   Fig. 5.   Dimensionless permeability K * as a function of  obtained by numerical analysis for the 3-dimensional array. Line represents the interpolation formula (Eq. 17).

View larger version (16K): [in this window] [in a new window]   Fig. 6.   Comparison of numerical results for the 3-dimensional array to theoretical formulas. See text for discussion.

At variance to the two-dimensional geometry, to our knowledge there are no other reports in the literature of direct determinations, either analytical or numerical, of the local flow fields in three-dimensional arrays of fibers. The available approaches to calculate K * are either semiempirical or based on the extrapolation to the three-dimensional geometry of solutions obtained for two-dimensional geometries. In the following, we compare our numerical results with the most widely used of these approaches.

According to the semiempirical equation of Carman-Kozeny (see Ref. 16 for review), the Darcy permeability of a porous medium is

(18)

where r_h is the hydraulic radius of the medium, defined as the ratio between the volume available to water and the wetted area, and G is the dimensionless Kozeny constant. Equation 18 has been tested experimentally for many porous media and gives accurate results only for values of  > 0.4-0.5. Some authors, however, claimed that in certain circumstances the equation can be extended to much lower values of  (3). For most porous media, the best agreement with experimental data is obtained with G  5, as reviewed in Ref. 16. To compare Eq. 18 with the results of our numerical simulations, we derived analytically the hydraulic radius of our three-dimensional periodic array, that is

(19)

We verified numerically that this expression is correct, calculating the volume available to water and the wetted area of the finite element meshwork. Combining Eqs. 19 and 2 with Eq. 18, and dividing by r², we obtain

(20)

Since the ratio r/L is uniquely determined by , K * in Eq. 20 depends only on  and G. To evaluate the agreement between Eq. 20 and our numerical analysis, we calculated the values of G that satisfy Eq. 20 for each entry of Table 2. As shown in Table 3, values of G are almost constant for  ranging from 0.4 to 0.7, indicating a good agreement with numerical results, but they increase importantly for lower values of . Thus, in line with previous observations obtained for other models of porous media (16), the Carman-Kozeny relation (with the constant G = 5.7) predicts quite accurately the Darcy permeability of an ordered array of fibers only for  > 0.4, whereas for lower values of  it sensibly overestimates K_D.

Some authors applied the Carman-Kozeny relation to the basement membrane (2, 3, 23, 26), assuming, at variance to our assumption, that the volume of fiber intersections is negligible. The hydraulic radius calculated with this assumption is r_h = 1/2r(1  )/; thus Eq. 18 gives the well-known expression

(21)

There are important deviations between the simplified and true value of r_h for our model (from 8 to 38% for  ranging from 0.1 to 0.7). We then calculated G values using Eq. 21 with values for  and K * reported in Table 2. As shown in Table 3, values of the constant G are not constant for the entire range of  values adopted, and G is lower than 4 for values of  > 0.4. We thus conclude that Eq. 21 does not fit well to our data.

Besides the Carman-Kozeny relation, other analytical approaches have been proposed to calculate K * for a three-dimensional fibrous matrix. Among them, we considered those of Iberall (19), Spielman and Goren (29), and Jackson and James (21) that are more often cited in the literature dealing with permeability of basement membrane. We applied the expressions proposed by these authors to our geometrical model and compared calculated K_D with our numerical results (see Fig. 6). At variance with the Carman-Kozeny relation, these approaches are expected to give good results only for low values of , since they are based on two simplifying assumptions that apply only to dilute fibrous media, as follows: 1) that the total drag exerted by a fluid is given by the sum of the drag on individual cylinders, neglecting their interactions; and 2) that the drag on a single cylinder is the linear combination of the drag exerted by the flow parallel and perpendicular to the fiber. On the basis of these assumptions, the dimensionless permeability K * for our model, with an arbitrary orientation of the average velocity U, is given by

(22)

where K*_t and K*_n are the dimensionless permeabilities of an array of parallel cylinders, with the same , with axis parallel and perpendicular to the mean flow direction, respectively.

Iberall (3, 19) proposed to calculate K*_t and K*_n by means of relations derived previously by Emersleben (11) and Lamb (22). In this way, assuming Re = 0, one obtains the so-called drag theory relation

(23)

As shown in Fig. 6, our numerical data agree with Eq. 23 only for very low values of  (<0.01) and differ importantly for higher values . The reason for this poor agreement is that the expression used for the drag due to perpendicular flow given in Ref. 22 is not accurate for low Re. In details, this drag tends to zero for Re  0, and only the contribution of the parallel fibers remains; that is, of course, unacceptable.

Spielman and Goren (29) determined the drag exerted on a cylinder for both parallel and perpendicular flow with an approach based on the Brinkmann equation and obtained implicit equations relating K*_t and K*_n to , which can be combined with Eq. 22 to give

(24)

where K₀ and K₁ are the modified Bessel functions of order zero and one. Data calculated using Eq. 24, as shown in Fig. 6, agree reasonably well with our numerical data for a wide interval of  values. In detail, for  < 0.5, the maximum deviation with numerical data is less than 18%.

The third analytical approach we have considered is that of Jackson and James (20), who determined K * for a cubic array of cylinders. As the fraction of perpendicular and parallel fibers in a cubic array and in our model is the same, these authors actually employed Eq. 22 to derive K*, and they used for K*_t and K*_n the dimensionless permeabilities of a two-dimensional square array of parallel cylinders calculated by Drummond and Tahir (6). We adapted this procedure to our tetrahedral array, considering K*_t and K*_n of a two-dimensional equilateral triangular array, also taken from Drummond and Tahir (6). Truncating the series in Ref. 6 to the second order in  we get²

(25)

(26)

Indeed our three-dimensional structure is more similar to the combination of four triangular arrays of parallel cylinders, with the axis oriented as the four cylinders in the unit element (see Fig. 1). Inserting Eqs. 25 and 26 into Eq. 22 and expanding in series for small , one obtains

(27)

Data obtained using Eq. 27 are shown in Fig. 6. The agreement with our numerical data is not satisfactory, as there are differences that vary from 20 to 50% as  increases. Considering, as in the original analysis of Jackson and James (21), only the first two terms of the series, these deviations are much higher. Furthermore, although the series coefficients appearing in the expressions of K*_t and K*_n given in Ref. 6 for equilateral triangular array are very close to those for square array, the coefficient 1.4975 in Eq. 27 differs sensibly from the coefficient 0.931 in the expression of K * for the cubic array reported by Jackson and James (21). To verify this discrepancy, we repeated their calculation and found that the value 0.931 reported in the original publication should be replaced by 1.476. Finally, it is interesting to note that Eqs. 24 and 27, although obtained by very different methods, have the same asymptotic expression for   0 

(28)

which apply well to our numerical data, whereas the drag theory relation (Eq. 23) has a different asymptotic expression.

On the basis of the above comparisons, we can conclude that none of the analytical expressions we considered predicts K * with very good accuracy in the whole range of  adopted in the present study. The best agreement is shown by Eq. 24 but is still unsatisfactory. Therefore, we propose to use our interpolating formula (Eq. 17) to express K * as a function of  for an ordered array of fibers such as that considered in our model. Whether Eq. 17 applies also to different spatial arrangements of fibers cannot be established on the basis of the present results and could eventually be decided with more complex sensitivity studies. In the case of basement membranes, since the actual ultrastructure is not known in detail, such studies would not be justified at the present moment. Our model is consistent with available limited observations of the GBM ultrastructure. Since the various analytical expressions considered do not well reproduce numerical data for such a model, we suggest to use Eq. 17, together with the definition of K * given in Eq. 3, to relate Darcy permeability of the GBM to the radius of the extracellular protein fibrils and the fractional solid volume.

As stated previously, experimental estimation of the dependence of K_D on r and  is technically difficult. There are only few reports on the dimensions of the protein fibers that compose the GBM. The available theoretical and experimental observations indicate that the radius of these fibers (r) should range from 0.8 to 1.0 nm (26). Much higher values have been obtained (r > 6 nm) using electron microscopy techniques (18, 24), but there are reasons to suspect that tissue processing for scanning electron microscopy may influence the dimensions of the fibrillar material observed within the GBM. Similarly direct measurements and theoretical predictions would indicate that fractional solid volume () of the GBM is about 0.1 (26). Using our interpolating formula, we calculated, assuming a value of  = 0.1, that K_D = 1.36 and 2.13 for r = 0.8 and 1.0 nm. These values are close to the values of the Darcy permeability of isolated rat GBM derived by Drumond and Deen (7), on the basis of the experimental results of Daniels et al. (5) and of Robinson and Walton (26), who reported values of K_D = 2.7 and 1.8 nm², respectively. Although there is great uncertainty on these parameters, it appears that our theoretical approach closely reflects independent estimations of the geometrical parameters of the GBM meshwork and its Darcy permeability.

The present analysis can be useful in estimating the dimensions of fibrous material that composes the GBM. Recent in vitro and in vivo studies (4, 5, 13, 26) reported estimates of K_D for the GBM, and there are evidences that changes in GBM hydraulic permeability develop in experimental and human glomerular injury (8, 10, 12, 13, 15, 32). Thus, it would be interesting to know whether these changes in permeability are related to changes in the spatial organization of the GBM fibrillar proteins or to dimensional changes of the fibrous or both. Since K_D depends on two parameters, however, estimates of K_D alone do not allow one to calculate both r and , but another independent relation is needed. A way to overcome this limitation could be the simultaneous consideration of the transport of water and macromolecules across the GBM. The permeability of a fibrous membrane to a molecule of given size depends on the geometrical parameters of the fibers. In detail, according to Ogston et al. (25), the permeability coefficients of fibrous membranes to neutral macromolecules can be expressed as a function of r and  (3). Combining this relationship with our expression of the Darcy permeability given by Eqs. 3 and 17, one can uniquely determine the values of  and r from direct measurements of hydraulic membrane permeability and filtration rates of neutral macromolecules.