COMPUTATIONAL MODELS FOR FLUID EXCHANGE BETWEEN MICROCIRCULATION AND TISSUE INTERSTITIUM

. The aim of this work is to develop a computational model able to capture the interplay between microcirculation and interstitial ﬂow. Such phenomena are at the basis of the exchange of nutrients, wastes and pharmacological agents between the cardiovascular system and the organs. They are particularly interesting for the study of eﬀective therapies to treat vascularized tumors with drugs. We develop a model applicable at the microscopic scale, where the capillaries and the interstitial volume can be described as independent structures capable to propagate ﬂow. We facilitate the analysis of complex capillary bed conﬁgurations, by representing the capillaries as a one-dimensional network, ending up with a heterogeneous system characterized by channels embedded into a porous medium. We use the immersed boundary method to couple the one-dimensional with the three-dimensional ﬂow through the network and the interstitial volume, respectively. The main idea consists in replacing the immersed network with an equivalent concentrated source term. After discussing the details for the implementation of a computational solver, we apply it to compare ﬂow within healthy and tumor tissue samples.

1. Introduction. The microcirculation is a fundamental part of the cardiovascular system [25], because it is responsible for mass transfer from blood to organs. Theoretical models, with a variety of different approaches, help to understand and quantify the main mechanisms at the basis of these phenomena. Mathematical models become even more significant when the attention is focused on the study of pathological conditions. A prominent case is the study of angiogenesis and the related flow and mass transport in tumors, because mass transport is one of the main vehicles to target tumors with pharmacological agents. Among the vast literature in this field, we mention the sequence of works by Baxter & Jain, [1,2,3,4], which represents a consistent theoretical treatment of transport in tumors. Pathological angiogenesis leads to a microvascular network with different characteristics than the healthy case [7]. In particular, tumor neovasculature is generally characterized by a tortuous geometry, leaky capillary walls and lack of lymphatic drainage network [7,23]. The combination of these effects, usually called enhanced permeability and retention, has a significant impact on the flow through vascularized tumors as well as on the ability of targeting tumors with drugs or nanoparticles [14].
The literature in this field is vast and heterogeneous, we don't have the ambition to provide an overview here. We rather focus on the available studies that have addressed mathematical models of microcirculation and interstitial flow. We identify two different approaches. On one hand, the role of microvessels on a macroscopic portion of tissue can be described by means of homogenization methods. In this case, the heterogeneous material combining the capillary bed and the interstitial tissue is described as a homogeneous continuum, whose effective transport properties are determined by means of averaging the solution of specific cell problems defined on a periodic unit cell. For the application of these techniques to transport in tumors we refer to [8,39]. On the other hand, several authors have addressed the same problem at the microscopic scale. Provided that the geometry of the microvessels is available, the interaction between the network of channels and the surrounding interstitial tissue can be modelled using the fundamental equations of flow. The peculiar nature of microcirculation introduces additional difficulties and allows for some simplifications at the same time. On one hand, the capillaries can be modelled as long and narrow cylindrical vessels. Then, asymptotic expansions that exploit the large aspect ratio of the channel can be derived to approximate the fluid exchange from one capillary to the surrounding tissue. This idea has been successfully exploited to study the microvascular flow in simple arrays of capillaries [5,15,16]. On the other hand, vascular networks are characterized by a complex, possibly irregular geometry. The previous semi-analytic methods may be hardly applied to realistic configurations. We believe that numerical methods may override this obstacle. For example, the method of Green's functions, which arises from the previous approach, is capable to handle complex vascular networks, because it substantially relies on numerical approximation methods. It has been extensively applied to the study of tumors [20,37,38]. The same methodology has also been recently used in [41].
Our idea stems from the latter family of methods. We aim to study the interplay of microcirculation and interstitial flow on a space scale that is sufficiently small to clearly separate the capillary bed from the interstitial tissue. Simultaneously, we aim to develop a computational model that is capable to handle realistic and extended tissue portions, where the configuration of the network is not oversimplified. Our idea is to replace the approximation method based on Green's functions with the finite element method, which will be applied to numerically solve the governing equations of the flow in the interstitial and capillary domains. By this way, we override some of the classical limitations of the previous methods, which are restricted to the solution of linear problems and to the ability of determining the exact expression of the fundamental solutions of a given differential operator.
To reduce the computational cost of the model, we apply the immersed finite element method [28,44]. More precisely, the capillary bed is modelled as a network of one-dimensional channels. Due to the natural leakage of capillaries, it acts as a concentrated source of flow immersed into the interstitial volume. This reduced modelling approach significantly simplifies the issues related to the simulation of the flow in the microvessels. The main methodological and theoretical aspects of the method have already been addressed in the works by D'Angelo [9,10,11]. However, its application to flow in tumors, accounting for realistic capillary configurations, is new. The main steps of the proposed methodology are summarized in Figure 1. In particular, the derivation of the equations on the basis of general fluid dynamics principles is addressed in Section 2 (see also Figure 1, panels (a, b)). The numerical approximation and algorithmic implementation, is addressed in Section 3 (see Figure 1, panel (c)). We finally demonstrate the ability of the model to describe microcirculation and interstitial flow in healthy and tumor tissue samples characterized by realistic transport properties (see Figure 1, panel (d)). Conclusions and future perspectives are addressed at the end of the manuscript.

