Hrev_master Veins and Lymphatics 2017; volume 6:6849 [Veins and Lymphatics 2017; 6:6849] [page 101] A time-dependent multi-layered mathematical model of filtration and solute exchange, the revised Starling principle and the Landis experiments Laura Facchini,1 Alberto Bellin,2 Eleuterio F. Toro3 1Department of Mathematics, University of Trento; 2Department of Civil, Environmental and Mechanical Engineering, University of Trento; 3Laboratory of Applied Mathematics. Department of Civil, Environmental and Mechanical Engineering, University of Trento, Italy Abstract Cell oxygenation and nutrition is vitally important for human and animal life. Oxygen and nutrients are transported by the blood stream and cross microvessel walls to penetrate the cell’s mem- brane. Pathological alterations in the transport of oxygen, and other nutrition elements, across microvessel walls may have seri- ous consequences to cell life, possibly leading to localized cell necrosis. We present a transient model of plasma filtration and solute transport across microvessel walls by coupling flow and transport equations, the latter being non-linear in solute concentra- tion. The microvessel wall is modeled through the superimposition of two or more membranes with different physical properties, rep- resenting key structural elements. With this model, the combined effect of the endothelial cells, the glycocalyx and other coating membranes specific of certain microvessels, can be analyzed. We investigate the role of transient external pressures in the study of trans-vascular filtration and solute exchange during the drop of blood capillary pressure due to the pathological decrease of blood volume called hypovolaemia, as well as hemorrhage. We discuss the advantage of using a multi-layered model, rather than a model considering the microvessel wall as a single and homogeneous membrane. Introduction All cells of the human body need constant oxygenation and sustenance, provided by blood circulation and an efficient microvascular exchange. In order to preserve volume homoeosta- sis, the plasma volume that globally filtrates through microvessels is removed from the interstitium. The traditional view, by now dis- proved, is that small veins and venules continuously re-absorb interstitial fluid. This seems no longer true, because, under steady- state conditions, capillary volumetric flow leads to a rise of the peri-capillary interstitial osmotic (oncotic) pressure, which impedes flow reversal.1 The modern view is that reabsorption is operated by the lymphatic system (see e.g. Levick, chapter 112). The volumetric flux per unit area (JV/A) of plasma through cap- illary wall is controlled by the Starling equation:3 JV/A = Lp[(pc–po)– σ(Πc–Πo)], (1) where p and Π indicate fluid pressure and osmotic pressure, respectively, while the subscript indicates where the pressures are considered, with c referring to the lumen of the capillary and o referring to the interstitium, outside the microvessel. Furthermore, Lp = ℓp ∆r is the mechanical filtration capacity, given by the product of the permeability ℓp and the thickness ∆r of the membrane. Finally, σ ≤ 1 is the reflection coefficient (see e.g. Katchalsky and Curran, chapter 10.34), which should be introduced to take into account the permeability of the microvessel wall to plasma pro- teins. According to equation (1), a reduction of pc might lead to reab- sorption of interstitial fluid, i.e. a negative JV, under the hypothesis that both Πo and po remain constant. This observation leads to the traditional interpretation that volume homoeostasis in the tissues is ensured by reabsorption of microvascular filtration by low-pres- sure capillaries and venules. This view has been recently corrected by accurate interpretations of experimental studies, including those of Landis.5,6 In these experiments, a sudden drop of pc resulted in a temporary absorption, which however is rapidly counterbalanced by the rise of Πo due to the increase of plasma proteins concentra- tion in the interstitium, which re-establishes a positive flux.1,7-9 Sustained absorption is possible only in specialized tissues, such as the intestinal mucosa, the renal cortex and lymph nodes, where the interstitium receives an external source of fluid compensating for the reabsorbed fluid, such that Πo remains constant, despite the reabsorption. The reasoning behind the demonstration that absorption cannot be maintained in low-pressure microvessels, is based on a static interpretation of Starling’s equation (1) through the construction of the absorption curve, i.e. a curve showing the evolution of the absorption pressure Πc −Πo against the filtration rate JV. In particu- lar, this curve is compared with Starling’s equation (1), which is depicted in the same graph as a straight line (see e.g. Levick, Figure 71). This analysis has the merit of evidencing the main mechanisms controlling filtration across microvessels, but cannot be used to actually model the transient behavior resulting from per- turbation of one of the four terms in equation (1). Before discussing existing models of microvascular filtration and formulate our objectives, we will describe the structure of Correspondence: Laura Facchini, Laboratory of Applied Mathematics, Department of Civil, Environmental and Mechanical Engineering, University of Trento, Italy. E-mail: lari317@yahoo.it Key words: Revised Starling principle; Landis experiments; filtration and solute transport; time-dependent multi-layer one-dimensional model. Funding: this work is partially funded by Fondazione Cassa di Risparmio di Trento e Rovereto (CARITRO, Italy), project no. 2011.0214. Contributions: LF, AB and EFT conceived the work and derived the model; LF performed the simulations; LF and AB wrote the paper. Conflict of interest: the authors declare no potential conflict of interest. Received for publication: 12 June 2017. Revision received: 12 September 2017. Accepted for publication: 13 September 2017. This work is licensed under a Creative Commons Attribution 4.0 License (by-nc 4.0). ©Copyright L. Facchini et al., 2017 Licensee PAGEPress, Italy Veins and Lymphatics 2017; 6:6849 doi:10.4081/vl.2017.6849 No n c om me rci al us e o nly Article [page 102] [Veins and Lymphatics 2017; 6:6849] microvessels, because it is crucial to the construction of the con- ceptual model. The extra-cranial capillary vessel wall is known to be com- posed by a one-cell-thick layer of endothelial cells, internally coat- ed by the surface glycocalyx, and currently recognized as crucial for microvascular wall homoeostasis. Measurements show that glycocalyx thickness varies between 50 and 500 nm,10 though they may underestimate the real thickness due to the dehydrating effect of electron microscopy fixation and processing.11 In fact, more recent measurements using cryo-EM12 show that glycocalyx may be as much as 10 to 20-times thicker than previously thought. The fiber matrix of the glycocalyx was recently described by Squire et al.13 as a quasi-periodic three-dimensional fibrous mesh work with a characteristic spacing of 20 nm and anchoring foci (forming a hexagonal array), emanating from the underlying cortical cytoskeleton and extending from the luminal side of the endothe- lial cell surface into outer regions of the intercellular clefts.10 The presence of the glycocalyx and its sieving effect was first observed by Michel7 and Weinbaum,8 who argued that volumetric flow and solute transport across microvessel wall are driven by drops in pressures between the lumen and the endothelial cleft, just outside the glycocalyx, instead of between the lumen and the inter- stitium, as commonly assumed before. The endothelial cells are separated by clefts, which are seen as long slits partially sealed by dynamic complexes of cell-cell junc- tional proteins. They provide anchorage and cell stability, as in the case of the adherent junctions, or form an almost impermeable bar- rier to plasma and macromolecules (tight/occludens junctions), in particular in the continuous capillaries composing the blood-brain barrier.2,14 A number of mathematical models have been developed to simulate trans-vascular transport, chiefly across the arterial wall in the study of atherosclerosis, such as the model proposed by Pertktold et al.,15 based on early works by Quarteroni et al.16 and Prosi et al.,17 and recently extended to retinal microcirculation by Causin et al.18 This model couples volumetric flux across the ves- sel wall with hemodynamics under the hypotheses that the flux across the microvessel wall is steady and that the osmotic pressure can be neglected such that the volumetric flux is given by the Darcy’s law. In addition, solute flux is modeled through the linear advection-diffusion-reaction equation (decoupled from the flow equation), while the Navier-Stokes equations are used to model hemodynamics in the lumen. The glycocalyx layer is lumped into an interface, together with the endothelium, and considered as a membrane through the non-linear algebraic Kedem-Katchalsky equations.19 This study highlighted difficulties in coupling Darcy’s with the Navier-Stokes equations.15 In the works by Khakpour and Vafai20-22 and Ai and Vafai,23 flow and transport equations are decoupled. A review of some attempts to couple flow and transport equations can be found in Khakpour and Vafai.24 Several one- and three-dimensional inter- pretations have been proposed since 1950s,10,25 showing that per- meability to water and hydrophilic solutes in microvessels is con- trolled by the geometry of inter-endothelial cleft pathways and junction strands, leading to discontinuous leakages, and by a sur- face fiber layer (the glycocalyx), coating internally the endothelial cells. In the present paper, we propose a fully coupled non-linear model of flow and transport of macromolecules across microvessel walls, under the hypothesis that the vessel is circular and that flow and transport are radial and transient. The microvessel wall is rep- resented as the superimposition of two membranes, one internal and the other external, representing the glycocalyx and the endothelium, respectively. The mathematical model Statement of the problem The microvessel wall is represented as two concentric circular hollow cylinders, both with homogeneous but contrasting proper- ties, representing the glycocalyx and the endothelial cells, respec- tively. Owing to the relatively small variation of both hydrostatic and osmotic pressures in the direction of vessel’s axis, flow and trans- port are assumed radial in the cross-section of the vessel. Furthermore, deformations due to membrane compressibility are assumed to be small enough to neglect their effect on the geometry of the computational domain. Blood is a non-Newtonian fluid, but since red blood cells concentrate at the center of the vessel only plasma flows across the microvessel wall. Plasma is well repre- sented by a Newtonian fluid with characteristics similar to those of water.10 One of the most reliable theories describing the molecular fil- ter of the trans-vascular pathway, confirmed by quantitative evi- dence of ferritin exclusion, is the fiber matrix theory.10 This theory describes the trans-vascular molecular filter through the assump- tion of glycocalyx extension from the luminal side of the endothe- lial cell surface into outer regions of intercellular clefts.10 This jus- tifies our assumption of smooth, though sharp, transition between the properties of the two membranes. Governing equations in dimensionless form Under the above assumptions, the volumetric flow (specific discharge) qV = qV(x,y,z,t) and the total solute flow qS = qS(x,y,z,t) are related to the hydrostatic p = p(x,y,z,t) and the osmotic Π = Π(x,y,z,t) pressures through the following phenomenological expressions, derived from the application of the Onsager law (see Katchalsky and Curran, chapters 8 and 104 and Facchini et al.26): (2) where ℓp is the hydraulic conductivity, σ is the membrane reflection coefficient and ℓd is the diffusional permeability. Osmotic pressure and macromolecules concentration c = c(x,y,z,t) are related through the van’t Hoff Law:2 Π = R T c, (3) where R is the gas constant and T is the absolute temperature. Mass conservation of blood plasma and of the solute (i.e., macromolecules such as plasma proteins), applied under the hypothesis that the solid matrix is at rest and that the solution is diluted such that the density of the mixture differ slightly from the density of the blood plasma, leads to the following two coupled governing equations:4 (4) where the volumetric and solute fluxes appearing on the left-hand side of equation (4) assume the expressions given in the system of No n c om me rci al us e o nly Article equations (2) and the left-hand terms take into account the varia- tions in time of the fluid volume and solute stored within the mem- brane. In equations (4), ρ is the fluid (blood plasma) density, g is the standard gravitational acceleration and SS is the specific stor- age, which measures the storage capacity of the vessel walls, and is given by the following expression: SS = ρg (βm + nβw), (5) where βm and n are blood vessel wall compressibility and porosity, respectively, and βw is fluid compressibility. The hypothesis here is that during the transient state changes of porosity are small such that changes in the stored blood plasma occur with the membrane total volume remaining almost constant. Consequently, according to the classical theory of porous media,27 the combined effect of membrane and blood plasma compressibility can be represented through the specific storage SS given by the equation (5). After these preparatory steps, system (4) can be rewritten as follows (after switching to cylindrical coordinates): (6) for r ∈ (0,1) and t > 0, where α = (βm + nβw) pR, ξ = rc/(ro − rc), A = ℓp (σ − 1) and B = ℓpσ − ℓd. Hereafter, all the quantities are dimen- sionless, if not explicitly stated differently. All pressures are made dimensionless with respect to the reference pressure pR given by the (dimensional) initial external hydrostatic pressure |p0o|, the dis- tances with respect to the thickness ro − rc of the vessel shifted by the internal radius of the vessel wall identified with the internal surface of the glycocalyx rc, and time with respect to reference time given by tR=(∆r)2/(ℓpHpR). Furthermore, the reference for the physiological parameters is the harmonic mean ℓpH of hydraulic conductivity, weighted by vessel wall thickness: (7) where ℓpG refers to glycocalyx i.e. for r ∈ (rc,rg) and ℓpW to the endothelial cells composing the vessel wall r ∈ (rg,ro). The volumetric JV and solute JS fluxes crossing the microvessel wall are obtained by integrating the respective flows along a circle of radius r as follows: (8) (9) where volumetric and solute fluxes are dimensionless with respect to ℓpH |po0| and ℓpH |po0|2/(RT), respectively. Physiological parameters The reflection coefficient is let to vary according to the follow- ing expression:26 (10) with the usual meaning of the subscripts. Similar expressions are used for ℓp and ℓd. The transition between the two layers is smooth for ε > 0, and becomes discontinuous for ε → 0. In the present work, we assume ε2 = 10–4, unless otherwise stated, which produces a smooth, yet sharp, transition in line with experimental observa- tions.10 We consider passive transport of albumin, which accounts for half the plasma protein mass and generates about two-thirds of the colloid osmotic pressure.2 Table 1 shows typical values of geometrical and physiological properties of an intact capillary, as well as osmotic and hydrostatic pressures within the lumen and in the external interstitial space, at initial (with superscript 0) and final state of the transient experi- ment (see the Reproduction of the Landis experiment section). Also the compressibility of the fluid and of the membrane composing the microvessel wall, as well as media porosity, are reported. Reproduction of the Landis experiment During its experiments,5,6 Landis found a linear relation between the rates of filtration/absorption JV and capillary pressure pc in a population of frog mesentery capillaries. The traditional interpretation of the Landis5 experiment assumes that both hydro- static and osmotic pressures in the interstitium remain constant during reabsorption induced by a sudden drop of hydrostatic pres- sure below the osmotic pressure in the lumen. As discussed in the Introduction, this leads to the wrong interpretation that at the lower end of capillary beds and in venules absorption might be perma- nent. However, more accurate interpretations by Levick1 and Michel7 showed that reabsorption can only occur temporarily. The main mechanism preventing reabsorption in steady state is the increase of protein concentration in the interstitium due to flow reversal. This causes Πo to increase, thereby opposing flow reversal and finally leading to a steady-state condition with a smaller than before, but still positive (i.e., from the lumen to the interstitium) volumetric flow.2,9 This mechanism has been demonstrated by Levick1 by using a quasi-static interpretation (as discussed in the Introduction) and here we show that this is implied in the structure of the model (6). Initial and boundary conditions To numerically solve the PDE system (6), proper initial condi- tions should be assigned to both hydrostatic p(r,0) = p0(r) and osmotic Π(r,0) = Π0(r) pressures. The initial conditions for p0(r) and Π0(r) are provided according to the analytical solutions of the steady-state counterpart of the model (6), obtained by Facchini et al.26 and reported in the Steady-state implicit exact solutions sec- tion. The boundary conditions of hydrostatic and osmotic pressures used in the steady-state solutions are reported in Table 1 with sub- scripts c and o referring to the lumen of the capillary and the inter- stitial space, respectively. The thickness of the two membranes representing the glycocalyx and the endothelial cells, as well as the respective phenomenological parameters (i.e., ℓp, ℓd and σ), are also [Veins and Lymphatics 2017; 6:6849] [page 103] No n c om me rci al us e o nly Article reported in Table 1. Furthermore, the maximum simulation time is set to T = 3 min. At = T/5 the hydrostatic pressure, pc0, at the lumen is suddenly reduced from 22 mmHg to 7 mmHg, while the osmotic pressure is kept constant and equal to Πc0 = 24 mmHg (Figure 1): (11) where H(x) is the Heaviside step function36 defined as: (12) In the following we assume two different conditions for the interstitial pressures. First, we assume that both hydrostatic and osmotic pressures remain constant and equal to po0 and Πo0, respec- tively, during the time interval t = [0, T] spanning the entire simu- lation period, as shown in the last two rows of Figure 1A: (13) This condition mimics the traditional interpretation, which implies persistent reabsorption if the hydrostatic pressure in the lumen (pc) drops below po + σ (Πc − Πo). In fact, according to equa- tion (1), this is the value of pc below which volumetric flow JV reverses becoming negative. Then we consider the following boundary conditions for the interstitial pressures (Figure 1B): Figure 1. Behavior of the boundary conditions for (A) constant external pressures as described in equations (11)-(13), and (B) pressures varying with time according to equations (11)-(14). [page 104] [Veins and Lymphatics 2017; 6:6849] A B Table 1. Typical values of the material properties of a capillary: σ is the reflection coefficient, ℓp is the hydraulic conductivity, ℓd is the diffusional permeability. The superscripts G and W indicate the glycocalyx and the endothelial layers, respectively. As described in Initial and boundary conditions section, the bound- ary conditions change following equations (11)-(13)-(14) with initial (at t = 0) and final (at t = T) values of pressures during the transient experiment are denoted by superscripts 0 and T, respec- tively. The coefficient βw refers to the compressibility of horse blood plasma, while βm is the maximum value of blood vessel wall compressibility measured in dog descending thoracic aorta. The porosity n was measured in arterial graft prosthesis. Parameter [unit] Value Reference rc [µm] 5 [28] rg [µm] 5.15 [29] ro [µm] 5.5 [28] Πc0 [mmHg] 25 [9] Πo0 [mmHg] 8 [9] pc0 [mmHg] 22 [9] po0 [mmHg] −1.3 [9] ΠoT [mmHg] 15 [9] pcT [mmHg] 7.36 [9] poT [mmHg] −3.6 [9] σG 0.9 [30] σW 0.1 [31] ℓpG [µm2sec−1mmHg−1] 0.6013 [32] ℓpW [µm2sec−1mmHg−1] 4.1520 [32] ℓdG [µm2sec−1mmHg−1] 0.5357 [26] ℓdW [µm2sec−1mmHg−1] 3.6995 [26] βw [mmHg−1] 5.45 · 10−8 [33] βm [mmHg−1] 9.21 · 10−6 [34] n 0.5 [35] No n c om me rci al us e o nly Article (14) with (15) where W(x) is the Lambert W function,37,38 also called omega func- tion or product logarithm, which is the inverse function of z → z ez. In addition, τ is the time at which pressures in the interstitium reach steady state. In the present example, we assumed τ = 0.6 T, such as to mimic the recorded behavior of the pressures in the interstitium during the Landis5 experiment, as reported in Fig. 3 of Levick and Michel.9 Results of the transient simulations We first consider a vessel wall composed by two nearly homo- geneous layers, with physiological parameters varying according to equation (10) with ε2 = 10−4 (same for ℓp and ℓd), assuming either constant or time-varying pressures at the interstitium, according to equations (13) or equations (14), respectively. The physiological parameters are therefore nearly constant inside the two layers and vary sharply, yet continuously at the interface between the two lay- ers. Table 2 shows volumetric and solute fluxes at initial (t = 0) steady-state condition, immediately before and after the drop of the luminal hydrostatic pressure at t=t̂, and at the successive time t = T. Volumetric and solute flows reverse, i.e. both become nega- tive, immediately after the drop of the luminal pressure (Figure 2). Successively, volumetric flux remains negative if the interstitium pressures are kept constant according to the scheme shown in Figure 1A, but it rebounds and becomes positive at successive times for the more realistic transient boundary conditions shown in Figure 1B. These boundary conditions represent better than the case of constant interstitium pressures the physiological response Table 2. Volumetric and solute fluxes at selected times obtained by solving PDE system (18). The numerical solutions for the pressures were computed using a spatial grid of 1000 points and Dt = 500·(Dx)2. The external pressures are assumed constant as in equations (13) and transient according to equations (14). The values of the fluxes are reported both for the two- and single-layer model with equivalent dimensionless parameters equal to σeq = 0.798802, ℓpeq = 0.831226 and ℓdeq = 0.578381, computed by following the procedure described by Facchini et al.26 Layers BCs t = 0 t = t̂ – t = t̂ + t = T Volumetric flux [10−3 µm/sec] 2 Constant 27.4867 27.4867 −25.6272 −14.5961 2 Transient 27.4867 27.4867 −25.6267 3.1622 1 Transient 27.4867 27.4867 −12.1853 4.7463 Solute flux [10−2 mol/µm/sec] 2 Constant 14.6849 14.6849 −8.52912 6.3217 2 Transient 14.6849 14.6849 −8.53067 8.08544 1 Transient 14.6849 14.6849 −1.95824 4.9169 [Veins and Lymphatics 2017; 6:6849] [page 105] Figure 2. Behavior of the volumetric (A) and solute (B) fluxes obtained by numerically solving the system of differential equations (22) subjected to the initial and boundary conditions described in Initial and boundary conditions section and reported in Figure 1. In par- ticular, the external pressures are assumed constant as in equations (13) (dashed curves) and varying with time following equations (14) (solid curves). The numerical solutions for the pressures are computed using a spatial grid of 1000 points and Dt = 500 · (Dx)2. A BNo n c om me rci al us e o nly Article to flow reversal, which causes a temporary reduction of plasma volume within the interstitium, accompanied by a rise of the osmotic pressure, due to the increase of solute concentration, thereby contrasting flow reversal. Comparison between the single- and two-layer models In our recent work, we showed that in order to properly model transport of proteins across microvessel walls and the loss of the barrier effect consequent to glycocalyx damage, a two-layer model should be used, with the internal layer representing the glycocalyx and the external one the endothelial cells.26,39 However, in applica- tions, a single-layer model is typically used, requiring the defini- tion of equivalent properties encapsulating the effects of both gly- cocalyx and of the surrounding endothelial cells. To properly address the single-layer case we suggest computing the reflection coefficient by means of the expression proposed by Sugihara-Seki and Fu:10 (16) while ℓpeq and ℓdeq are obtained by imposing that at the initial steady-state condition volumetric and solute fluxes are the same in the two-layer model, considered here with discontinuous (step) transition of the physical properties, and in the single-layer homo- geneous equivalent model: (17) where (pc,Πc) and (po,Πo) are the initial values of the pressures at r = 0 and within the interstitium at r = 1, respectively. Figure 3 shows the comparison between volumetric and solute fluxes obtained with the single- and the two-layer models. Before the drop of the pressure at t=t̂, the two fluxes are the same in the two-layer and single-layer equivalent model as imposed to com- pute ℓpeq and ℓdeq. At larger times, the single-layer model produces higher volumetric, but lower solute fluxes during the transitory and the following new steady-state condition (Table 2). At the beginning of the transient phase following the drop of the pressure pc, solute flux reverses and assumes the following values: JS = −8.53067 and −1.95824 for the two- and single-layer model, respectively (Table 2). In agreement with what we already observed for the stationary case,26 which here has been assumed as initial con- dition, most of the excess of pressure is dissipated within the glyco- calyx. Due to the sudden drop in the lumen, the hydrostatic pressure reduces sharply, though smoothly, assuming values smaller than po (the value in the interstitium) in the last portion of the glycocalyx and across the endothelial cells, where it slightly rises while approaching the interstitium (Figure 4A). This causes a temporary reversal volumetric flow, which induces an increase of solute con- centration and therefore of osmotic pressure in the outer portion of the glycocalyx and in the clefts (Figure 4B). In fact, due to the rela- tively high velocity of the plasma in the glycocalyx and the solute slow motion in the endothelial clefts, which is the consequence of the small gradient of the osmotic pressure in this region, plasma pro- teins are expelled to the interstitium while flow prevents them from going back by diffusion, thus creating a zone at the contact between the glycocalyx and the endothelial cells in which their concentration (and thus the osmotic pressure) is lower than in the surroundings (Figure 4B). This dynamic, due to the coupling between flow and transport in a composite membrane, contrasts the reversal of the vol- umetric flow, which after a relatively short transient behavior becomes again positive, though much smaller than in the initial steady-state conditions (Figure 3A). This behavior is in line with the physiological explanation provided by Levick1 and discussed in the Introduction. The single-layer model does not reproduce correctly the behavior of the hydrostatic and osmotic pressures across the gly- cocalyx and the endothelial cells, though the transient regime of both Figure 3. Comparison of volumetric (A) and solute (B) fluxes obtained with the single- and two-layer models with boundary conditions given by equations (14). The numerical solutions for the pressures are computed using a spatial grid of 1000 points, Dt = 500 · (Dx)2 and physiological parameters σ(r), ℓp(r) and ℓd(r) varying according to equation (10) for σ and a similar expression for the other two parameters. For the single-layer case, the equivalent dimensionless physiological parameters are: σeq = 0.798802, ℓpeq = 0.831226 and ℓdeq = 0.578381. [page 106] [Veins and Lymphatics 2017; 6:6849] A BNo n c om me rci al us e o nly Article the volumetric and solute fluxes is somewhat captured. In particular, the single-layer model is unable to reproduce the dilution effect at the contact between the glycocalyx and the endothelial cells (com- pare Figures 4 and 5), and the solute flux is underestimated during the transient behavior and the following new steady-state condition (Figure 3B). Concluding remarks We have presented and discussed a two-layer model of tran- sient flow and transport of solutes across microvessel walls. The dynamic model allows computing the evolution in time of volu- metric and solute fluxes resulting from perturbations of state vari- ables, such as the hydrostatic and osmotic pressures. We have applied the model to represent the transient behavior occurring after a sudden drop in the luminal pressure, as done in the famous experiment by Landis.5 Our model is able to dynami- cally reproduce the evolution in time of pressures and fluxes and confirms the emergence of a transient state of reabsorption as observed by Landis5 and successively justified by phenomenolog- ical analyses.1,7-9 This transient state ends with the development of a new steady-state condition in which reabsorption disappears since all fluxes are positive, thereby ruling out the interpretation that reabsorption of interstitial fluid occurs at the level of low-pres- sure capillaries and venules, in agreement with the most recent physiological studies.9 In the present work, we have attempted to reproduce the values of the fluxes in Landis experiments at our best, due to the lack of exact values for the physiological parame- ters of the microvessels involved in the experiment. Nonetheless, we obtained an analogous behavior and plausible values of pres- sures and fluxes. One advantage of our model is that it allows sim- ulating the consequence of glycocalyx damage in all pathologies in which this is deemed important. Figure 4. Hydrostatic (A) and osmotic (B) pressures in the two-layer model as a function of distance from the internal glycocalyx surface for external pressures varying according to equations (14), and for several times. PDE system (18) has been numerically solved with a spatial grid of 1000 points and Dt = 500 · (Dx)2 and physiological parameters σ(r), ℓp(r) and ℓd(r) varying according to equation (10) for σ and a similar expression for the other two parameters. Figure 5. Hydrostatic (A) and osmotic (B) pressures in the single-layer model as a function of distance from the internal glycocalyx sur- face for external pressures varying according to equations (14), and for several times. PDE system (18) has been numerically solved with a spatial grid of 1000 points and Dt = 500 · (Dx)2. The equivalent dimensionless physiological parameters are σeq = 0.798802, ℓpeq= 0.831226 and ℓdeq = 0.578381. [Veins and Lymphatics 2017; 6:6849] [page 107] A B A BNo n c om me rci al us e o nly Article [page 108] [Veins and Lymphatics 2017; 6:6849] In a broader context, the present model may be extended and used in conjunction with global models for the dynamics of the human circulation.40-42 In particular, it has been hypothesized that extracranial venous strictures43 that hamper brain venous return may induce cerebral venous hypertension.42,44,45 It would be of interest, for example, to study the effect of cerebral venous hyper- tension on brain perfusion. Numerical solution of the governing equations Numerical scheme To simplify the notation of the numerical scheme, we rewrite the non-linear system of PDEs (6) as follows: (18) for r ∈ (0,1) and t > 0, in the unknown p and Π, both functions of r and t, where the auxiliary functions are defined as: (19) To numerically solve system (18), we apply the implicit Crank- Nicolson finite-difference scheme proposed by Freeze,46 for a con- stant spacing ∆x, which leads to the following nonlinear system of algebraic equations: (20) with the following source terms: (21) Notice that the governing equations are written in a radial ref- erence system, such that ∆x is the spacing along the radial coordi- nate. The algebraic non-linear system (20) is solved by using Newton method. In equations (20) and (21), ∆x is the mesh spacing of the spatial grid {ri}i=0..N+1 and ∆t denotes the time step between two adjacent temporal grid points {tn}n=0..M+1, while hydrostatic and osmotic pressures are approximated by pin ≈ p(ri,tn) and Πin ≈ Π(ri,tn), respectively. Furthermore, In order to solve system (6), suitable initial and boundary con- ditions, the latter at all times t ∈ [0,T], should be defined. In addi- tion, the grid is designed such that rg falls between two nodes, of a grid composed of 1000 nodes, unless otherwise stated; the toler- ance for stopping the iterations of the Newton method is set to 10−10, unless otherwise stated. Validation of the numerical scheme: convergence to steady state As discussed in the Numerical scheme section, the coupled non-linear transient flow and transport model (6) is solved by using the numerical scheme proposed by Freeze46 and shown by Facchini et al.26 to compare well with the analytical solution of the stationary counterpart of model (6). Therefore, we first analyze the convergence of numerical solutions of hydrostatic and osmotic pressures to the steady-state solution obtained by Facchini et al.26 The numerical solutions are obtained for a two-layer membrane with discontinuous (i.e., ε = 0) physiological parameters at the interface. The initial conditions are a linear distribution of both pressures across microvessel wall, while constant boundary conditions, equal to the initial values reported in Table 1, are applied at the luminal (r = rc) and interstitial (r = ro) interfaces. Since this distri- bution of the pressures is not consistent with the steady-state behavior, transient conditions will develop, finally converging to the steady-state solution, for which an analytical solution is avail- able.26 The dynamics of hydrostatic and osmotic pressures varies sig- nificantly during the transient generated by imposing an unrealistic linear decline across the microvessel wall. The hydrostatic pres- sure declines rapidly showing a significant difference with respect to the initial condition already at t = 10−6, while no variations are observed in the osmotic pressure (Figure 6). A longer time should be waited for, i.e. t = 10−2, in order to observe a similar variation of the osmotic pressure, when hydrostatic pressure is nearly in steady-state conditions. Both reach steady state at t = 10−1 and the numerical solution approximates very well the analytical solution obtained by Facchini et al.26 In Figure 7, the L2-norm of the difference between the transient numerical solution and the steady-state analytical solution of both osmotic and hydrostatic pressure are depicted with respect to time t. Note that the convergence is quick at the beginning (until around t = 10−4), but then it gradually slows down and convergence is obtained for t ≃10−1 (notice the logarithmic scale of the figure). Time-dependent model Since analytical solutions of equations (6) are not available, we check the accuracy of the numerical scheme as follows. First, we assume that the solutions for both hydrostatic and osmotic pres- sures are given by the following arbitrary expressions: No n c om me rci al us e o nly Article [Veins and Lymphatics 2017; 6:6849] [page 109] (22) The arbitrary increase of 5 for the osmotic pressure is neces- sary to assure that osmotic pressure remains always positive. Given their arbitrary choice, the functions (22) do not satisfy the system of equations (6), but the following modified system of differential equations, where the source terms Ŝ 1 (r,t) and Ŝ 2 (r,t) are obtained by substituting p (r, t) with p̂ (r,t) and Π(r, t) with Π̂ (r, t) into the equations (6): (23) The solution of system (23) with initial conditions p̂ (r, 0), = sin(2 π r), Π̂ (r, 0) = sin(2 π r)+5 and boundary conditions p̂ (0, t) = p̂ (1, t) =sin(2 π t), Π̂ (0, t) = Π̂ (1, t)=sin(2 π t) + 5 should be equal to the imposed functions (22). Notice that the only change in the numerical implementation lies in the right-hand sides of the discretized equations (20), which become now equal to S1(i,n)+4(Δx)2 (ri+ξ) Ŝ1 (ri,tn+1) and S2(i,n) +8(Δx)2 (ri+ξ) Ŝ 2 (ri,tn+1), respectively, where S1(i,n) and S2(i,n) are the source terms (21) of the discretized equations (20). Since equations (22) are periodic in t with period equal to 1, we evaluate the accuracy of the numerical method by comparing the exact solutions (solid curves) and the numerical solutions (circles) after one period, as depicted in Figure 8. The imposed and the numerical solutions of the modified system of equations (23) are barely distinguishable, thereby confirming that the implicit Crank-Nicolson finite-difference scheme pro- posed by Freeze46 works well also when applied to transient equa- tions (6). A more comprehensive analysis of the numerical error is provided in Table 3, which shows the most used norms of the error at time t = 1, defined as the difference between the numerical solu- tion and the imposed functions (22) at the grid nodes. Table 4 shows the empirical order of convergence of the numerical scheme to the exact solution. As expected, it approaches the theoretical order of accuracy of Crank-Nicholson method, Table 3. Lp-norm errors at t = 1 of the two-layer model with Dt = (Dx)2 and ε2 = 10−2, numerically solved by varying the number of grid nodes as reported in the first column. Mesh Dx Dt L∞-error L1-error L2-error 7 0.166667 0.0277778 0.437544 0.284086 0.276584 13 0.0833333 0.00694444 0.116974 0.0636364 0.0682671 25 0.0416667 0.00173611 0.0299723 0.0175038 0.0178583 49 0.0208333 0.000434028 0.00734658 0.00463064 0.00446781 97 0.0104167 0.000108507 0.0018257 0.00118912 0.00112179 Table 4. Empirical orders of accuracy at time t = 1 for the two-layer model with Dt = (Dx)2 and ε2=10−2, according to the Lp-norm errors reported in Table 3. Mesh Dx Dt L∞-error L1-error L2-error 13 0.0833333 0.00694444 1.90324 2.1584 2.01846 25 0.0416667 0.00173611 1.96448 1.86218 1.9346 49 0.0208333 0.000434028 2.02849 1.91838 1.99895 97 0.0104167 0.000108507 2.00862 1.96132 1.99377 Figure 6. Hydrostatic (A) and osmotic (B) pressures for the two-layer model at given times t = {0,10−6,10−4,10−2,10−1}, obtained numer- ically with a spatial grid of 1000 points, Dt = (Dx)2, ε = 0 and tolerance equal to 10−9. The analytical steady-state solutions (thick curves) for both hydrostatic and osmotic pressures obtained by Facchini et al.26 are also shown and compared with the numerical solution at time t = 10−1 (black bullets) when steady state is reached. A B No n c om me rci al us e o nly Article which is second order. This occurs also for ∆t = k(∆x)2 with k larger than 1, but smaller than k = 1000. It can be shown that, for k ≥ 1000, the numerical solutions are no longer convergent to their analytical counterparts. Steady-state implicit exact solutions Single-layer solutions For a single-layer vessel wall with uniform physiological parameters (ℓp, ℓd, σ) and constant pressures in the lumen and the interstitium, the solution of the stationary counterpart of the system of equations (6), obtained by letting ∂p/∂t = 0 and ∂Π/∂t = 0, leads to the following expressions:26 (24) for the hydrostatic pressure and (25) for the osmotic pressure, respectively. In equations (24) and (25), the integration constants assume the following expressions: (26) (27) k4=fe f ξδ (28) where f = (k1/k2)(σ − 1)Πc − 1 and δ=k12 (σ-1)2/ [k2 (ℓpσ2 − ℓd)]. By imposing the osmotic pressure equal to Πc at r = rc, we obtain the following expression in the only unknown k2, provided that k1 is given by equation (26): (29) Two-layer solutions In this case, the physiological parameters (ℓp, σ, ℓd) and the integration constants (k1, k2, k3 and k4) are layer-specific. The need to respect continuity in both hydrostatic and osmotic pressures at the inner and outer surfaces of the microvessel walls, leads to the following expressions k3G, k4G and k4G, k4W to be substituted into equa- tions (24) and (25):26 (30) where (31) By imposing the continuity conditions at the interface between the two layers at r = rg, we obtain the following two equations for the two remaining unknowns (k1 and k2) in equations (24) and (25):26 (32) Figure 7. L2-norm of the difference between the transient numer- ical solution (obtained numerically with a spatial grid of 1000 points, Dt = (Dx)2, ε = 0 and tolerance equal to 10−9) and the steady-state analytical solution for the two-layer model of both the osmotic (dashed curve) and the hydrostatic pressure (solid curve). [page 110] [Veins and Lymphatics 2017; 6:6849] No n c om me rci al us e o nly Article [Veins and Lymphatics 2017; 6:6849] [page 111] (33) where g1 and g2 assume the following expressions (34) (35) with σG ≠ σW and (36) References 1. Levick JR. Capillary filtration–absorption balance reconsidered in light of dynamic extravascular factors. Exper Physiol 1991;76:825-57. 2. Levick JR. An introduction to cardiovascular physiology. Boca Raton, FL: CRC Press; 2010. 3. Starling EH. On the absorption of fluids from the connective tissue spaces. J Physiol 1896;19:312-26. 4. Katchalsky A, Curran PF. Non-equilibrium thermodynamics in biophysics. Cambridge, MA: Harvard Univ. Press; 1965. 5. Landis EM. Micro-injection studies of capillary permeability. II. The relation between capillary pressure and the rate at which fluid passes through the walls of single capillaries. Am J Physiol 1927;82:217-38. 6. Landis EM. Micro-injection studies of capillary blood pressure in human skin. Heart 1930;15:209-23. 7. Michel CC. Starling: the formulation of his hypothesis of microvascular fluid exchange and its significance after 100 years. Exper Physiol 1997;82:1-30. 8. Weinbaum S. 1997 Whitaker distinguished lecture: models to solve mysteries in biomechanics at the cellular level; a new view of fiber matrix layers. Ann Biomed Engine 1998;26:627- 43. 9. Levick JR, Michel CC. Microvascular fluid exchange and the revised Starling principle. Cardiovasc Res 2010;87:198-210. 10. Sugihara-Seki M, Fu BM. Blood flow and permeability in microvessels. Fluid Dyn Res 2005;37:82-132. 11. Weinbaum S, Tarbell JM, Damiano ER. The structure and func- tion of the endothelial glycocalyx layer. Ann Rev Biomed Engine 2007;9:121-67. 12. Ebong EE, Macaluso FP, Spray JM, Tarbell DC. Imaging the endothelial glycocalyx in vitro by rapid freezing/freeze substi- tution transmission electron microscopy. Arterioscler Thromb Vasc Biol 2011;31:1908-15. 13. Squire JM, Chew M, Nneji G, et al. Quasi-periodic substructure in the microvessel endothelial glycocalyx: a possible explana- tion for molecular filtering? J Struct Biol 2001;136:239-55. 14. Hawkins BT, Davis TP. The blood-brain barrier/neurovascular unit in health and disease. Pharmacol Rev 2005;57:173-85. 15. Formaggia L, Quarteroni A, Veneziani A. Cardiovascular math- ematics. Milan: Springer; 2009. 16. Quarteroni A, Veneziani A, Zunino P. Mathematical and numer- ical modeling of solute dynamics in blood flow and arterial walls. SIAM J Numer Analys 2002;39:1488-511. 17. Prosi M, Zunino P, Perktold K, Quarteroni A. Mathematical and numerical models for transfer of low-density lipoproteins through the arterial walls: a new methodology for the model set up with applications to the study of disturbed lumenal flow. J Biomechan 2005;38:903-17. 18. Causin P, Guidoboni G, Malgaroli F, et al. Blood flow mechan- ics and oxygen transport and delivery in the retinal microcircu- lation: multiscale mathematical modeling and numerical simu- lation. Biomechan Model Mechanobiol 2016;15:525-42. 19. Kedem O, Katchalsky A. Thermodynamic analysis of the per- Figure 8. Hydrostatic (A) and osmotic (B) pressures of the modified problem for the two-layer model after one period (t = 1), obtained numerically with a spatial grid of 49 points, Dt = 5 · 10−3 and ε2 = 10−2 (circles), compared to the analytical solutions (solid curves). A B No n c om me rci al us e o nly Article [page 112] [Veins and Lymphatics 2017; 6:6849] meability of biological membranes to non-electrolytes. Biochim Biophys Acta 1958;27:229-46. 20. Khakpour M, Vafai K. A comprehensive analytical solution of macromolecular transport within an artery. Int J Heat Mass Transfer 2008;51:2905-13. 21. Khakpour M, Vafai K. Analysis of transport phenomena within PEM fuel cells - an analytical solution. Int J Heat Mass Transfer 2008;51:3712-23. 22. Khakpour M, Vafai K. Effects of gender–related geometrical characteristics of aorta-iliac bifurcation on hemodynamics and macromolecule concentration distribution. In J Heat Mass Transfer 2008;51:5542-51. 23. Ai L, Vafai K. A coupling model for macromolecule transport in a stenosed arterial wall. Int J Heat Mass Transfer 2006;49:1568-91. 24. Khakpour M, Vafai K. Critical assessment of arterial transport models. Int J Heat Mass Transfer 2008;51:807-22. 25. Michel CC, Curry FE. Microvascular permeability. Physiol Rev 1999;79:703-61. 26. Facchini L, Bellin A, Toro EF. A mathematical model for filtra- tion and macromolecule transport across capillary walls. Microvasc Res 2014;94:52-63. 27. Bear J. Hydraulics of groundwater. New York, NY: McGraw- Hill International Book Company; 1979. 28. Charm SE, Kurland GS. Blood flow and microcirculation. New York: Wiley; 1974. 29. Adamson RH, Lenz JF, Zhang X, et al. Oncotic pressures opposing filtration across non–fenestrated rat microvessels. J Physiol 2004;557:889-907. 30. Michel CC, Phillips ME. Steady-state fluid filtration at different capillary pressures in perfused frog mesenteric capillaries. J Physiol 1987;388:421-35. 31. Hu X, Weinbaum S. A new view of Starling’s hypothesis at the microstructural level. Microvasc Res 1999;58:281-304. 32. Speziale S, Tenti G, Sivaloganathan S. A poroelastic model of transcapillary flow in normal tissue. Microvasc Res 2008;75:285-95. 33. Urick RJ. A sound velocity method for determining the com- pressibility of finely divided substances. J Appl Phys 1947;18:983-7. 34. Carew TE, Vaishnav RN, Patel DJ. Compressibility of the arte- rial wall. Circulation Res 1968;23:61-8. 35. Robinson TC. Arterial graft prosthesis. US Patent 4731073; 1988. Available from: http://www.patentlens.net/ patentlens/patent/US_4731073/en/ 36. Kawamoto A, Matsumori T, Yamasaki S, et al. Heaviside pro- jection based topology optimization by a PDE–filtered scalar function. Struct Multidisciplin Optimizat 2011;44:19-24. 37. Corless RM, Gonnet GH, Hare DEG, et al. On the Lambert W function. Adv Computat Math 1996;5:329-59. 38. Barry DA, Parlange JY, Li L, et al. Analytical approximations for real values of the Lambert W-function. Math Comput Simulat 2000;53:95-103. 39. Facchini L, Bellin A, Toro EF. Modeling loss of microvascular wall homeostasis during glycocalyx deterioration and hyperten- sion that impacts plasma filtration and solute exchange. Curr Neurovasc Res 2016;13:147-55. 40. Müller LO, Toro EF. A global multi-scale model for the human circulation with emphasis on the venous system. Int J Numer Methods Biomed Engine 2013;30:681-725. 41. Müller LO, Toro EF. Enhanced global mathematical model for studying cerebral venous blood flow. J Biomechan 2014;47: 3361-72. 42. Müller LO, Toro EF, Haacke EM, Utriainen D. Impact of CCSVI on cerebral haemodynamics: a mathematical study using MRI angiographic and flow data. Phlebol 2016;31:305- 24. 43. Zamboni P, Galeotti R, Menegatti E, et al. Chronic cere- brospinal venous insufficiency in patients with multiple sclero- sis. J Neurol Neurosurg Psychiatr 2009;80:392-9. 44. Talbert DG. Raised venous pressure as a factor in multiple scle- rosis. Med Hypothes 2008;70:1112-7. 45. Toro EF. Brain venous haemodynamics, neurological diseases and mathematical modelling. A review. Appl Math Computat 2016;272:542-79. 46. Freeze RA. A stochastic-conceptual analysis of one-dimension- al groundwater flow in nonuniform homogeneous media. Water Resour Res 1975;11:725-41. No n c om me rci al us e o nly