A new model of filtration and macromolecules transport across capillary walls
AA NEW MODEL OF FILTRATION AND MACROMOLECULESTRANSPORT ACROSS CAPILLARY WALLS
LAURA FACCHINI, ALBERTO BELLIN, AND ELEUTERIO F. TORO
Abstract.
Metabolic substrates, such as oxygen and glucose, are rapidly de-livered to the cells of large organisms through filtration across microvesselswalls. Modelling this important process is complicated by the strong couplingbetween flow and transport equations, which are linked through the osmoticpressure induced by the colloidal plasma proteins. The microvessel wall isa composite media with the internal glycocalyx layer exerting a remarkablesieving effect on macromolecules, with respect to the external layer composedby the endothelial cells. The physiological structure of the microvessel is rep-resented as the superimposition of two membranes with different properties;the inner membrane represents the glycocalyx, while the outer membrane rep-resents the surrounding endothelial cells. Application of the mass conserva-tion principle and thermodynamic considerations lead to a model composedby two coupled second-order partial differential equations in the hydrostaticand osmotic pressures, one expressing volumetric mass conservation and theother, which is non-linear in the unknown osmotic pressure, expressing macro-molecules mass conservation. Despite the complexity of the system, the as-sumption that the properties of the layers are piece-wise constant allows us toobtain analytical solutions for the two pressures. This solution is in agreementwith experimental observations, which contrary to common belief, show thatflow reversal cannot occur in steady-state conditions unless the hydrostaticpressure in the lumen drops below physiologically plausible values. The ob-served variations of the volumetric flux and the solute mass flux in case of asignificant reduction of the hydrostatic pressure at the lumen are in qualita-tive agreement with observed variations during detailed experiments reportedin the literature. On the other hand, homogenising the microvessel wall intoa single-layer membrane with equivalent properties leads to a very differentdistribution of pressure across the microvessel walls, not consistent with ob-servations. Introduction
Transcapillary flow occurring in small and large capillaries plays a decisive rolein human physiology by ensuring an endless flow of oxygen and other electrolytesneeded to sustain cell metabolism. Altogether, the vessel wall operates as a semiper-meable membrane, which is selective with respect to the size of the molecules, suchthat water and electrolytes pass through the wall much more easily than the pro-teins. This leads to an ultrafiltrate, the interstitial fluid, with a substantially re-duced protein content. The interstitial fluid transfers oxygen and other nutrients
Mathematics Subject Classification.
Key words and phrases.
Ultrafiltration; Starling’s law; Capillary wall; Nonlinear transport ofmacromolecules.This work is partially funded by CARITRO (
Fondazione Cassa di Risparmio di Trento eRovereto , Italy). a r X i v : . [ q - b i o . S C ] A ug LAURA FACCHINI, ALBERTO BELLIN, AND ELEUTERIO F. TORO to the cell and receives carbon dioxide and other waste products before draininginto the lymphatic system and return to the bloodstream in the venous part of thecirculatory system.As evidenced by Starling (1896), volumetric flow through the microvessel wall iscontrolled by the net imbalance between the osmotic pressure of the plasma pro-teins and the capillary pressure generated by the heart beat. Both pressures canchange to exert a regulatory action on the filtration, such as for example duringexercise when an increased filtration is triggered by a larger capillary pressure andthe plasma volume reduces by up 20%. On the other hand, an increased filtrationoccurs during cardiac failure, which causes excess water accumulation in the tis-sues (oedema). Substantial movement of fluids occurs during the rapid swelling ofacutely inflamed tissues, while a rapid absorption of interstitial fluid into the bloodstream follows an acute hemorrhage.The microvessel wall is specialised to the function and the compartment of theorganism in which it operates. It is typically composed of a single layer of endothe-lial cells, which are internally coated with a − nm thick hydrated gel, calledglycocalyx. The glycocalyx protrudes into the lumen in hairy tufts, forming a size-and charge-selective molecular sieve to plasma proteins, while being permeable towater and small solutes (oxygen and other nutrients) (Levick, 2010, Ch. 9). Theendothelial cells are separated by inter-cellular clefts, which can be partially closedby junctional strands, thereby increasing the selectivity of the whole membrane.For example, in order to impede neurotoxic agents contained in the blood streamreaching the interstitial fluid, the clefts of cerebral capillaries are closed by multiplejunctional strands with no gaps. In addition, the external surface of the vessel iscoated by the pericytes, encircled by the basement membrane, in turn wrappedby astrocyte feet (Li et al, 2010), which introduces two additional layers, therebyreducing further permeability to macromolecules and forming the blood-brain bar-rier (see e.g., Levick, 2010). The breakdown of the blood-brain barrier with theassociated increase of vessel permeability has been observed in many brain dis-eases. Examples include stroke, traumatic head injury, Alzheimer’s disease, AIDS,brain cancer, meningitis etc. (see Li et al (2010)). In addition, blood-brain barrierrupture has been associated with multiple sclerosis, as discussed for instance byZamboni et al (2009), Singh and Zamboni (2009) and Haacke et al (2005).This two-layer structure of most microvessels has been evidenced by electronmicrograph after perfusion with cationized ferritin (Turner et al, 1983) and reflectsmorphometric measurements performed later (Hu et al, 2000; Adamson et al, 2004;Levick, 2010). Explicit modelling of the effect on glycocalyx and clefts at junctionsbetween the endothelial cells has been performed by solving the Navier-Stokes (NS)equations at the micro scale (see e.g., Sugihara-Seki and Fu, 2005; Sugihara-Sekiet al, 2008). The main drawback of this modelling approach, besides the highcomputational burden, is in the difficulty to model the interaction between macro-molecules and the fibre cells composing the glycocalyx, which feedbacks to thevolumetric flow through the osmotic pressure. To overcome this difficulty, hybridmethods have been used in which the glycocalyx has been modelled as a membrane(porous media), while the flow through the clefts has been modelled by solving theNS equations (Sugihara-Seki et al, 2008). A similar approach has been used byProsi et al (2005) and Formaggia et al (2009) to model mass transfer across arterialwalls in patients affected by atherosclerosis. MODEL OF PLASMA FILTRATION AND MACROMOLECULES TRANSPORT 3 (cid:0) (cid:0) yzx Figure 1.
Schematic of a capillary, whose wall is composed byfolded endothelial cells with the glycocalyx coated at their luminalside.A further simplified, yet effective, way to represent the vessel wall is by the su-perimposition of two membranes with different properties. The external membranemimics the effect of the monolayer of endothelial cells joined edge to edge along seg-ments forming an irregular pattern of connections, in a crazy-paving resemblance,without representing explicitly the structure of the clefts. The connections are par-tially closed by tight junctions, Figure 1. The internal membrane represents theglycocalyx coating the layer of endothelial cells (Levick, 2010). Considering thatthe single layer of endothelial cells is folded to form an annular semipermeable bar-rier around the blood stream, the transcapillary flow can be assumed as mainlyradial and orthogonal to the blood flow direction z (Figure 1).This simplified computational domain agrees with morphometric measurements(see e.g., Hu et al, 2000; Adamson et al, 2004; Levick, 2010), but differs from existingstudies in the way the main structural elements are combined.The classical way to model flow and transport across a microvessel is to representthe glycocalyx and the clefts as a homogeneous membrane, with equivalent prop-erties. The Starling’s law is then applied to this homogenised composite, such thatcapillary filtration rate can be written as proportional to the difference betweenthe hydrostatic and osmotic pressure drops between the blood and the interstitialfluid (Fu et al, 1994; Zhang et al, 2006). This simple conceptual model has beenshown to be unable to interpret the experiments conducted by Landis (1932) andsuccessively by Adamson et al (2004) and Hu et al (2000). In one of his experi-ments, Landis (1932) showed that at steady state, fluid exchange in perfused singlecapillaries of frog mesentery did not invert direction, leading to absorption, whenhydrostatic pressure inside the lumen was lowered below the limit value that theStarling’s law indicates for inversion (Levick, 2010). As a possible interpretationof the apparent breakdown of the Landis law, Michel (1997) and Weinbaum (1998)argued that the filtration rate may be controlled by the drop of osmotic pressurebetween the lumen and a position in the cleft at the contact with the glycocalyx, LAURA FACCHINI, ALBERTO BELLIN, AND ELEUTERIO F. TORO rather than the interstitium. This leads to an important change in the conceptualmodel, ruling out models with a single equivalent homogeneous membrane lumpingthe effect of both glycocalyx and the clefts at the junctions of the endothelial cells.To verify this conceptual model, and avoid complex micro-scale modelling in viewof applications at a larger scale, we propose to represent the vessel wall as the su-perimposition of two membranes with different properties. The internal membranerepresents the glycocalyx, while the external membrane is introduced to mimic theeffect of the endothelial cells.The specific objective of the present work is to solve analytically the coupledflow and transport equations, with the latter being non linear, resulting from theapplication of general physical and thermodynamic principles to the compositemembrane forming the vessel wall. The equations are written in radial coordinates,to take advantage of the radial symmetry of the vessels.The paper is organised as follows. The model we use to describe the physiologicalprocesses controlling filtration and macro-molecules transport across the vessel wallis described in Section 2. The analytical solution of the system of two differentialequations is reported in Section 2.6 in case of a discontinuous variation of physicalproperties between the two membranes. The case of smooth transition is exploredin Section 3 by using a suitable numerical scheme, while in Section 4 we comparethe results of our model with others from the literature. Finally, conclusions arereported in Section 5. 2.
The mathematical model
Statement of the problem.
We idealise the microvessel as two concentrichollow cylinders representing (from the lumen outward) the glycocalyx and thesurrounding endothelial cells. The two hollow cylinders are considered rigid, owingto the small compliance of microvessels, including venules (Levick, 2010). Theresulting computational domain is shown in Figure 2 with the dimensions of thetwo membranes reported in Table 1. Blood flow is along the longitudinal axis ofthe microvessel and we assume that the variation of the target macro-moleculeconcentration is small along the flow direction (Intaglietta et al, 1996).A widely accepted rheological model of blood flowing in vessel considers an inter-nal Red Blood Cells (RBCs) rich inner core surrounded by a relatively thin plasmalayer, which can be well approximated as a Newtonian fluid (Sriram et al, 2011).With the further assumption that the two cylindrical layers are homogeneous, flowacross the microvessel wall is radial and at a first approximation controlled by thelocal hydrostatic and osmotic pressures. In addition, we consider the case of a singlenot reacting molecule and isothermal conditions (Katchalsky and Curran, 1965).2.2.
Governing equations.
Under the above hypotheses, the solvent flow q v andthe diffusional macromolecule flow q d through the microvessel wall are coupled andgiven by the following phenomenological equations (Katchalsky and Curran, 1965): q v = − (cid:96) p ( ∇ p − σ ∇ Π) ,q d = σ(cid:96) p ∇ p − (cid:96) d ∇ Π , (1)where p is the hydrostatic pressure and Π is the osmotic pressure, which emergesbecause the size of the macro-molecule is comparable with the size of the aperturesin the glycocalyx and the endothelial cells. In addition, (cid:96) p = k/µ is the ratiobetween the hydraulic permeability of the membrane and the solvent viscosity, MODEL OF PLASMA FILTRATION AND MACROMOLECULES TRANSPORT 5 (cid:0) (cid:0) x rzθ r T r G r E Figure 2.
Sketch of the domain: a long hollow circular cylindercomposed by two homogeneous porous membranes representing theglycocalyx for r ∈ ( r G , r E ) and the endothelium for r ∈ ( r E , r T ) . σ ∈ [0 , is the reflection coefficient of the membrane and (cid:96) d is the diffusionalpermeability (Michel and Curry, 1999). Equations (1) are written for a singlemacromolecule. In case of two or more macromolecules the terms involving theosmotic pressure should be summed over all the relevant macromolecules. Theosmotic pressure depends on the solute (macromolecule) concentration c , throughthe following expression (Levick, 2010): Π =
R T c, (2)where R is the gas constant and T is the absolute temperature.The reflection coefficient σ in Equations (1) reflects the impediment exerted bythe pore to the free movement of the macromolecule and approaches zero as thecharacteristic size of the pore is much larger than the characteristic size of themacromolecule. In this situation, which is typical of small molecules, the effectof osmotic pressure tends to zero and the diffusion coefficient tends to the freediffusion coefficient, which depends only on the characteristics of the molecule andthe temperature.The total flux q s of macromolecule is given by the sum of the convective anddiffusive components: q s = c ( q v + q d ) . (3)Mass conservation of the flowing solvent and of the macromolecules under steady-state conditions leads to the following governing equations: (cid:26) div q v = 0 , div q s = 0 , (4) LAURA FACCHINI, ALBERTO BELLIN, AND ELEUTERIO F. TORO (cid:0) r G r E r T ( p G , Π G ) ( p T , Π T ) Figure 3.
Sketch of the domain indicating the relevant geometricelements: the internal radius r G at the lumen side of the microves-sel, the radius of the interface between the glycocalyx and theendothelial cells, r E , and the external radius r T . In addition, p G and Π G are the hydrostatic and osmotic pressures, respectively,within the lumen, while p T and Π T are the same quantities in theexternal interstitial space.which written in the radial coordinate system assume the following form: dd r (cid:18) r(cid:96) p d p d r (cid:19) − dd r (cid:18) r(cid:96) p σ dΠd r (cid:19) = 0 , dd r (cid:20) r(cid:96) p ( σ −
1) Π RT d p d r (cid:21) + dd r (cid:20) r ( (cid:96) p σ − (cid:96) d ) Π RT dΠd r (cid:21) = 0 , (5)which are defined in the interval r ∈ ( r G , r T ) from the lumen side of the glycocalyxto the external surface of the endothelial cells.In the present work we seek the analytical solution of this system of two non-linear coupled equations subjected to the following boundary conditions (Figure3): p G = p ( r G ) , p T = p ( r T ) , Π G = Π( r G ) , Π T = Π( r T ) , (6)where the subscripts G and T indicate the internal surface of the glycocalyx andthe external surface of the endothelial cells, respectively.2.3. Material properties.
For the analytical solution, we consider the propertiesof the membranes as piece-wise constant with a discontinuous (abrupt) change at r = r E , the interface between the glycocalyx and the endothelial cells. This abrupttransition is convenient for obtaining the analytical solution, but not necessarilyrepresents the real transition of the physical properties. As a possible alternativewe consider the following model of smooth transition (see Figure 3 for the meaningof the symbols): σ ( r ) = [1 − w ( r − r E )] σ G + [1 + w ( r − r E )] σ W ,(cid:96) p ( r ) = [1 − w ( r − r E )] (cid:96) Gp + [1 + w ( r − r E )] (cid:96) Wp ,(cid:96) d ( r ) = [1 − w ( r − r E )] (cid:96) Gd + [1 + w ( r − r E )] (cid:96) Wd , (7) MODEL OF PLASMA FILTRATION AND MACROMOLECULES TRANSPORT 7 where w = w ( r ) is the smoothing function defined as follows:(8) w ( r ) = r √ ε + r , where both sub- and super-scripts G and W indicate the properties of glycocalyxand the endothelial cells, respectively. With this function we can control how prop-erties vary at the interface between the two layers, with a discontinuous transitionoccurring for ε → . With ε > the transition becomes progressively smootherto simulate possible gradual transitions, with different degrees of smoothness asindicated by Sugihara-Seki and Fu (2005).2.4. Dimensionless flow and transport equations.
To facilitate the analysis,it is convenient to make the above steady-state flow and transport Equations (5)dimensionless with respect to the following reference quantities: the vessel wallthickness, ∆ r = r T − r G , for the length, the interstitial hydrostatic pressure, | p T | ,for both hydrostatic and osmotic pressures and finally (cid:96) Hp for both (cid:96) p and (cid:96) d , where (cid:96) Hp is the weighted harmonic mean for the hydraulic conductivity, of the two layerscomposing the microvessel wall: (cid:96) Hp = r T − r G r E − r G (cid:96) Gp + r T − r E (cid:96) Wp , (9)with (cid:96) Gp and (cid:96) Wp indicating the hydraulic conductivity of the glycocalyx and theendothelial cells layers, respectively. In addition, the dimensionless radius is definedas follows: r = ( r ∗ − r G ) / ∆ r = r ∗ / ∆ r − ξ , where r ∗ indicates the dimensionalradius varying from r ∗ = r G to r ∗ = r T and ξ = r G / ∆ r . With this definition thedimensionless radius r lies between 0 and 1.After these preliminary steps, system (5) assumes the following dimensionlessform (hereafter all the quantities are considered dimensionless, unless otherwisestated): dd r (cid:18) F d p d r (cid:19) + dd r (cid:18) G dΠd r (cid:19) = 0 , Π (cid:20) dd r (cid:18) H d p d r (cid:19) + dd r (cid:18) L dΠd r (cid:19)(cid:21) + dΠd r (cid:18) H d p d r + L dΠd r (cid:19) = 0 , (10)where the auxiliary functions F , G , H and L are defined as follows: F ( r ) = ( r + ξ ) (cid:96) p ( r ) , G ( r ) = − ( r + ξ ) (cid:96) p ( r ) σ ( r ) , H ( r ) = ( r + ξ ) (cid:96) p ( r ) [ σ ( r ) − , L ( r ) = ( r + ξ ) [ (cid:96) p ( r ) σ ( r ) − (cid:96) d ( r )] . (11)In addition, we consider the following boundary conditions, which are the di-mensionless counterpart of Equations (6): p ( r = 0) = p G , p ( r = 1) = p T , Π( r = 0) = Π G , Π( r = 1) = Π T . (12) LAURA FACCHINI, ALBERTO BELLIN, AND ELEUTERIO F. TORO
Discharge and flux reconstruction.
Solvent (volume) and solute fluxes aregiven by the following expressions: q v = − (cid:96) p (cid:18) d p d r − σ dΠd r (cid:19) ,q s = Π (cid:20) (cid:96) p ( σ −
1) d p d r + ( (cid:96) p σ − (cid:96) d ) dΠd r (cid:21) , (13)which are written in the dimensionless form with respect to the quantities intro-duced in Section 2.4. Finally, the total fluxes of solvent and solute crossing themicrovessel wall are given by:(14) J v = (cid:90) π q v ( r ) ( r + ξ ) dθ = − π ( r + ξ ) (cid:96) p (cid:18) d p d r − σ dΠd r (cid:19) and(15) J s = (cid:90) π q s ( r ) ( r + ξ ) dθ = 2 π ( r + ξ )Π (cid:20) (cid:96) p ( σ −
1) d p d r + ( (cid:96) p σ − (cid:96) d ) dΠd r (cid:21) , where θ is the angle as depicted in Figure 2, while the volume and solute fluxes aredimensionless with respect to (cid:96) Hp | p T | and (cid:96) Hp | p T | / ( RT ) , respectively. Notice thatin these two expressions all the quantities are dimensional.With the material properties defined through the dimensionless counterpart ofEquations (7) the solutions of the hydrostatic and osmotic pressures are differen-tiable everywhere, while transition between the material properties of the two layersis controlled by the parameter ε .2.6. Analytical solution of the single- and two-layer model.
The two differ-ential equations composing system (10) can be integrated with respect to r leadingto the following expressions for the fluxes: k = ( r + ξ ) (cid:96) p (cid:18) d p d r − σ dΠd r (cid:19) ,k = ( r + ξ )Π (cid:20) (cid:96) p ( σ −
1) d p d r + ( (cid:96) p σ − (cid:96) d ) dΠd r (cid:21) . (16)After some algebraic manipulations, system (16) can be written in the followingform: k = ( r + ξ ) (cid:96) p (cid:18) d p d r − σ dΠd r (cid:19) ,k = Π (cid:20) ( σ − k + ( r + ξ )( (cid:96) p σ − (cid:96) d ) dΠd r (cid:21) . (17)2.6.1. The single homogeneous layer solution.
The second equation of system (17)contains only the osmotic pressure Π as unknown and therefore can be solvedanalytically with the boundary conditions (12) for the case of a single equivalentlayer, and with additional conditions at the interface between the layers in case ofmultiple layers.We consider first the case of a single layer with equivalent membrane properties.In this case the solution of the system (17) is:(18) Π( r ) = k k σ − k ( r + ξ ) − ( k ) k ( σ − (cid:96) p σ − (cid:96) d , MODEL OF PLASMA FILTRATION AND MACROMOLECULES TRANSPORT 9 for the osmotic pressure Π and(19) p ( r ) = k (cid:96) p ln( r + ξ ) + k + σ Π( r ) == k (cid:96) p ln( r + ξ ) + k + k k σσ − k ( r + ξ ) − ( k ) k ( σ − (cid:96) p σ − (cid:96) d , for the hydrostatic pressure, when both k and k are non-zero and the materialproperties (cid:96) p , (cid:96) d and σ are equivalent parameters. In Equations (18) and (19) W( z ) is the Lambert W function (see e.g., Corless et al, 1996; Barry et al, 2000), alsocalled omega function or product logarithm, which is the solution of the followingalgebraic non-linear equation:(20) z = W( z ) e W( z ) . The solutions (18) and (19) require that four constants k , k , k and k be eval-uated by imposing boundary conditions on the hydrostatic and osmotic pressuresat the lumen and external surfaces of the microvessel. This leads to the followingexplicit expressions for k , k and k : k = (cid:96) p ( p G − p T ) − σ (Π G − Π T )ln( ξ ) − ln(1 + ξ ) , (21) k = ( p T − σ Π T ) ln( ξ ) − ( p G − σ Π G ) ln(1 + ξ )ln( ξ ) − ln(1 + ξ ) , (22) k = f e f ξ δ , (23)where f = ( k /k )( σ − G − and δ = k ( σ − / (cid:2) k ( (cid:96) p σ − (cid:96) d ) (cid:3) .By imposing that the osmotic pressure be equal to Π G at r = r G we obtain, aftera few manipulations, the following expression in the only unknown k , provided that k is given by Equation (21):(24) (cid:20) k k ( σ − T − (cid:21) e k k ( σ − T − − f e f (cid:18) ξ ξ (cid:19) ( k ) k ( σ − (cid:96) p σ − (cid:96) d = 0 . Equation (24) can be solved by using the Newton-Raphson method (Hildebrand,1987) with the following initial guess: k = Π G + Π T (cid:96) p ( σ − p G − p T ) + ( (cid:96) p σ − (cid:96) d )(Π G − Π T )ln( r G + ξ ) − ln( r T + ξ ) , (25)which is the exact solution of Equation (24) in the special case of k = 0 .2.6.2. The multiple layer solution.
The expressions (18) and (19) for the osmoticand hydrostatic pressures can be applied to all the layers of a multi-layer microves-sel, provided that the constants k and k , are layer-specific and that the properties (cid:96) p , (cid:96) d and σ are discontinuous across the boundary between adjacent layers. Onthe other hand, k and k are global quantities since they are equal to the volumeand solute mass fluxes (divided by ∓ π ) crossing the microvessel wall. Therefore,in addition to the 4 boundary conditions at the inner and outer surfaces, continuity of the hydrostatic and osmotic pressures as well as constancy of the volumetric andmacromolecule fluxes should be imposed at each interface. This results in a systemof n + 2 equations, in the same number of unknown consisting in the n valuesof both k and k , which are layer-specific quantities in addition to the two globalquantities k and k . In particular, for the two-layer model besides the bound-ary conditions (12) the following continuity conditions should be imposed at theinterface between the first and second layer at r = r E :(26) p ( r E ) = p ( r E ); Π ( r E ) = Π ( r E ) , where the subscripts “1” and “2” refers to the solution within the first (glycocalyx)and the second (endothelial cells) layer.The case of a smooth transition of these properties between the adjacent layerswill be discussed subsequently with the help of numerical solutions.By substituting the layer-specific quantities k i and k i reported in the Appen-dix A into the Equations (26) we obtain the following two equations in the twounknowns k and k : g e g − f G e f G (cid:18) ξr E + ξ (cid:19) δ G = 0 , (27) g e g − f W e f W (cid:18) ξr E + ξ (cid:19) δ W = 0 , (28)where g and g assume the following expressions g = k k σ G − σ G − σ W [( p T − σ W Π T ) − ( p G − σ G Π G ) + k β ] − , (29) g = σ W − σ G − g ) − , (30)with σ G (cid:54) = σ W and β = (cid:96) Wp ln( ξ ) + ( (cid:96) Gp − (cid:96) Wp ) ln( r E + ξ ) − (cid:96) Gp ln(1 + ξ ) (cid:96) Gp (cid:96) Wp . (31)In case the reflection coefficients are the same in both layers (i.e. σ G = σ W ),equations (27)-(28) become f G e f G (cid:18) ξr E + ξ (cid:19) δ G − f W e f W (cid:18) ξr E + ξ (cid:19) δ W = 0 , (32)with k = (cid:96) Gp (cid:96) Wp [( p G − p T ) − σ (Π G − Π T )] (cid:96) Wp ln( ξ ) + ( (cid:96) Gp − (cid:96) Wp ) ln( r E + ξ ) − (cid:96) Gp ln(1 + ξ ) , (33)where σ := σ G = σ W .For the two-layer model considered here it is more convenient to indicate with G quantities referring to the first layer (glycocalyx) and with W quantities referringto the second layer (the endothelial cells).Equations (27) and (28) can be solved by using the Newton-Raphson methodwith the initial guess for the unknowns k and k obtained by computing volumetricand solute mass fluxes for the case in which the interstitial pressures are appliedto the external surface of the glycocalyx at r = r E . These fluxes can be obtainedanalytically exactly for k and as a first-order approximation for k , as follows: MODEL OF PLASMA FILTRATION AND MACROMOLECULES TRANSPORT 11 k (0)1 = (cid:96) Gp ( p G − p T ) − σ G (Π G − Π T )ln( ξ ) − ln( r E + ξ ) , (34) k (0)2 = Π T k ( σ G − . (35)2.7. Parameters of two-layer model.
Table 1 shows typical values of the geo-metrical and physiological characteristics of an intact vessel, as well as the valuesof the osmotic and hydrostatic pressures within the lumen and in the externalinterstitial space utilised in the present study.The analytical solution for the two-layer case presented in Section 2.6 has beenobtained assuming a discontinuous transition of properties at the interface betweenthe glycocalyx and the endothelial cells. Smoother transitions are also possible andwill be analysed successively, by using a suitable numerical solution.Thermodynamic considerations on the phenomenological Equations (1) discussedin Katchalsky and Curran (1965), lead to the following constraint:(36) Pe < σ , where Pe = (cid:96) p /(cid:96) d is the Péclet number, which represents the reciprocal strength ofadvective and diffusive transport processes: when Pe is high advection dominatesover diffusion and vice-versa when Pe is small.Parameter [unit] Value Reference r G [ µ m] 5 Charm and Kurland (1974) r E [ µ m] 5.15 Adamson et al (2004) r T [ µ m] 5.5 Charm and Kurland (1974) Π G [mmHg] 25 Levick (1991) Π T [mmHg] 12 Levick (1991) p G [mmHg] 20 Levick (1991) p T [mmHg] − Levick (1991) α σ G σ W (cid:96) Gp [ µ m sec − mmHg − ] 0.601854 Speziale et al (2008) (cid:96) Wp [ µ m sec − mmHg − ] 4.15203 Speziale et al (2008) (cid:96) Gd [ µ m sec − mmHg − ] 0.536252 (cid:96) Wd [ µ m sec − mmHg − ] 3.69946 Table 1.
Typical values of the material properties of a microves-sel: σ is the reflection coefficient, (cid:96) p is the hydraulic conductivity, (cid:96) d is the diffusional permeability. The coefficient α depends on thePéclet number and is defined such as to respect condition (36). Thesuperscripts G and W indicate the glycocalyx and the endotheliallayers, respectively. In order to respect the condition (36) everywhere within the computational do-main, including the transition zone, we choose (cid:96) d as follows:(37) (cid:96) d ( r ) := α (cid:20) max r G
The analytical solutions described in the Section 2.6 are valid for a two-layermodel with discontinuous properties at the interface between the glycocalyx andthe endothelial cells. However, the case of smooth transition between the two lay-ers cannot be solved analytically, for which it is necessary to resort to numericalsolutions. In this section we describe numerical methods and we assess conver-gence properties, accuracy and efficiency for the case of a vessel composed of twohomogeneous layers with different material properties.3.1.
Description of numerical schemes.
Among the possible numerical schemes,which can be used to solve the Boundary Value Problem (BVP) (10)-(12), in thepresent work we consider a classical Finite Difference (FD) scheme and a Runge-Kutta shooting scheme. The domain [0 , is discretised by a regular mesh r i = ih ,for i = 0 . . . N + 1 , where h = 1 / ( N + 1) is the mesh spacing. The grid is designedin such a way that the interface point r E lies between two adjacent grid points.The unknowns are the functions p ( r ) and Π( r ) , for which we seek approximations p i ≈ p ( r i ) and Π i ≈ Π( r i ) . Following Freeze (1975) the diffusion terms of theEquations (10) are approximated as follows: dd r (cid:20) K ( r ) d f d r (cid:21) r = r i ≈ h (cid:18) K i + K i +1 f i +1 − f i h − K i − + K i f i − f i − h (cid:19) , (38)where K ( r ) is substituted with the functions F ( r ) , G ( r ) , H ( r ) or L ( r ) , in therespective diffusion terms in Equations (10). The application of this numericalscheme leads, after imposing the following boundary conditions: p = p G , p N +1 = p T , Π = Π G , Π N +1 = Π T , (39)to a sparse non-linear algebraic system of N equations in N unknowns, whichcan be solved by using the Newton method. The second strategy makes use ofthe shooting method, which is based on converting the BVP (10)-(12) to an initialvalue problem for the following augmented system in the unknowns y ( r ) = Π( r ) , y ( r ) = Π (cid:48) ( r ) , y ( r ) = p ( r ) , y ( r ) = p (cid:48) ( r ) :(40) y (cid:48) = y ,y (cid:48) = − FL − GH (cid:20) ( FH (cid:48) − F (cid:48) H ) y + ( FL (cid:48) − G (cid:48) H ) y + F ( H y + L y ) y y (cid:21) ,y (cid:48) = y ,y (cid:48) = 1 FL − GH (cid:20) ( GH (cid:48) − F (cid:48) L ) y + ( GL (cid:48) − G (cid:48) L ) y + G ( H y + L y ) y y (cid:21) , with the following initial conditions:(41) p (0) = p G , p (cid:48) (0) = dp G , Π(0) = Π G , Π (cid:48) (0) = d Π G . MODEL OF PLASMA FILTRATION AND MACROMOLECULES TRANSPORT 13
Initial value problem (40)-(41) is solved many times by using Runge-Kuttaschemes of orders from RK = 1 to RK = 4 , until convergence is achieved (Hilde-brand, 1987).The initial conditions are changed according to the boundary conditions at r = 1 and the procedure is stopped when the boundary values p N +1 and Π N +1 con-verge to p T and Π T , respectively. Given the initial slopes ( d Π ( k − G , dp ( k − G ) and ( d Π ( k ) G , dp ( k ) G ) , the initial conditions are updated according to the following linearinterpolation scheme: d Π ( k +1) G = d Π ( k ) G + d Π ( k ) G − d Π ( k − G Π ( k ) N +1 − Π ( k − N +1 (Π T − Π ( k ) N +1 ) ,dp ( k +1) G = dp ( k ) G + dp ( k ) G − dp ( k − G p ( k ) N +1 − p ( k − N +1 ( p T − p ( k ) N +1 ) , (42)until || Π T − Π ( k ) N +1 , p T − p ( k ) N +1 || is smaller than a given tolerance, where Π ( k ) N +1 and p ( k ) N +1 are the numerical solutions at the boundary r = 1 obtained by solving theinitial value problem with the initial conditions at the stage k .For the first two stages we used(43) (cid:40) d Π (0) G = Π T − Π G ,dp (0) G = p T − p G , (cid:40) d Π (1) G = α d Π (0) G ,dp (1) G = α dp (0) G , with α = 2 and α = 0 . .3.2. Assessment of the computational methods.
Preliminary simulations wereconducted by increasing the number of grid nodes by a factor of 2 at each refine-ment level k according to the following expression: N k = 2 k ( N −
1) + 1 , where N = 199 is the initial number of grid nodes. Empirical convergence is evaluatedthrough the following L -norm:(44) E k := || (Π , p ) k − ( ˆΠ , ˆ p ) || , which is normalised by the mesh spacing h k = 1 / (cid:2) k ( N − (cid:3) . Notice that the L -norm (44) is computed with reference to the vector of dimension N k , containingthe pairs (Π k , p k ) at each level k of discretisation. The reference solution ( ˆΠ , ˆ p ) isobtained numerically by using N = 50689 nodes after observing that the L -normof the difference between the numerical solutions obtained with the shooting RK method at the discretisation levels k = 8 and k = 7 was of the order of − .The convergence of the numerical solution with the number N k of grid nodesis shown in Figure 4. The rate of convergence increases with the order of thenumerical scheme, thereby confirming the applicability of the proposed methods tothe numerical solution of the non-linear system of differential Equations (10).For a sufficiently small error and a given mesh size the error reduces as the orderof the numerical scheme increases. Notice that RK and Freeze, being both ofsecond order of accuracy, show the same convergence rate (the slope of the curves),but the latter is affected by a smaller error. In addition, Figure 4 shows that thetarget error of − (see solid horizontal line) is attained with a mesh of N = 25345 points with the Freeze’s method, while for RK and RK shooting methods it isattained with coarser grids of N = 6337 and N = 1585 points, respectively. − − − − − Number of grid points N k E rr o r i np r e ss u r e E k FreezeShoot-RK1Shoot-RK2Shoot-RK3Shoot-RK4
Figure 4.
Error E k = || (Π , p ) k − ( ˆΠ , ˆ p ) || as function of grid pointsfor Freeze’s finite difference method and Runge-Kutta shootingmethods of orders 1 to 4.However, if a larger error is admitted, for example − , Freeze’s scheme reachesthe target at a coarser grid than the shooting methods up to the third order.The main indication provided by the Figure 4 is that from a computationalpoint of view it is more effective to improve accuracy of the solution by increasingthe order of accuracy of a scheme rather than refining the mesh. For relativelycoarse, yet acceptable error, the situation reverses and Freeze’s method becomesmore effective than the other schemes.4. Results
Validation of the numerical scheme.
In view of its applicability to the caseof smooth transition of the membrane properties at the interface between the gly-cocalyx and the endothelial cells, we compare the numerical solution obtained fromFreeze’s method with the analytical solution for discontinuous membrane propertiesdiscussed in Section 2.6. Preliminary simulations showed that a satisfying agree-ment between numerical and analytical solutions for the osmotic and hydrostaticpressures can be obtained with only 25 grid nodes. However, in order to obtain agood agreement also for the fluxes, the number of nodes should be increased sig-nificantly. Figure 5 shows the following relative differences: ∆ J i = | ( J N k i − J i ) /J i | , i = v, s , where J N k i is the volumetric flux (for i = v ) or the solute mass flux (for i = s ) computed numerically with N k grid nodes and J i is the corresponding flux ob-tained with the analytical solution, i.e. J v = − π k and J s = 2 π k , respectively.The relative differences of the two fluxes decline oscillating around at what appearsto be a common power law function of the number of grid nodes (notice the log-log scale used in the Figure 5). For N k = 18433 , the relative difference is smallerthan − and − for the volumetric and the solute mass fluxes, respectively.Therefore, in the following, if not explicitly stated, the numerical simulations areperformed with the Freeze’s scheme by using N k = 18433 grid nodes, which ensuresgood accuracy at a reasonable computational cost.Figure 6 compares the hydrostatic and osmotic pressures across the microvesselswall, obtained by solving numerically Equations (10) with the analytical solutions MODEL OF PLASMA FILTRATION AND MACROMOLECULES TRANSPORT 15 − − − − − − Number of grid points F l u x E rr o r Volume (cid:29)ux relative errorSolute (cid:29)ux relative error
Figure 5.
Relative differences between numerical and analyticalsolutions of the volumetric and solute mass fluxes, as a function ofthe number of grid nodes. . . . . . . . . . r Π ( r ) Analytical, − L Linearised, − L Analytical, − L Numerical, − L ( a ) . . . . . . . . . − r p ( r ) Analytical, − L Linearised, − L Analytical, − L Numerical, − L ( b ) Figure 6.
Comparison between the numerical solution of the os-motic (a) and the hydrostatic (b) pressures for the two-layer case,obtained with the Freeze’s scheme by using grid nodes, andthe corresponding analytical solutions (Analytical, 2-L). The singlelayer analytical solutions (Analytical, 1-L) are also shown togetherwith the linearised analytical solutions (Linearised, 1-L) presentedin the Appendix B. In all cases ε = 0 and for ease of represen-tation the numerical solution is shown only at a few grid points.The properties of the two layers are reported in the Table 1 to-gether with the Dirichlet boundary conditions at the lumen andinterstitial sides of the microvessel wall.(18) and (19), respectively. The difference between the analytical and the numericalsolutions is negligible with an error equal to .
12 10 − and .
15 10 − for the osmoticand hydrostatic pressures, respectively.4.2. Comparison between the two- and single-layer models.
In most ap-plications the microvessel is considered homogeneous, under the assumption thathomogeneous equivalent properties can be obtained that mimic the combined effectof the glycocalyx and the endothelial cells. Equivalent parameters can be defined as the parameters that when used into the solutions for a homogeneous media lead tothe same volumetric and solute fluxes of the heterogeneous (two-layer) media. How-ever, given the non-linearity of the governing Equations (10) equivalent parametersvalid for any choice of boundary conditions cannot be defined, since they depend ofthe structure of the governing equations and the boundary conditions as well (Mil-ton, 2002). Exploring this issue in depth would require a detailed analysis, which isbeyond the objectives of the present work, we therefore limit ourselves to computethe equivalent parameters for the boundary conditions and media properties of thebase case reported in the Table 1.The equivalent parameters to be used in the analytical solutions for a singlehomogeneous layer can be obtained by imposing that the fluxes are conserved, i.e.by imposing the following conditions:(45) k H = k , and k H = k , where the superscript H indicates that the flux is evaluated with the single-layermodel, while k and k are the fluxes of the heterogeneous two-layer model. All thefluxes, that in the Equations (45) are divided by ∓ π , are obtained as describedin Section 2.6 with the parameters and the boundary conditions showed in Table1. Since the equivalent parameters to be defined are three ( (cid:96) eqp , (cid:96) eqd , σ eq ), whilethe conditions imposed by Equations (45) are two, we set the equivalent reflectioncoefficient by using the following expression suggested by Sugihara-Seki and Fu(2005):(46) σ eq = (cid:96) Gd (cid:96) Wd (cid:96) Gd + (cid:96) Wd (cid:18) σ G (cid:96) Gd + σ W (cid:96) Wd (cid:19) , which for the media properties shown in Table 1 assumes the following value: σ eq =0 . . With this value of σ eq the two Equations (45) can be solved obtaining (cid:96) eqp =0 . , (cid:96) eqd = 0 . , which substituted into Equations (18) and (19) provide thebehaviour of the osmotic and hydrostatic pressures, respectively, for the equivalenthomogeneous single-layer media.The analytical solution of the osmotic pressures is shown in Figure 6a. Theosmotic pressure declines rapidly across the glycocalyx, reaches a minimum at theinterface with the endothelial cells and then it increases again to the value imposedas boundary condition at the external surface of the microvessel. This behaviour,and in particular the minimum of the osmotic pressure at the external surface ofthe glycocalyx, is in agreement with a recent reinterpretation of the Starling lawproposed independently by Michel (1997) and Weinbaum (1998), which providesan improved interpretation of the classic experiments conducted by Landis (1927);see also Levick and Michel (2010) for a complete review. Dilution in the clefts justoutside the glycocalyx is an important physiological mechanism, which has beenindicated by Michel (1997) and Weinbaum (1998) as the cause preventing reversalsteady-state flow (absorption) when capillary hydrostatic pressure was lowered to − cmH O ( . − . mmHg ) in the Landis experiment. A similar behaviouris shown in Figure 6b for the hydrostatic pressure with a strong reduction across theglycocalyx followed by a mild reduction across the endothelial cells. Hydrostaticpressure is not differentiable at the interface between the two layers, which is dueto the discontinuity in the media properties, but the pressure gradient does notreverse across the endothelial cells, as for the osmotic pressure. The importantresult shown in the Figures 6a and 6b is that most of the pressures drop between MODEL OF PLASMA FILTRATION AND MACROMOLECULES TRANSPORT 17 . . . . . . . . . r Π ( r ) Numerical, ε = 10 − Numerical, ε = 10 − Analytical, ε = 0 ( a ) . . . . . . . . . − r p ( r ) Numerical, ε = 10 − Numerical, ε = 10 − Analytical, ε = 0 ( b ) Figure 7.
Osmotic (a) and hydrostatic (b) pressures across themicrovessel wall for ε = 0 , − , − . The case ε = 0 is obtainedby means of the analytical solutions, while the cases with ε =10 − , − are obtained numerically with the Freeze’s scheme byusing 18433 grid nodes.the lumen and the interstitium occurs in the glycocalyx, confirming the importanceof this hydrated gel in controlling flow and solute mass exchange (see e.g., Levick,2010).A striking difference can be observed in the Figures 6a and 6b between the single-and two-layer models, with the latter showing a smooth but steep decline of both thehydrostatic and osmotic pressures within the glycocalyx. This is due to the strongsieving effect that glycocalyx exerts on macromolecules, such that only a very smallfraction of them reached the clefts. In the two-layer model this effect is reproducedby using a large reflection coefficient. Notice that, due to the larger aperture ofthe clefts, macromolecules move with small to negligible impediment, as soon asthey have crossed the glycocalyx. In the membrane model adopted in this work thealmost free movement of the macromolecules in the clefts is represented by adoptinga small reflection coefficient ( σ = 0 . ). Because of the high selectivity of theglycocalyx, concentration of macromolecules is small at the interface between theglycocalyx and the endothelial cells, resulting in an osmotic pressure smaller thanin the interstitium. This feeds back to the hydrostatic pressure, which also showsa strong decline within the glycocalyx, as discussed above. This behaviour, whichis consistent with the observation that flow cannot be reversed by simply reducingthe hydrostatic pressure in the lumen as discussed by Levick and Mortimer (1999)and Levick and Michel (2010), is not captured by the single-layer model, whichinstead predicts much higher pressures at the interface between the glycocalyx andthe endothelial cells and a gradual decline of the pressures across the microvesselwall, with a gradient that increases with the distance to account for the progressiveincrease of the surface crossed by the flows.An important consequence of this different behaviour of the pressures is that thesingle-layer model is unable to capture the effect on volumetric and solute massfluxes of glycocalyx deterioration, which being located in the lumen side of themicrovessel is more prone to be damaged, than the endothelial cells.Figures 7a and 7b show the distribution of the osmotic and hydrostatic pres-sures, respectively, across the microvessel wall for the following three values of theparameter controlling property variations at the interface: ε = 0 , − , − . The solution for the discontinuous transition, i.e. for ε = 0 is the analytical solutiondiscussed in Section 2.6, while for ε > the solutions are numerical and obtainedwith the Freeze’s scheme by using grid nodes. A smooth, yet sharp, tran-sition in the material properties eliminates the discontinuity in the first derivativeof the pressures at the interface between the two layers and contemporaneouslypressures within the glycocalyx are steeper and with a smaller curvature than inthe discontinuous case. The pressures at the position where the interface is locatedin the discontinuous case show negligible variations such as the distribution of thepressures within the endothelial cells. The progressive increase of the gradientsof the hydrostatic and osmotic pressures within the glycocalyx occurring when ε increases, which is accompanied by the reduction of the reflection coefficient closeto the interface, leads to an increase of both the volume and solute mass fluxes(see Table 2). The relative increase of J v is negligible ( . ) for ε = 10 − , butit increases rapidly with ε , reaching for ε = 10 − . J s is more sensitive tovariations of ε , with an increase of . and . for ε = 10 − and − , re-spectively, with respect to the solute mass flux obtained with a sharp transition( ε = 0 ) of the material properties.Description J v J s Analytical ε = 0 545 .
586 2802 . Numerical ε = 0 545 .
607 2802 . Numerical ε = 10 − .
354 3110 . Numerical ε = 10 − .
809 3945 . Table 2.
Volumetric and solute mass fluxes for the following val-ues of ε , the parameter controlling the smoothness of the proper-ties transition between the two layers: ε = 0 , − , − .Figures 8a and 8b show the effect of changes in the osmotic pressure at the lu-men side (of the blood) on the behaviour of the osmotic and hydrostatic pressures,respectively. A variation of the lumen osmotic pressure, with respect to the ref-erence value of Π G = 25 mmHg shown in Table 1, causes a variation of the samesign, but smaller, in the osmotic pressure at the interface between the two layers.An opposite behaviour is observed for the hydrostatic pressure at the interface,which as shown in the Figure 8b reduces as the lumen osmotic pressure increases.Interestingly, the change in the osmotic pressure drop across the microvessel wallfeeds back to the hydrostatic pressure distribution through the coupling with thenon-linear transport equation.This effect cannot be reproduced with linearised models decoupling flow andtransport processes, such as that presented by Speziale et al (2008).The impact of osmotic pressure variations in the lumen is shown in the Figure9. The most relevant information contained in the figure is the opposite behaviourof the two fluxes; an increase of the lumen osmotic pressure with respect to thereference case, with all the other quantities remaining the same, leads to a reductionof the volumetric flux and a contemporaneous increase in the solute mass flux. Theopposite occurs, when the Π G is reduced below the reference case: the volumetricflux increases, while the solute mass flux reduces, as an effect of the reduction inthe osmotic pressure drop across the microvessel wall. MODEL OF PLASMA FILTRATION AND MACROMOLECULES TRANSPORT 19 . . . . . . . . . r Π ( r ) Π G = 28Π G = 24Π G = 20Π G = 16Π G = 12Π G = 8Π G = 4 ( a ) . . . . . . . . . − r p ( r ) Π G = 28Π G = 24Π G = 20Π G = 16Π G = 12Π G = 8Π G = 4 ( b ) Figure 8.
Behaviour of the osmotic and hydrostatic pressuresacross the microvessel wall for several values of the osmotic pres-sure Π G in the lumen. , , , ,
600 Π G J v , , , , , , J s Figure 9.
Volumetric and solute mass fluxes, relative to the ref-erence case with the boundary conditions and material propertiesshown in Table 1, as a function of the blood (lumen) osmotic pres-sure.Figures 10a and 10b show the behaviour of the osmotic and hydrostatic pressuresacross the microvessel wall for several values of the blood hydrostatic pressure p G within the lumen. The reduction of the osmotic pressure at the interface betweenthe two layers, with respect to the interstitial value, becomes progressively smalleras the hydrostatic pressure within the lumen reduces, and it vanishes at p G (cid:39) .For smaller values of p G the osmotic pressure at the interface between the two layersremains higher than the external osmotic pressure. At the interface the hydrostaticpressure is higher for higher lumen hydrostatic pressures, but to a lesser extentwith respect to the increase in the lumen. This leads to a higher pressure drop forhigher lumen hydrostatic pressures. A similar behaviour is shown by the osmoticpressure, but with a smaller variations in the pressure drop, due to the fact thatthe osmotic pressure at the lumen does not change. As shown in Figure 11, bothfluxes increase with the hydrostatic pressure at the lumen. However, the volumetricflux reduces to zero as the lumen hydrostatic pressure reduces to mmHg . Thisresults is consistent with the Landis (1932) experiments showing no volumetric flux . . . . . . . . . r Π ( r ) p G = 40 p G = 35 p G = 30 p G = 25 p G = 20 p G = 15 p G = 10 ( a ) . . . . . . . . . r p ( r ) p G = 40 p G = 35 p G = 30 p G = 25 p G = 20 p G = 15 p G = 10 ( b ) Figure 10.
Behaviour of the osmotic (a) and hydrostatic (b) pres-sures across the microvessel wall for several values of the hydro-static pressure p G within the lumen.
10 15 20 25 30 35 40 − , , , p G J v
10 15 20 25 30 35 401 , , , , , , J s Figure 11.
Volumetric and solute mass fluxes, with the boundaryconditions and material properties shown in Table 1, as a functionof the blood (lumen) hydrostatic pressure.inversion at steady state at pressures as low as cm H O ( . mm Hg ), which hadbeen considered the basis for revisiting the Starling’s law (Michel, 1997; Weinbaum,1998). 5. Conclusions
We have presented and discussed a new model of flow and transport of macro-molecules (proteins) across the composite wall of a microvessel. The microvessel isrepresented as a two-layer hollow cylinder. The inner layer represents the glycoca-lyx, an hydrated membrane exerting a remarkable sieving effect on macromolecules,and the external layer representing the endothelial cells, which are folded and con-nected along clefts spiralling in an irregular manner along the longitudinal microves-sel axes. The clefts are partially closed by the tight junctions. We represent thiscomposite media as two membranes of different thickness and properties. Flow andnon-linear transport equations are coupled through the osmotic pressure, which isassumed proportional to the concentration of macromolecules in the plasma. Weshow that, by assuming radial symmetry, this model can be solved analytically for
MODEL OF PLASMA FILTRATION AND MACROMOLECULES TRANSPORT 21 the general case of n − layers. The solution is consistent with the mechanistic revis-itation of the classical Starling law proposed independently by Michel (1997) andWeinbaum (1998). In particular, it well represents the dilution occurring in the cleftspace at the external surface of the glycocalyx, with the corresponding reduction ofthe osmotic pressure to values smaller than in the external tissues, which is in linewith recent observations (Adamson et al, 2004) and claimed as the main mecha-nism preventing flow inversion at low hydrostatic pressures. Our model differs fromother published models in several aspects. Differently from Speziale et al (2008) wesolve the full system of coupled differential equations for flow and transport withoutlinearising the transport equation in a n − layer setup, which allows us to handlespecialised microvessels. For simplicity, the application was limited to a two-layermicrovessel, which is the most common type of microvessel in humans and othermammals. However, the extension to four layers, typical of brain microvessels, canbe obtained at the cost of a more complicated structure of the solution, due tothe need to impose conservation of hydrostatic and osmotic pressures across thethree interfaces, while conservation of volumetric and mass fluxes are obtained byimposing that the coefficients k and k are the same in the three layers. The appli-cation of the model to an homogenised microvessel, representing the combined effectof glycocalyx and endothelial cells with a single layer membrane characterised bysomewhat equivalent properties, as suggested by Speziale et al (2008) for example,evidenced a strikingly different distribution of the pressures within the microvesselwall, which are significantly higher than those of the two-layer model, in partic-ular at the interface between the two layers. A better match may be obtained ifboundary conditions are applied to the external surface of the glycocalyx, therebyneglecting the effect of the endothelial cells and the dilution occurring in the cleftat the contact with the external surface of the glycocalyx, which has been indicatedas an important physiological mechanism controlling the volumetric flux (see e.g.,Levick and Michel, 2010).To summarise, our solution of the n − layer model of microvessel has a level ofcomplexity comparable to existing homogenised single-layer models (Speziale et al,2008), but showed to be much more accurate in describing the combined effectof glycocalyx and endothelial cells, including the dilution occurring in the cleft atthe contact with the external surface of the glycocalyx, on controlling volumetricflow and solute mass transport across the microvessel wall. Our model is compu-tationally much more effective than micro-scale approaches, such as that proposedby Sugihara-Seki et al (2008); with a moderate effort it can be implemented intolarge-scale models representing blood circulation in the human body. This is animportant feature that full micro-scale models, resorting to sophisticated numericalmethods and requiring parallel computing, cannot enjoy. In addition, the betterreproduction of the hydrostatic pressure across the microvessel wall, with respectto the homogeneous single layer model, makes this approach appealing for appli-cations dealing with the mechanical response of the microvessel to changes of theinternal hydrostatic pressure. Appendix A. Parameters for a two-layer microvessel
The imposition of the pressures at the inner surface of the glycocalyx equal tothe pressures in the blood stream and of the pressures at the outer surface of theendothelial cells equal to the pressures in the interstitial leads to the following expressions for k i and k i : k G = ( p G − σ G Π G ) − k (cid:96) Gp ln( ξ ) ,k W = ( p T − σ W Π T ) − k (cid:96) Wp ln(1 + ξ ) ,k G = f G e f G ξ δ G ,k W = f W e f W (1 + ξ ) δ W , (47)where f G = k k ( σ G − G − ,f W = k k ( σ W − T − ,δ G = ( k ) k ( σ G − (cid:96) Gp σ G − (cid:96) Gd ,δ W = ( k ) k ( σ W − (cid:96) Wp σ W − (cid:96) Wd . (48) Appendix B. Approximate analytical solutions for a single-layermicrovessel
The first equation of system (10) can be written as a function of the net pressure P ( r ) = p ( r ) − σ Π( r ) as follows: ( r + ξ ) (cid:96) p d P d r = c , (49)then the linearised version of the second equation of system (17), obtained bydecomposing the osmotic pressures in a mean value Π m plus a perturbation (cid:15) ( r ) and neglecting the terms of the second order and higher in the perturbation, assumethe following form: ( σ − c d (cid:15) d r + ( (cid:96) p σ − (cid:96) d ) Π m (cid:20) d (cid:15) d r + ( r + ξ ) d (cid:15) d r (cid:21) = 0 , (50)in the unknown function (cid:15) .Equations (49)-(50) with boundary conditions P (0) = P G = p G − σ Π G ,P (1) = P T = p T − σ Π T ,(cid:15) (0) = Π G − Π m = Π G − Π T ,(cid:15) (1) = Π T − Π m = − Π G − Π T , (51)can be separately solved analytically to obtain: P ( r ) = c ln( r + ξ ) + c , Π( r ) = c ( r + ξ ) − β + c ,p ( r ) = P ( r ) + σ Π( r ) = c ln( r + ξ ) + σc ( r + ξ ) − β + c + σc , (52) MODEL OF PLASMA FILTRATION AND MACROMOLECULES TRANSPORT 23 where the coefficients assume the following expressions: c = ( p G − p T ) − σ (Π G − Π T )ln( ξ ) − ln(1 + ξ ) ,c = Π G − Π T ξ − β − (1 + ξ ) − β ,c = [ p T ln( ξ ) − p G ln(1 + ξ )] − σ [Π T ln( ξ ) − Π G ln(1 + ξ )]ln( ξ ) − ln(1 + ξ ) ,c = Π G ξ β − Π T (1 + ξ ) β ξ β − (1 + ξ ) β , (53)for(54) β = 2( σ − (cid:96) p σ − (cid:96) d ( p G − p T ) − σ (Π G − Π T )[ln( ξ ) − ln(1 + ξ )](Π G + Π T ) (cid:54) = 0 . References
Adamson R, Lenz J, Zhang X, Adamson G, Weinbaum S, Curry F (2004) Oncoticpressures opposing filtration across non-fenestrated rat microvessels. Journal ofPhysiology 557(3):889–907Barry D, Parlange J, Li L, Prommer H, Cunningham C, Stagnitti F (2000) Ana-lytical approximations for real values of the Lambert W-function. Mathematicsand Computers in Simulation 53:95–103Charm S, Kurland G (1974) Blood Flow and Microcirculation. Wiley, New YorkCorless R, Gonnet G, Hare D, Jeffrey D, Knuth D (1996) On the Lambert Wfunction. Advances in Computational Mathematics 5:329–359Formaggia L, Quarteroni A, Veneziani A (2009) Cardiovascular Mathematics.Springer, MilanFreeze R (1975) A stochastic-conceptual analysis of one-dimensional groundwaterflow in nonuniform homogeneous media. Water Resources Research 11(5):725–741Fu B, Weinbaum S, Tsay R, Curry F (1994) A junction-orifice-fiber entrance layermodel for capillary permeability: Application to frog mesenteric capillaries. Jour-nal of Biomechanical Engineering 116:502–513Haacke E, Cheng N, House M, Liu Q, Neelavalli J, Ogg R, Khan A, Ayaz M, KirschW, Obenaus A (2005) Imaging iron stores in the brain using magnetic resonanceimaging. Magnetic Resonance Imaging 23(1):1–25Hildebrand F (1987) Introduction to Numerical Analysis. Dover PublicationsHu X, Weinbaum S (1999) A new view of Starling’s hypothesis at the microstruc-tural level. Microvascular Research 58:281–304Hu X, Adamson R, Liu B, Curry F, Weinbaum S (2000) Starling forces that opposefiltration after tissue oncotic pressure is increased. American Journal of Physiol-ogy – Heart and Circulatory Physiology 279(4):H1724–H1736Intaglietta M, Johnson P, Winslow R (1996) Microvascular and tissue oxygen dis-tribution. Cardiovascular Research 32:632–643Katchalsky A, Curran P (1965) Non-equilibrium Thermodynamics in Biophysics.Harvard Univ. PressLandis E (1927) Micro-injection studies of capillary permeability. II. The relationbetween capillary pressure and the rate at which fluid passes through the wallsof single capillaries. American Journal of Physiology 82(2):217–238
Landis E (1932) Factors controlling the movement of fluid through the humancapillary wall. Yale Journal of Biology and Medicine 5(3):201–225Levick J (1991) Capillary filtration-absorption balance reconsidered in light of dy-namic extravascular factors. Experimental Physiology 76:825–857Levick J (2010) An Introduction To Cardiovascular Physiology. Hodder ArnoldLevick J, Michel C (2010) Microvascular fluid exchange and the revised Starlingprinciple. Cardiovascular Research 87(2):198–210Levick J, Mortimer P (1999) Fluid ‘balance’ between microcirculation and inter-stitium in skin and other tissues: Revision of the classical filtration-reabsorptionscheme. In: Messmer K (ed) Microcirculation in chronic venous insufficiency.Progress in Applied Microcirculation, vol 23, Karger S., pp 42–62Li G, Yuan W, Fu B (2010) A model for the blood-brain barrier permeability towater and small solutes. Journal of Biomechanics 43(11):2133–2140Michel C (1997) Starling: The formulation of his hypothesis of microvascular fluidexchange and its significance after 100 years. Experimental Physiology 82(1):1–30Michel C, Curry F (1999) Microvascular permeability. Physiological reviews79(3):703–761Michel C, Phillips M (1987) Steady-state fluid filtration at different capillary pres-sures in perfused frog mesenteric capillaries. Journal of Physiology 388:421–435Milton G (2002) Theories of Composites. Cambridge University Press, CambridgeProsi M, Zunino P, Perktold K, Quarteroni A (2005) Mathematical and numericalmodels for transfer of low-density lipoproteins through the arterial walls: A newmethodology for the model set up with applications to the study of disturbedlumenal flow. Journal of Biomechanics 38:903–917Singh A, Zamboni P (2009) Anomalous venous blood flow and iron deposition inmultiple sclerosis. Journal of Cerebral Blood Flow and Metabolism 29(12):1867–1878Speziale S, Tenti G, Sivaloganathan S (2008) A poroelastic model of transcapillaryflow in normal tissue. Microvascular Research 75(2):285–295Sriram K, Vazquez B, Yalcin O, Johnson P, Intaglietta M, Tartakovsky D (2011)The effect of small changes in hematocrit on nitric oxide transport in arterioles.Antioxidants and Redox Signaling 14:175–185Starling E (1896) On the absorption of fluids from the connective tissue spaces.Journal of Physiology 19(4):312–326Sugihara-Seki M, Fu B (2005) Blood flow and permeability in microvessels. FluidDynamics Research 37:82–132Sugihara-Seki M, Akinaga T, Itano T (2008) Flow across microvessel walls throughthe endothelial surface glycocalyx and the interendothelial cleft. Journal of FluidMechanics 601:229–252Turner M, Clough G, Michel C (1983) The effects of cationised ferritin and nativeferritin upon the filtration coefficient of single frog capillaries. Evidence that pro-teins in the endothelial cell coat influence permeability. Microvascular Research25(2):205–222Weinbaum S (1998) 1997 Whitaker distinguished lecture: Models to solve mysteriesin biomechanics at the cellular level; A new view of fiber matrix layers. Annalsof Biomedical Engineering 26(4):627–643Zamboni P, Galeotti R, Menegatti E, Malagoni A, Tacconi G, Dall’Ara S, Bar-tolomei I, Salvi F (2009) Chronic cerebrospinal venous insufficiency in patients
MODEL OF PLASMA FILTRATION AND MACROMOLECULES TRANSPORT 25 with multiple sclerosis. Journal of Neurology, Neurosurgery, and Psychiatry80:392–399Zhang X, Adamson R, Curry F, Weinbaum S (2006) A 1-d model to explore theeffects of tissue loading and tissue concentration gradients in the revised Starlingprinciple. American Journal of Physiology - Heart and Circulatory Physiology291(6):H2950–H2964
Department of Mathematics, University of Trento (Italy), via Sommarive 14,38123 Trento (TN)
E-mail address : [email protected] Department of Civil, Environmental and Mechanical Engineering, University ofTrento (Italy), via Mesiano 77, 38123 Trento (TN)
E-mail address : [email protected] Laboratory of Applied Mathematics. DICAM, University of Trento (Italy), viaMesiano 77, 38123 Trento (TN)
E-mail address ::