2.
Model set up. We define a mathematical model for fluid transport in a permeable biological tissue perfused by a capillary network. We consider a domain Ω that is composed by two parts, Ω v and Ω t , the capillary bed and the tumor interstitium, respectively. Assuming that the capillaries can be described as cylindrical vessels, we denote with Γ the outer surface of Ω v , with R its radius and with Λ the centerline of the capillary network. Figure 1 (panels (a, b)) shows a description of the geometrical model. The radius of the vessels could be in general a function of the arc length along Λ. At this stage, any physical quantity of interest, such as the blood pressure p and the blood velocity u, is a function of space (being x ∈ Ω the spatial coordinates) and time t. These quantities obey different balance laws, depending on the portion of the domain of interest and, in general, they are not continuous at the interface between subdomains. We first address the fluid transport in each portion of Ω, then we discuss the proper interface conditions in order to close the resulting coupled differential problem. We consider the tumor interstitium Ω t as an isotropic porous medium, such as the Darcy's law applies, while we start assuming a Newtonian model for the blood flow in the capillaries. The rheology of blood is analyzed in detail in [34], where it is pointed out that non-Newtonian models may be more appropriate to describe blood behavior in particular conditions. Microcirculation is an extreme case where the size of vessels is the smallest and the effect of blood pulsation is almost negligible. We will discuss in the next sections how the blood flow model could adapted to these special conditions. An essential effect for the applications we have in mind is the lymphatic drainage. The lymphatic vessels consist of one way endothelium conduits from the peripheral tissues to the blood circulation. Excess of fluid extravasated from the blood circulation is drained by lymphatic vessels and returned to the blood stream. A functional lymphatic network rapidly removes fluid and this results in lower interstitial fluid pressure and biochemical concentration levels [2]. For this reason, interstitial flow and lymphatic drainage should be considered together. Unlike the capillary network, we don't have a geometrical description of the lymphatic vessels, so we can't directly define the interaction between the lymphatics and the tissue. Following the work by Soltani and Chen, [40], we model the lymphatic drainage as a sink term in the equation for the interstitial flow. More precisely, we assume that the volumetric flow rate due to lymphatic vessels, Φ LF , is proportional to the pressure difference between the interstitium and the lymphatics, namely Φ LF (p t ) = L LF p S V (p t − p L ), where L LF p is the hydraulic conductivity of the lymphatic wall, S V is the surface area of lymphatic vessels per unit volume of tissue and p L is the hydrostatic pressure within the lymphatic channels. As a consequence of all the modelling assumptions, the fluid problem in the entire domain Ω reads as follows: where µ and k denote the dynamic blood viscosity and the constant tumor permeability, respectively, and ρ is the blood density. At the interface Γ = ∂Ω v ∩ ∂Ω t we impose continuity of the flow: where n is the outward unit vector normal to the capillary surface. The fluid flux across the capillary wall can be obtained on the basis of linear non-equilibrium thermodynamic arguments, originally developed by Kedem and Katchalsky. In particular L p is the hydraulic conductivity of the vessel wall, R g is the universal gas constant and T is the absolute temperature. Because of osmosis, the pressure drop across the capillary wall is affected by the difference in the concentration of chemicals, namely c v − c t , where c v and c t denote the concentration in the capillaries and in the interstitium, respectively. The osmotic pressure is modulated by the reflection coefficient σ that quantifies the departure of a semi-permeable membrane from the ideal permeability (where any molecule is able to travel across the membrane without resistance). Finally, to be uniquely solvable, problem (1) must be complemented by boundary conditions on ∂Ω t and ∂Ω v . The prescription of these conditions significantly depends on the particular features of the problem at hand, as well as on the available data. For this reason, we postpone any further consideration on boundary conditions to Section 4, where we will discuss the numerical simulations and the related results.
2.1. Coupling microcirculation with interstitial flow. The previous fully three-dimensional model is able to capture the phenomena we are interested in. However, two relevant simplifications may be applied without significant loss of accuracy. At the modelling level, a quasi-static flow model can be replaced to the Navier-Stokes equations in deformable domains. More importantly, we aim to override the technical difficulties that arise in the numerical approximation of the coupling between a complex network with the surrounding volume. To this purpose, we adopt the multiscale approach developed in [9,10,11], which is inspired to the immersed boundary method.
2.1.1. An immersed boundary method for the interaction of a network with a bulk volume. The concept of the immersed boundary (IB) method applied to our case can be outlined as follows. To avoid resolving the complex three-dimensional (3D) geometry of the capillary network, we exploit the IB method combined with the assumption of large aspect ratio between vessel radius and capillary axial length. This approach is represented in Figure 1, panels (a, b). More precisely, we apply a suitable rescaling of the equations and let the capillary radius go to zero (R → 0). By this way, we replace the immersed interface and the related interface conditions with an equivalent mass source.
We denote with f the flux released by the surface Γ, which is a flux per unit area. The definition of f comes from the interface conditions that we have prescribed above. At the interface between the capillary network and the tissue we require that where u t (t, x) with x ∈ Γ is the volume averaged interstitial filtration velocity in the tissue and f (p t/v (t, x)) is a point-wise constitutive law for the capillary leakage in terms of the fluid pressure, denoted here with the shorthand symbol p t/v (t, x). The immersed boundary method is able to represent the action of f on Γ as an equivalent source term, F , distributed on the entire domain Ω. More precisely, where v plays the role of a test function in the variational formulation that will be defined with more details later on. Hence, we use the notation F (p t/v (t, x)) = f (p t/v (t, x))δ Γ , meaning that F is the Dirac measure concentrated on Γ, having density f . Proceeding along the lines of [9], when R → 0 we aim to replace the mass flux per unit area by an equivalent mass flux per unit length, distributed on the centerline Λ of the capillary network. To start with, we recall the assumption that the vessels can be represented as cylinders originated by a given mean line Λ. Let γ(s) be the intersection of Γ with a plane orthogonal to Λ, located at s and let (s, θ) be the local axial and angular coordinates on the cylindrical surface generated by Λ with radius R. We apply the mean value theorem to represent the action of F on v in (3) by means of an integral with respect to the arc length on Λ. More precisely, there existsθ ∈ [0, 2π] such that Then, we exploit the fact that capillaries are narrow with respect to the characteristic dimension of the surrounding volume. Namely, we assume that R |Ω t | 1/d where d = 2, 3 is the number of space dimensions of the model. Provided that f is a linear function or operator, we conclude that We observe that, while v(s, θ) is a smooth function that can be evaluated on Λ, the solution of the problem may not be regular enough to define the point-wise value of p t/v | Λ . For this reason, the average operator on γ(s) is still applied to p t/v , even in the limit case when R → 0. In conclusion, substituting the previous formula to equation (4) we recover 2.1.2. Models for microvascular flow. The IB method described above is naturally combined with the one-dimensional (1D) model for blood flow and transport in the cardiovascular system. The derivation of such model from the full Navier-Stokes equations is a vivid field of research. We refer the interested reader to [18,31] for an introduction and for instance to [6,17,33] for more advanced studies. For microcirculation, however, the derivation of a reduced flow model is significantly simpler than in the general case. To develop the model we rely on the following assumptions: (i) the displacement of the capillary walls can be neglected, because the pressure pulsation at the level of capillaries is small; (ii) the convective effects can be neglected, because the flow in each capillary is slow; (iii) the flow almost instantaneously adapts to the changes in pressure at the network boundaries, because the resistance of the network is large with respect to its inductance. This means that the quasi-static approximation is acceptable. As a result of that, the blood flow along each branch of the capillary network can be described by means of Poiseuille's law for laminar stationary flow of incompressible viscous fluid through a cylindrical tube with radius R. Let us decompose the network Λ into individual branches Λ i , i = 1, . . . , N . We denote with λ i an arbitrary orientation of each branch that defines the increasing direction of the arc length s i . Let λ, s be the same quantities referring to the entire network Λ. According to Poiseuille's flow, conservation of mass and momentum become, where g i is a generic source term. The governing flow equation on Λ is obtained by summing (7) over the index i. In conclusion, we now represent the blood flow in the capillary bed on its centerline Λ. The coupled problem for microcirculation and interstitial flow consists to find the pressure fields p t , p v and the velocity fields where the term f (p t/v ) accounts for the blood flow leakage from vessels to tissue and it has to be understood as the Dirac measure concentrated on Λ and having line density f . A schematic visualization of this flow is reported in Figure 1, panel (d), where black arrows depict the fluid exchange between the capillaries and the interstitial tissue. The expression of f (p t/v ) is provided using the Kedem-Katchalsky equation (2). Since in this work we focus on flow, we drop the effects of the concentration gradients across the capillary walls. In this case, the constitutive law for the leakage of the capillary walls reduces to the Starling's law of filtration, We notice that f is not a simple function, but rather an integral operator, as it includes the computation of the mean value of the interstitial pressure p t . Since the capillary bed is now approximated with its centerline, the average p v (s) coincides with the pointwise value p v (s). We will discuss later on how to approximate p t (t, s) by means of quadrature rule. We finally observe that in problem (8) the distinction between the subregion Ω t and the entire domain Ω is no longer meaningful, because Λ has null measure in R d . For notational convenience, in what follows we will then identify Ω t with Ω and Ω v with Λ.

Dimensional analysis.
Writing the equations in dimensionless form is essential to put into evidence the most significant mechanisms governing the flow through healthy and tumor tissue. We first identify the characteristic dimensions of our problem. We choose length, velocity and pressure as primary variables for the analysis. The characteristic length, d, is the average spacing between capillary vessels, the characteristic velocity, U , is the average velocity in the capillary bed and the characteristic pressure, P , is the average pressure in the interstitial space. Estimates of these values are reported in [2] for healthy tissue. The dimensionless form of (8) is In the Poiseuille's equation we use the non dimensional radius R = R/d. The dimensionless groups affecting our equations are the following, which represent the hydraulic conductivity of the tissue, the non dimensional lymphatic drainage, the hydraulic conductivity of the capillary walls and the hydraulic conductivity of the capillary bed, respectively. We refer to section 4 for an estimate of these dimensionless groups magnitude and the related discussion.
3. Numerical approximation. For complex geometrical configurations explicit solutions of problem (10) are not available. Numerical simulations are the only way of applying the model to real cases. Besides applications, the study of numerical approximation methods for problem (10) requires first to address existence, uniqueness and regularity of the exact solutions and then to analyse the accuracy of the proposed scheme. These topics, already addressed for a similar problem setting in [10,11], are particularly relevant in this case because they inform us about the ability of the scheme to approximate the quantities of interest for applications. The solution of problem (10) does not satisfy standard regularity estimates, because the forcing term of equation (10)(a) is a Dirac measure. To characterize the regularity of the trial and test spaces we resort to weighted Sobolev spaces. More precisely, let us denote with L 2 α (Ω) the Hilbert space of measurable functions such that is the space of functions whose generalized derivatives up to order m belong to L 2 α (Ω) and we adopt the standard notation for the related weighted norms and semi-norms. Let now V t = H 1 α (Ω) be the natural trial space for the problem in the interstitium and V v,0 be the subspace of H 1 (Λ) of functions which vanish on the boundaries of Λ. We also introduce the Kontrat'ev type weighted spaces H m α (Ω) that characterize functions such that The discretization of problem (10) is achieved by means of the finite element method that arises from the variational formulation of the problem combined with a partition of the domain into small elements, as shown in Figure 1, panel (c). We multiply the first equation by a test function q t and integrate over Ω. The Laplace operator is treated using integration by parts combined, for the sake of simplicity, with homogeneous Neumann conditions on ∂Ω, while we use (6) to write We proceed similarly for the governing equation on the capillary bed. We denote with m − i , m + i the extrema of Λ i oriented along the arc length and with λ ± i (m i ) the reference outgoing unit vectors of Λ i at those points. After integration by parts on each branch Λ i separately, we obtain the following equation, for any i = 1, . . . , N where p v,0 denotes the lifting of nonhomogeneous Dirichlet boundary data for the capillary network. Let i ∈ J be the indices that identify the branches with the common node m j . The flow balance at m j implies that Reminding that the test functions for the pressure field on the capillary bed are continuous on the entire network, namely q h v ∈ C 0 (Λ) because V v,0 ⊂ C 0 (Λ) on 1D manifolds, summing up the previous equations with respect to the number of branches, we obtain Then the weak formulation of (10) requires to find p t ∈ V t and p v ∈ V v,0 such that, with the following bilinear forms and right hand sides, At the discrete level, one of the advantages of our problem formulation is that the partition of the domains Ω and Λ into elements are completely independent. We denote with T h t an admissible family of partitions of Ω into tetrahedrons K ∈ T h t , where the apex h denotes the mesh characteristic size. We also assume that Ω has a simple shape, such that it can be exactly represented by a collection of elements.
Since natural boundary conditions will be applied on ∂Ω, we do not enforce any constraint on the degrees of freedom of V h t located at the boundary, but the definition of a t (·, ·) may be subject to some modifications.
For the discretization of the capillary bed, each branch Λ i is partitioned into a sufficiently large number of linear segments E, whose collection is Λ h i , which represents a finite element mesh on a one-dimensional manifold. Then, we will solve our equations on Λ h : We observe that the continuity of the discrete pressure at the junctions of the network is enforced by construction, by means of the approximation space. More precisely, we will use V h v,0 , that is the restriction of V h v to functions that vanish on the boundary of Λ, to enforce essential boundary conditions on the pressure, at the inflow and outflow sections of the capillary bed. The mesh characteristic size is denoted with a single parameter h, because we will proportionally refine both finite element spaces V h t , V h v . The discrete problem arising from (10) requires to find where the bilinear forms a t (·, ·), a v (·, ·), b Λ (·, ·) are the same as before, with the only difference that b Λ h (·, ·) is now defined over the discrete representation of the network Λ h . The solution of the problem (12) is characterized by a low regularity, namely α ∈ (0, 1). In other words, V t / ∈ H 1 (Ω). For this reason, studying the convergence properties of (12) to (13) is a challenging task. As proved in [10], the optimal convergence of the finite element method is observed when the computational mesh T h t is progressively refined as it approaches to Λ. Let h be the mesh characteristic size away from the singularity, let x K be the center of mass of the element K and let µ ∈ (0, 1] be a parameter that depends on the regularity of the solution. The desired refinement is obtained assuming that the local mesh size h K scales as d(x K , Λ) 1−µ up to the minimum size h K h 1/µ in the neighbourhood of Λ. To measure the approximation error, we need the following norm, Then, in Theorem 4.1 of [10], the following properties are proved. Problems (12) and (13) admit unique solutions in V t ×V v,0 and V h t ×V h v,0 respectively. The index α characterizing the regularity of the solution is positive and strictly less than unity, namely 0 < α ≤ t < 1. Furthermore, if the exact solution enjoys the additional regularity p t ∈ H 2 2+ (Ω) and p v ∈ H 2 (Λ) the following error estimate holds true when the numerical solution is calculated using linear finite elements, where the constants C t and C v are independent of the mesh characteristic size and the convergence rate of the scheme is optimal, namely p = 1 when a suitably refined mesh is used, that is µ ≤ (α − ). If this restriction on the mesh density is not satisfied, the convergence rate is restricted to p = α − . We conclude by observing that (14) measures the error with a very strong norm that incorporates the gradients of the solution over the fictitious region at a distance less than R from the network Λ. For our applications, we are mostly interested to control the error over the fluid flow exchanged across the capillary wall, precisely Q p v − p t L 2 (Λ) . The derivation of a specific error estimate for this quantity, possibly responding to less restrictive conditions on the mesh density than (14), is the subject of current investigation.
3.1. Algebraic formulation. We aim to study the matrix form of the variational problem (13 Figure 2. A basis function of the finite element space for the interstitium, ψ i t , is depicted on the left. This is a classical piecewise linear and continuous hat function. On the right, we show the basis functions ψ i v used for the approximation on the capillary network. Note that, by construction, those functions are continuous at the junctions of network branches.
3D and 1D meshes do not conform. A sketch of the approximation functions used on the interstitium and on the capillary bed is provided in Figure 2.
where ψ j t is the average of ψ j t , according to (5). The above equations form the following linear system, where the components of the right hand side vectors are respectively Matrices B tt , B tv , B vt and B vv are built using interpolation and average operators. In particular, we define a discrete operator able to extract the mean value of ψ i t and another one able to interpolate between V h t and V h v . For every node s k ∈ Λ h we define T h γ (s k ) as the discretization of the perimeter of the vessel γ(s k ). For simplicity, we assume that γ(s k ) is a circle of radius R defined on the orthogonal plane to Λ h at point s k . This set of points is used to interpolate the basis functions ψ i t . Let us introduce a local discrete interpolation matrix Π γ (s k ) which returns the values of each test function ψ i t on the set of points belonging to T h γ (s k ). Then, we consider the average operator π vt : V h t → V h v such that q t = π vt q t . The matrix that corresponds to this operator belongs to R N h v ×N h t and it is constructed such that each row is defined as, where w are the weights of the quadrature formula used to approximate q t = 1 2πR 2π 0 q t (s)Rdθ on the nodes belonging to T h γ (s k ). The discrete interpolation operator π vt : V h t → V h v returns the value of each basis function belonging to V h t in correspondence of nodes of V h v . In algebraic form it corresponds to an interpolation matrix Π vt ∈ R N h v ×N h t . Using these tools we obtain, For the assembly of (15) we use a code developed in GetFem++, a general purpose C++ finite element library [32]. To solve system (15) we apply the GMRES method with incomplete-LU preconditioning. We perform an analysis of the computational cost of the different parts of the algorithm when the characteristic size of both the 3D and the 1D computational meshes is proportionally decreased, obtaining the data reported in Figure 3. It shows that, when the size h becomes small, the major computational time is taken from the construction of the interpolation matrices Π vt and Π vt . This is a very interesting observation, because it reveals that computational issues may arise when dealing with the interaction of two geometrical structures, such as the mesh for the network and the bulk volume. Although this is not a severe limitation in our case, because we use moderately small domains and meshes, it may become a problem of paramount importance for larger computations. To override these drawbacks, the development of the finite element solver should be complemented by expertise in the field of computational geometry, in particular for the application of efficient search algorithms and data structures to build the interpolation matrices defined above.

Preliminary validation.
Let us consider Darcy's equation for the interstitium in a slightly simpler setting than problem (1), where in particular the lymphatic term is neglected. An alternative way to set up and solve the problem has been been studied in [8,38,41], starting from the Green's representation of the Laplace equation. It is possible to define pressure solution by means of boundary potentials, using the representation formula: where G(x − y) is the fundamental solution of the Laplace equation, p 0 is a farfield interstitial pressure, in particular p t → p 0 as |x| → ∞, and Γ is the external capillary surface. The two boundary integrals are called single and double layer potential, respectively. Substituting the boundary condition on the flux (2) into the single layer potential and manipulating the double layer potential, we obtain an integro-differential representation formula for the exact solution, The solution of this model in a simple configuration featuring the flow along a single linear capillary with prescribed boundary conditions is studied in [8]. Our goal is to compare that solution, reported in [8], with the one obtained with our numerical method using the same geometry and parameters of [8]. Prescribed boundary conditions are a pressure drop along the capillary and an imposed pressure value on the outer tissue domain, that is, p v (s = in) = 1, p v (s = out) = 0.5 and p t = p 0 = 0 on ∂Ω. We solve problem (10) neglecting the lymphatic term and we represent the capillary pressure as a function of the arclength s, for different values of the vascular conductivityL p , whereL p = µL p L and L is the domain dimension. The results shown in Figure 4 are in excellent agreement with what is reported in Fig.7 of [8]. This allows us to conclude that our numerical solver is correct and it represents a valid alternative to other solution strategies.
4. Application to microcirculation and interstitial flow. Fluid and mass transport within a tumor are governed by a subtle interplay of sinks and sources, such as the leakage of the capillary bed, the lymphatic drainage and the exchange of  . From left to right, rat93 geometry, its random perturbation and rat98 geometry. fluid with the exterior volume. The aim of this section is to apply the computational model (10) to investigate how these effects influence significant and measurable quantities characterizing the flow into a vascularized tumor mass. 4.1. Available data. We aim to analyse fluid transport through tumor tissue invivo. For the reconstruction of the geometrical model of the capillary bed, we use the available data for a R3230AC mammary carcinoma in rat dorsal skin flap preparation, available in [36]. We consider two datasets, obtained with independent experiments. The first one, labelled as rat93 shows the microvascular structure over a region with overall dimensions 250 × 370 × 200µm. We refer the interested reader to [37] for more details about this experiment. The average radius of the capillary vessels is assumed to be constant and set to R = 7.64 µm. The characteristic length of the problem is chosen as the average spacing between the capillaries, d = 50µm, according to what reported in [27]. The second case consists of the vascular network on a wider sample of dimensions 550 × 520 × 230µm and it is labelled as rat98. The characteristic size D of the considered tissue samples is thus in the range of 500µm. The details of the preparation are reported in [38].
According to [2] and [22], a vascularized healthy tissue is characterized by an average interstitial pressure P = 1mmHg and by a characteristic flow speed in the capillary bed of U = 100µm/s. In the healthy case, the parameters that characterize the transport properties of the tissue are the hydraulic conductivity of the interstitium, k = 10 −18 m 2 , the hydraulic conductivity of the capillary walls, L p = 10 −12 m 2 s/kg and the plasma viscosity µ = 4 × 10 −3 kg/(ms). The magnitude of the lymphatic drainage, modeled as a distributed sink term, is estimated in [2] to be L LF p S V = 0.5 (mmHg h) −1 . Given these data, we quantify the magnitude of the dimensionless groups reported in (11). We obtain the following values, κ t = 2 × 10 −5 , Q LF = 5.2088 × 10 −5 , Q = 9.6007 × 10 −7 , κ v = 2.6759. (18) Since κ v , the dimensionless conductivity of the capillary bed, is significantly larger than the other quantities, we infer that, as expected, the transport in the coupled capillary/interstitial medium is dominated by the flow in the vascular network. More interestingly, we observe that the other dimensionless numbers lay in a similar range, namely 10 −6 Q, κ t , Q LF 10 −5 . This suggests that the interstitial flow, the leakage of the capillary bed and the lymphatic drainage have comparable effects on the interstitial flow and pressure. The significance of model (10) is the ability to capture the interplay between these phenomena. The aim of the forthcoming sections is indeed to use the model to analyze these effects in different conditions, representing for instance healthy and tumor tissue.

4.2.
Influence of the boundary conditions. The samples rat93 and rat98 represent microscopic regions separated from the surrounding tumor mass by artificial planar sections. An appropriate modelling of boundary conditions is required.
For the capillary flow, we aim to enforce a suitable pressure gradient along the network. Observing that the inflow and outflow sections of the network lay on the lateral side of the tissue slab, see Figure 1 (d), we enforce a given pressure p in on two adjacent faces and a pressure p out on the opposite ones. By this way, the pressure drop p in − p out is enforced at the tips of the network. To estimate the magnitude of the pressure drop, we use Poiseuille's law to fit a given value of the average blood velocity through healthy microvascular network, namelyū v = 0.2 mm/s. More precisely, using equation (8) which provides a pressure drop equal to p in − p out = 0.4056 mmHg for rat93 and p in − p out = 1.2522 mmHg for rat98.
For the interstitial flow, we aim to model the in-vivo configuration, where the available tumor sample is embedded into a similar environment. To represent this case, we believe that the most flexible option is to use Robin-type boundary conditions for the interstitial pressure, In equation (19), p 0 represents far field pressure value, while β can be interpreted as an effective conductivity accounting for layers of tissue surrounding the considered sample. Assuming that the interstitial pressure decays from p t to p 0 over a distance comparable to the sample characteristic size, D, dimensional analysis shows that a rough estimate of the conductivity is β = κ t /D. The specific aim of this section is to test the sensitivity of the model to variations of the parameter β over a few orders of magnitude around the reference value κ t /D, which is equivalent to 10 −6 in dimensionless form.
The results of the simulations obtained using the values β = {0, 10 −2 , 1, 10 2 } × (κ t /D) and the rat98 geometry are reported in Figure 6. The analysis of the interstitial pressure field p t shows that the results obtained using β = 0, that is homogeneous Neumann boundary conditions, or β = 10 −2 κ t /D are almost equivalent (not shown). As confirmed by the top panel of Figure 6, more significant differences are observed when β is further increased. In those cases, the boundary of the domain clearly feels the influence of the reference pressure p 0 , which is weakly enforced in proximity of the boundary. Anyway, this analysis leads us to conclude that, in the case of healthy tissue, boundary conditions mildly affect the interstitial pressure distribution. We claim that the sensitivity with respect to boundary conditions is mitigated by the presence of the uniformly distributed lymphatic drainage effect, which removes the fluid in excess released by the leaky capillary network. To confirm this hypothesis, the same set of numerical experiments has been performed for a modified model where the lymphatic drainage has been switched off. This is, in fact, the assumption that is usually adopted to model tumors. The results are reported on the bottom panel of Figure 6. A remarkable difference is observed with respect to the previous case. Now, the physiological capillary leakage can only be balanced by the flow exchanged through the external boundary. As a consequence of that, the interstitial pressure field is completely saturated when homogeneous Neumann conditions are enforced on the boundary, see Figure 6 (bottom-left). For this reason, it is essential to correctly capture the fluid flow through the artificial boundaries of the domain. According to the results reported here and in the forthcoming section, we believe that the range β ∈ (1, 10 2 )×(κ t /D) is the most adequate for this purpose.

4.3.
Comparison of flow in healthy and tumor tissue models. The natural application of the model is a comparison of flow in the healthy and pathological conditions. To pursue this aim, we compare the following cases: A, healthy tissue. With respect to our model, this case is defined by: -a normal capillary network configuration. To match this condition, we use the network rat93 available from [36], which shows a smooth and regular ramification of capillaries. -a tortuous capillary network configuration that is obtained in our case by means of a random perturbation of the points between the segments of the rat93 geometry, shown in Figure 5 (central geometry). To obtain a significant difference, the amplitude of the perturbations is adjusted such that the total length of the network almost doubles with respect to the healthy case.  In addition to these cases, we consider some intermediate configurations that will help us to highlight the competing effects of enhanced permeability and lymphatic drainage. The first one, labelled as case C below, represents the properties of the tumor treated with a vascular re-normalization therapy: C, tumor after vascular re-normalization therapy. The main characteristic of the model are reported from [24]: -for the re-normalized capillary bed geometry we use the rat93 data.
-a tumor capillary phenotype is assumed to be normal after the therapy.
-absence of lymphatic drainage function.
Keeping in mind that it does not correspond to an observed physiological state, it will be interesting to compare Case C with the dual one, which arises from the tumor model, where only the lymphatic drainage is restored to the healthy state. Furthermore, to achieve a more direct comparison with case C, we use the smooth vascular geometry in this case too. More precisely, the configuration, labelled as D, is defined as follows: D, tumor with active lymphatic drainage. This idealized case is obtained by: -the same capillary network used for the healthy case.
-a tumor capillary phenotype and corresponding wall conductivity.
-presence of healthy lymphatic drainage function.
An extensive comparison of the flow indicators for the test cases A, B, C and D, is reported in Table 1 below. The numerical experiments are also repeated for different types of boundary conditions on the artificial sections ∂Ω, namely, we vary the parameter β as β ∈ {0, 1, 10 2 }κ t /D.

4.3.1.
Analysis of blood flow in the capillary bed. Blood flow in the capillary bed is basically represented by the average value of the velocity in the network, that is directly computed from the pressure field in the network as follows, ∂p v ∂s λ · λds .
With respect to this quantity, from Table 1 we observe that cases A, C and D, which are characterized by the same network geometry, are almost equivalent for all numerical experiments, while the same quantity for case B is basically halved. These results provide a strong evidence that the blood filtration in the capillary network is inversely proportional to the total length of the network, as shown by the expression above, while it is almost insensitive to all the other variables of the problem. This conclusion is confirmed by the analysis of the dimensionless groups characterizing the problem, namely (18). Since κ v κ t , Q, Q LF , the flow problem in the capillary bed is decoupled from the one in the interstitial tissue. In other words, the feedback of the interstitial fluid pressure on the capillary network is almost negligible. Because of leakage, the network substantially acts as a source term on the tissue.
Looking at the spatial variation of the blood velocity on Λ, we observe that it is almost constant over the network, as a consequence of the fact that the pressure p v is linearly decreasing from p in to p out . This behavior is due to the fact that equation (8) c is linear and the leakage effect is small. Although this model is very popular for microcirculation, see for instance [1,8,10,35,39,41], it is affected by some limitations. Using a linear model implies that the capillary flow is not sensitive to the tortuosity of the network, which could be quantified in our case by the magnitude of the angles between the individual branches. Another limitation is that the presence of red blood cells is only indirectly accounted for, by suitably tuning the viscosity of the fluid. Although the full three-dimensional resolution of the fluid particle interaction, addressed for instance in [13,29,30], would be too demanding for our purposes, other reduced microcirculation models, such as [12,26], should be in future compared to the present approach.

4.3.2.
Quantitative indicators of interstitial flow. Interstitial flow directly affects how efficiently nutrients, drugs and wastes are transported to and from cells among the interstitial tissue. To analyze these effects we introduce two quantitative indicators: the net fluid flux from the capillary network to the interstitial volume and the equivalent conductivity of the tissue construct.
The local fluid exchange rate is defined as f (p t/v ) in equation (9). In Table 1 we report its total value over the network, This expression shows that f is affected by the hydraulic conductivity of the capillary walls, L p , as well as by the interstitial fluid pressure p t and because of the negative sign these quantities have a competitive role in determining the flux.
The other indicator is the norm of the diagonal hydraulic tensor K, defined according to Darcy's law for an isotropic porous construct, where u and δp are mean quantities computed over the entire tissue construct, where δp is defined as the difference from the local pressure and the basal pressure, namely δp(t, x) = p(t, x) − p 0 . The hydraulic conductivity tensor represents the ease with which a fluid can move through the medium and, according to equations (20), it is determined by the pressure drop in the construct, which in turn is affected by the conductivity of the capillaries as well as by the lymphatic drainage into the tissue. Since we are actually dealing with two different subregions, namely Ω v and Ω t , to compute proper values of u and δp we use the following definitions: |Ωt| Ωt u t dσ, Furthermore, observing that the sample rat93 features an almost planar network in x and y directions, the final form of equation (20) results in: Given K xx and K yy from the equations above, we compute a representative value of the tensor K using the Frobenius norm: The computed values for K F are reported in Table 1. High values of K F indicate that the construct could be well perfused, conversely low values of K F mean that the construct is impervious. The total fluid flux f and the hydraulic conductivity indicator K F are affected by both the capillary leakage and the interstitial fluid pressure. Understanding which of these two factors dominates in different conditions will be the objective of the forthcoming section.

4.3.3.
Analysis of capillary leakage and interstitial pressure. The effect of enhanced permeability and retention. We discuss the results of the simulations case by case, for the test cases A, B, C and D defined before. Case A. Table 1 and Figure 7 show that f is rather insensitive to the outer boundary conditions on the interstitial volume for the healthy tissue model. This is due to the lymphatic drainage effect, which removes the excess of fluid in the interstitial volume no matter of how much fluid is exchanged across the artificial sections ∂Ω.
In other words, the retention effect is absent for the healthy tissue. This allows the capillary leakage to reach its physiological range. The presence of lymphatic drainage causes the pressure field in the tissue to be very low, since the interstitial pressure nearly approaches the minimum value in all the domain, as shown in Figure 8 (top-left). The value of the K F indicator is higher than all other cases for the entire range of β and it slightly increases with β, due to the fact that the artificial boundaries become more permeable to flow.
Case B. The situation is radically different for the tumor case, because the interstitial fluid pressure is highly sensitive to the conditions that regulate the fluid exchange with the external volume. We analyze the variation of f by progressively increasing the values of β. The first dataset of Table 1 corresponds to homogeneous Neumann boundary conditions for the external boundary ∂Ω. In absence of lymphatic drainage and fluid exchange with the exterior, the interstitial fluid pressure becomes completely saturated. More precisely, we observe that the average interstitial fluid pressure on Ω, denoted as p t , reaches the average capillary pressure over the network, which is not reported but could easily be quantified as (p in +p out )/2 0.2 mmHg, since we know that p v linearly varies over Λ. As a result of that, the pressure gradient p v − p t is significantly lower than the physiological value and the corresponding flux practically vanishes, namely f = 1.0004 × 10 −25 cm 3 /s. At the same time, the local fluid exchange is not completely negligible, owing to the variations of the pressure in the capillary bed. In the worst case β = 0, the maximal flow rate is indeed equal to max Λ f (p t/v ) = 8.7531 × 10 −14 cm 2 /s. In conclusion, the fluid retention effect combined with the increased conductivity of the tumor capillary phenotype leads to interstitial fluid pressure saturation and significant reduction of the interstitial flow. As confirmed by Table 1 and Figure 7, the situation progressively improves when β increases, because the exchange of fluid with the exterior region decreases p t and restores more natural values of the capillary transmural pressure gradient p v − p t . As a consequence of that, the transmural capillary flux increases with β and it attains values larger than in Case A. This suggests that the augmented capillary hydraulic conductivity becomes the dominant factor for extravasation. Regarding the value of K F , we observe that it is lower than all the other cases, confirming that the overall interstitial flow is less than in the case of healthy tissue. Case C. The interstitial fluid pressure is again highly sensitive to the boundary condition that regulates the fluid exchange with the external volume, since the lymphatic system is absent. When β = 0 the value of f is similar to the value reached in Case B, f = 4.3887×10 −26 cm 3 /s and the mean interstitial fluid pressure becomes completely saturated reaching the mean value p t = 0.222084 mmHg. On the contrary, the behavior becomes more similar to Case A when the value of β increases, namely β ∈ {1, 10 2 }κ t /D, because the boundary conditions contrast the absence of the lymphatic drainage. Figure 7 (top) shows that the fluid flux f reaches values similar to Case A, because the conductivity of the capillary walls is equal to the one of healthy tissue. The two cases are comparable also with respect to the value of K F , although the healthy tissue model is slightly more permeable. Figure 7 (bottom) shows that K F increases together with the exchange of fluid with the exterior region, namely the parameter β, because it induces a low pressure field within the tissue, as we observe in Figure 8 (bottom-left). Case D. This last case is characterized by all the factors that increase the leakage from capillaries to tissue. Indeed we are considering the presence of the lymphatic network and we choose the high value of the vessels conductivity, equal to the tumor case. On one hand, thanks to the first effect, the interstitial fluid pressure is insensitive to the boundary conditions on the interstitial volume and there is no retention effect, as it happens for Case A. On the other hand, the high value of the hydraulic conductivity leads to a considerable interstitial fluid flux for each values of β. The inspection of the interstitial pressure field, shown in Figure 8 (bottomright), suggests that high conductivity and lymphatics drainage play against each other, causing a significant transversal pressure gradient in the neighbourhood of the capillary bed. As a result of that, the equivalent conductivity indicator K F is rather small if compared to the previous cases, see in particular Figure 7 (bottom).
In conclusion, these results highlight the importance of avoiding fluid retention in order to facilitate interstitial flow. In absence of lymphatic drainage, an appropriate fluid exchange from the tumor mass to the exterior could increase the ability of releasing therapeutic agents from the capillary network. In conditions where the drainage is limited, but not completely absent, the vascular renormalization therapy [24] has the potential to restore the physiological flow conditions.

5.
Conclusions and future perspectives. We have developed a computational model able to capture the flow through a heterogeneous system characterized by a network of leaky channels embedded into a porous medium. A model reduction technique, inspired to the immersed boundary method, allows us to achieve simulations of non trivial network geometries with a moderate computational cost. The application of the model to the flow through vascularized solid tumors has been extensively discussed. The efficiency of the method is tested by using the realistic vascular geometries reported in [37,38]. This feature suggests that the model may be successfully coupled with a dynamic model for angiogenesis, such as the one recently appeared in [43].
This work is part of a more ambitious research project, where flow models will be combined with mass transport, in order to analyze the distribution of nutrients or therapeutic agents such as drugs or pharmacologically active nanoparticles. The application of this modelling framework to the transport of small molecules, such as oxygen, has already been addressed in [37]. The treatment of vascular diseases using locally delivered nanoparticles is on the edge of biomedical research [14]. Theoretical models describing these phenomena into small arteries are currently being developed from both computational and analytical standpoints, see for example [21] and [42], respectively. We believe that the extension of those works to the level of microcirculation has a great potential with respect to biomedical applications. We also see a fertile ground for application to a variety of pathologies. For example, very similar models have been applied to study the microcirculation in the retina and in particular the correlation between ocular hemodynamics and glaucoma [19].