Modeling tissue perfusion in terms of 1d-3d embedded mixed-dimension coupled problems with distributed sources

TL;DR

This paper introduces a novel mixed-dimension modeling approach for tissue perfusion that improves accuracy and convergence in simulating fluid exchange between microvasculature and surrounding tissue, enabling larger-scale simulations.

Contribution

The method models vascular networks with embedded 1D segments and distributed sources, enhancing flux approximation and convergence over existing mixed-dimension techniques.

Findings

Error in exchange flux approximation reduced by a factor of 3 for coarse grids.

Improved convergence rates due to increased regularity of pressure solutions.

Potential for larger domain simulations with limited computational resources.

Abstract

We present a new method for modeling tissue perfusion on the capillary scale. The microvasculature is represented by a network of one-dimensional vessel segments embedded in the extra-vascular space. Vascular and extra-vascular space exchange fluid over the vessel walls. This exchange is modeled by distributed sources using smooth kernel functions for the extra-vascular domain. It is shown that the proposed method may significantly improve the approximation of the exchange flux, in comparison with existing methods for mixed-dimension embedded problems. Furthermore, the method exhibits better convergence rates of the relevant quantities due to the increased regularity of the extra-vascular pressure solution. Numerical experiments with a vascular network from the rat cortex show that the error in the approximation of the exchange flux for coarse grid resolution may be decreased by a factor of $3$. This may open the way for computing on larger network domains, where a fine grid resolution cannot be achieved in practical simulations due to constraints in computational resources, for example in the context of uncertainty quantification.

Tables

Table 1. Table 1 : Vessel specifications for the numerical experiment in Section   4.3 . The vertex indices and coordinates (in µ ​ m micrometer \mathrm{\SIUnitSymbolMicro m} ) are listed in the first two columns. Boundary conditions are given, if the vertex is a boundary vertex. The 4th and 5th columns list the segment connectivity and the corresponding segment radius (in µ ​ m micrometer \mathrm{\SIUnitSymbolMicro m} ).

index	coordinates ($\mu m$)	boundary condition	segment	radius ($\mu m$)
0	(20, 20, -10)	$v=1.0\mathrm{mm}{s}^{-1}$	(0, 1)	3.5
1	(20, 20, 50)	-	(1, 2)	3.5
2	(20, 20, 110)	${\tilde{p}}_{v,\text{out}}=2300\mathrm{Pa}$	(3, 4)	3.0
3	(60, -10, 20)	$v=0.8\mathrm{mm}{s}^{-1}$	(4, 5)	3.0
4	(60, 20, 20)	-	(5, 6)	3.0
5	(70, 40, 20)	-	(6, 7)	3.0
6	(80, 60, 20)	-	(8, 9)	3.0
7	(110, 60, 20)	${\tilde{p}}_{v,\text{out}}=3400\mathrm{Pa}$	(9, 10)	3.0
8	(80, 110, 80)	$v=0.8\mathrm{mm}{s}^{-1}$	(10, 11)	3.0
9	(80, 80, 80)	-	(11, 12)	3.0
10	(80, 70, 60)	-	(13, 14)	3.0
11	(20, 80, 20)	-	(14, 15)	3.0
12	(20, 80, -10)	${\tilde{p}}_{v,\text{out}}=2300\mathrm{Pa}$	(15, 16)	3.0
13	(60, -10, 80)	${\tilde{p}}_{v,\text{out}}=3400\mathrm{Pa}$	(16, 17)	3.0
14	(60, 20, 80)	-
15	(40, 70, 80)	-
16	(20, 80, 80)	-
17	(-10, 80, 80)	$v=0.8\mathrm{mm}{s}^{-1}$

Equations

$${\hat{q}}_m={\rho }_I{L}_{p}S\left[({\tilde{p}}_v-{\tilde{p}}_{t,\mathbb{W}}-\sigma \Delta \Pi \right],$$

$$-\frac{\partial }{\partial s}\left({\rho }_B\frac{\pi R^4}{8{\mu }_B}\frac{\partial p_v}{\partial s}\right)$$

$$-\nabla \cdot \left(\frac{{\rho }_I}{{\mu }_I}\kappa \nabla p_t\right)$$

$${\hat{q}}_m$$

$${\int }_{\Omega }f{\delta }_{\Lambda }\text{d}x={\int }_{\Lambda }f\text{d}s\forall f$$

$$-\frac{\partial }{\partial s}\left({\rho }_B\frac{\pi R^4}{8{\mu }_B}\frac{\partial p_v}{\partial s}\right)$$

$$-\nabla \cdot \left(\frac{{\rho }_I}{{\mu }_I}\kappa \nabla p_t\right)$$

$${\hat{q}}_m$$

$$\int {\text{displaylimits}}_0^{l_i}\int {\text{displaylimits}}_0^{2\pi }\int {\text{displaylimits}}_0^{\varrho (s)}{\Phi }_{{\Lambda }_i}r\text{d}r\text{d}\theta \text{d}s=l_i,$$

$$\Xi =\frac{(p_v-p_{t,\mathbb{W}})}{(p_v-p_{t,0})},$$

$${\Phi }^{\text{const}}(r)=\{\begin{array}{ll}\frac{1}{\pi {\varrho }^2}& r\le \varrho ,\\ 0& r>\varrho .\end{array}$$

$$-\frac{1}{r}\frac{\partial }{\partial r}\left(r\frac{\partial p_t}{\partial r}\right)$$

$$p_t(r)=p_{t,\mathbb{W}}-\frac{{\hat{q}}_m{\mu }_I}{2\pi {\rho }_I\kappa }\ln \left(\frac{r}{R}\right),r>\varrho ,$$

$$\frac{\partial p_t}{\partial r}$$

$$p_t(r)=-\frac{{\hat{q}}_m{\mu }_I}{2\pi {\rho }_I\kappa }\frac{r^2}{2{\varrho }^2}+C_0,r\le \varrho .$$

$$p_t(r)=\{\begin{array}{ll}{p}_{t,\mathbb{W}}-\frac{{\hat{q}}_m{\mu }_I}{2\pi {\rho }_I\kappa }\left[\frac{r^2}{2{\varrho }^2}+\ln \left(\frac{\varrho }{R}\right)-\frac{1}{2}\right]& r\le \varrho ,\\ p_{t,\mathbb{W}}-\frac{{\hat{q}}_m{\mu }_I}{2\pi {\rho }_I\kappa }\ln \left(\frac{r}{R}\right)& r>\varrho .\end{array}$$

$$p_{t,0}=p_{t,\mathbb{W}}-\frac{L_{p}R{\mu }_I}{\kappa }(p_v-p_{t,0})\Xi \left[\ln \left(\frac{\varrho }{R}\right)-\frac{1}{2}\right],$$

$$\Xi$$

$$={\left\{1+\frac{R{L}_p{\mu }_I}{\kappa }\left[\ln \left(\frac{\varrho }{R}\right)-\frac{1}{2}\right]\right\}}^{-1}(\text{for}{\Phi }^{\text{const}}).$$

$${\Phi }^{\text{cubic}}(r)=\{\begin{array}{ll}\frac{10}{3\pi {\varrho }^2}\left[\frac{2r^3}{{\varrho }^3}-\frac{3r^2}{{\varrho }^2}+1\right]& r\le \varrho ,\\ 0& r>\varrho ,\end{array}$$

$$p_t(r)=\{\begin{array}{ll}{p}_{t,\mathbb{W}}-\frac{{\hat{q}}_m{\mu }_I}{2\pi {\rho }_I\kappa }\left[\frac{r^2}{{\varrho }^2}\left(\frac{8r^3}{15{\varrho }^3}-\frac{5r^2}{4{\varrho }^2}+\frac{5}{3}\right)+\ln \left(\frac{\varrho }{R}\right)-\frac{19}{20}\right]& r\le \varrho ,\\ p_{t,\mathbb{W}}-\frac{{\hat{q}}_m{\mu }_I}{2\pi {\rho }_I\kappa }\ln \left(\frac{r}{R}\right)& r>\varrho ,\end{array}$$

$$\Xi ({\Phi }^{\text{cubic}})={\left\{1+\frac{R{L}_p{\mu }_I}{\kappa }\left[\ln \left(\frac{\varrho }{R}\right)-\frac{19}{20}\right]\right\}}^{-1}.$$

$$p_t(\mathbf{x})=\sum _{i=1}^N{p}_{t,i}(\mathbf{x}),$$

$${\hat{q}}_{m,i}=2\pi R_i{\rho }_I{L}_p(p_v-p_t({\mathbf{x}}_i)){\Xi }_i,$$

$$p_{t,j}({\mathbf{x}}_i)=\frac{1}{2\pi }{\int }_0^{2\pi }{p_{t,j}|}_{R_i}\text{d}\theta =:p_{t,j,{\mathbb{W}}_i}=:C_{j,i}$$

$$p_{t,i}({\mathbf{x}}_i)-p_{t,i,{\mathbb{W}}_i}=p_{t,i}({\mathbf{x}}_i)+\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N{C}_{j,i}-(p_{t,i,{\mathbb{W}}_i}+\sum _{\begin{array}{c}j=1\\ j\ne i\end{array}}^N{C}_{j,i})=p_t({\mathbf{x}}_i)-p_{t,{\mathbb{W}}_i},$$

$$p_{t,{\mathbb{W}}_i}:=\frac{1}{2\pi }{\int }_0^{2\pi }{p_t|}_{R_i}\text{d}\theta .$$

$${\Xi }_i({\Phi }^{\text{const}})=\frac{p_{v,i}-p_{t,{\mathbb{W}}_i}}{p_{v,i}-p_t({\mathbf{x}}_i)}={\left\{1+\frac{R_i{L}_p{\mu }_I}{\kappa }\left[\ln \left(\frac{{\varrho }_i}{R_i}\right)-\frac{1}{2}\right]\right\}}^{-1}.$$

$$-{\int }_{\partial K_{\Omega }}\left[\frac{{\rho }_I}{{\mu }_I}\kappa \nabla p_t\right]\cdot {\mathbf{n}}_{K_{\Omega },\sigma }\text{d}A={\int }_{K_{\Omega }}{\hat{q}}_m{\Phi }_{\Lambda }\text{d}x,$$

$$F_{K_{\Omega },\sigma }\approx -{\int }_{\sigma }\left[\frac{{\rho }_I}{{\mu }_I}\kappa \nabla p_t\right]\cdot {\mathbf{n}}_{K_{\Omega },\sigma }.$$

$$Q_{K_{\Omega }}={\int }_{K_{\Omega }}{\hat{q}}_m{\Phi }_{\Lambda }\text{d}x=\sum _{{\Phi }_{\Gamma }\in {\Phi }_{\Lambda }}{Q}_{\Gamma }\frac{1}{|\Gamma |}{\int }_{K_{\Omega }\cap \mathcal{S}({\Phi }_{\Gamma })}{\Phi }_{\Gamma }\text{d}x,$$

$$Q_{\Gamma }\approx {\int }_{\Gamma }{\rho }_I{L}_{p}S(p_v-p_{t,0})\Xi \text{d}s.$$

$$\sum _{\sigma \subset \partial K_{\Omega }}{F}_{K_{\Omega },\sigma }=Q_{K_{\Omega }},\forall K_{\Omega }\in {\Omega }_h.$$

Full text

1 Introduction

Most biological tissues contain embedded blood vessels that supply the cells in the extra-vascular space with oxygen and nutrients. Mass transfer to the extra-vascular space is essential for, e.g., the function of the liver, the clearance of substances in the blood by the kidneys, or the proliferation of therapeutic agents. Moreover, fluid exchange as well as transport of molecules over the capillary walls are of particular importance, if the blood-brain-barrier (BBB) is impaired, e.g. due to brain diseases like multiple sclerosis [1, 2, 3, 4] or glioblastoma [3, 5]. The accurate description of the fluid exchange is essential to analyze, e.g., advanced MR imaging techniques for the brain [6, 7]. Theoretical models can help to understand and quantify the relevant flow and transport processes in biological tissue [8, 4, 9].

Due to the large number of vessels in the microvasculature, fully three-dimensionally resolved models are only feasible for small tissue samples. Herein, we consider a class of models, for which the microvasculature is represented by a network of vessel segments embedded in the extra-vascular space, and the flow in the vessels is modeled by one-dimensional partial differential equations (PDEs), whereas the extra-vascular space is modeled as a homogenized porous medium by a three-dimensional PDE. The PDEs are coupled via source terms. Such models have been used to study the proliferation of cancer drugs [10, 11, 12], nano-particle transport in hypothermia therapy [13], oxygen proliferation [14, 15, 16, 17, 18], and contrast agent perfusion [9, 4].

This work focuses on the application of biological tissue perfusion. However, the presented method can also be applied to model other network-like structured embedded in porous media, such as plant root systems [19, 20, 21] and wells in petroleum engineering or geothermal applications [22, 23, 24, 25].

1.1 Mixed-dimension embedded model with line sources

We consider a region of biological tissue, $\Omega\subset\mathbb{R}^{3}$, comprising cells, the extra-cellular matrix, and interstitial fluid. Embedded in $\Omega$, a network structure of tubular vessel segments, $\Lambda\subset\mathbb{R}^{3}$, models the microvasculature. The domains are conceptually shown in Fig. 1. Blood is assumed to behave like a Newtonian fluid with constant density, $\rho_{B}$, and dynamic viscosity, $\mu_{B}$. Blood flow through the vessel lumen is governed by the Hagen-Poiseuille law. These strong assumptions merely simplify the following analysis, but are not required for the derivation of the presented method. The extra-vascular space is modeled as a porous medium with porosity $\phi$ and isotropic intrinsic permeability $\kappa$. The interstitial fluid flow is described by Darcy’s law. The two compartments exchange mass over the vessel walls, $\mathbb{W}$. This transmural mass exchange is commonly modeled by Starling’s law [26] as

[TABLE]

where transmural mass flux, $\hat{q}_{m}$ ($\mathrm{kg}\text{\,}{\mathrm{s}}^{-1}\text{\,}{\mathrm{m}}^{2}$), is formulated in terms of the filtration coefficient, $L_{p}$ ($\mathrm{m}\text{\,}{\mathrm{Pa}}^{-1}\text{\,}\mathrm{s}$), the cross-sectional vessel wall surface $S=2\pi R$, the mean hydraulic pressures in the lumen, $\tilde{p}_{v}$, and on the outer vessel wall surface, $\tilde{p}_{t,\mathbb{W}}$, the oncotic pressure difference, $\Delta\Pi=\Pi_{v}-\Pi_{t}$, and the osmotic reflection coefficient $\sigma\in\left[0,1\right]$. The oncotic pressure is exerted by large proteins in the blood plasma such as albumin. The difference in oncotic pressure between the vascular and the extra-vascular space effectively pulls fluid into the vessels. It can be assumed constant in the capillary bed, with $\Delta\Pi=$2633\text{,}\mathrm{Pa}$$ [27]. Furthermore, we choose $\sigma=1$, which corresponds to the vessel wall modeled as a perfect selectively-permeable membrane. Thus, the problem can be reformulated in terms of effective pressures, $p_{\alpha}=\tilde{p}_{\alpha}-\Pi_{\alpha},\;\alpha\in\{v,t\}$.

Flow in the given tissue region can be described by a coupled system of partial differential equations, [cf. 10, 11, 4],

[TABLE]

where $s$ is the local coordinate in axial direction of the vessel segments, $R$ denotes the equivalent vessel radius corresponding to the cross-sectional area of the vessel, $A_{v}=\pi R^{2}$, and $\rho_{I}$ and $\mu_{I}$ are density and dynamic viscosity of the interstitial fluid, respectively. The Dirac delta function $\delta_{\Lambda}$ with the property

[TABLE]

restricts $\hat{q}_{m}$ in Eq. 1.2b on the vessel center-line. The derivation of the formulation with line sources, Eq. 1.2, is discussed in [11] and based on scaling the vessel radius to the zero limit under the constraint of flux equality. However, solutions to Eq. 1.2 exhibit a singularity at $r=0$, which is both nonphysical and challenging for numerical methods. The author of [11, 28] proves convergence of standard finite element methods in weighted norms. The authors of [29] can show optimal convergence up to a log-factor in classical but local norms excluding a small neighborhood around the singularity. It is argued that the solution is only physically meaningful for $r>R$. However, numerical experiments with standard finite element methods show that optimal convergence rates can only be achieved, if the cell size $h_{\Omega}$ of the grid discretizing $\Omega$, is in the order of magnitude of the vessel radius $R$ and smaller. In [17], Green’s functions are used to solve a system of equations similar to Eq. 1.2. The authors exploit the superposition property of the Laplace operator. The vessel network is represented by a collection of point sources along the vessel center-lines. The solution in $\Omega$ is constructed by adding up the contributions of all point sources and a smooth correction function to satisfy the boundary conditions. However, due to global interaction of those contributions, the numerical solution of the point source strengths results in dense system matrices. In [30], the support of the line source contributions is narrowed by a smooth cut-off function and the correction function is numerically approximated, yielding sparse system matrices. This method is also known as subtraction method from the field of electroencephalography (EEG) source reconstruction [31]. A local version of such methods has been well-studied in petroleum engineering and is known as the Peaceman well model [32, 33] in its simplest form. Local approximations are difficult to construct for arbitrarily oriented wells [23, 24]. Moreover, these well models are derived in a discrete setting.

1.2 Other mixed-dimension embedded models

A different approach is taken in [34]. Instead of using line sources, the exchange term between vessels and extra-vascular space is considered on the actual surface of the cylindrical vessel tubes, thus increasing the dimension of the source term by one. As a consequence, the singularity is replaced by a smoother continuation of the pressure function for $r<R$, leading to better numerical properties. However, good approximations of the exchange flux $\hat{q}_{m}$ still requires very fine grid resolutions. There is the necessity for a method with a good approximation of the exchange term, even if a very fine grid resolution is not computationally feasible. Moreover, we are seeking a method that is robust and efficient regarding large local variations of vessel radii, which readily occur in complex vessel networks.

1.3 Scope of this work

In the following, we discuss a different formulation of Eq. 1.2, where the source term in the extra-vascular space is approximated as a distribution over a three-dimensional neighborhood of the vessel using a smoothing kernel. This novel approach for modeling tissue perfusion, is introduced in a continuous setting in Section 2. The model is discretized using a cell-centered finite volume method in Section 3. The robustness of the method and a comparison with existing methods is demonstrated and discussed in a series of numerical experiments in Section 4.

2 Mixed-dimension embedded model with distributed sources

Here, we consider the same set of equations as in Section 1.1 with a modified source term,

[TABLE]

where $\Phi_{\Lambda}$ is a set of kernel functions $\Phi_{\Lambda_{i}}$ that distribute $\hat{q}_{m}$ over a small tubular support region, $\mathcal{S}(\Phi_{\Lambda_{i}})$, with radius $\varrho(s)$, around a vessel segment $i$, such that $\Phi_{\Lambda_{i}}=0$ outside the support region. The flux $\hat{q}_{m}$ in Eq. 2.1c is now formulated in terms of the extra-vascular pressure $p_{t,0}$, evaluated at the vessel center-line. The function $\Xi=\Xi(\varrho(s),\Phi(s),\cdots)$ is a flux scaling factor depending on the support radius, and the chosen kernel function. We choose kernel functions $\Phi_{\Lambda_{i}}(s)$ with the property

[TABLE]

where $r$, $\theta$, $s$ are the radial, angular, and axial coordinate in a segment-local cylinder coordinate system, and $l_{i}$ is the length of segment $i$. On close observation, one can see that if the scaling factor is chosen as

[TABLE]

then in the limit of $\varrho\rightarrow 0$, Eq. 2.1 is equal to Eq. 1.2. Hence, the mixed-dimension embedded method with line sources is a special case of the presented method, where $\Phi_{\Lambda}:=\delta_{\Lambda}$. The formulation of [34] is obtained by choosing kernel functions that are $S^{-1}=(2\pi R)^{-1}$ on the vessel surface, and zero elsewhere. A similar formulation to Eq. 2.1 is introduced in [35] for an embedded discrete fracture model (2d-3d), but not further analyzed or numerically exploited.

2.1 Flux scaling factor $\Xi$ (single straight vessel)

For non-zero support, $\varrho>0$, smooth kernel functions regularize the pressure solution in the kernel support region. Hence, the solutions to Eq. 2.1 have no singularities. The scaling factor $\Xi$ can be analytically derived in terms of the kernel support radius $\varrho$ and the constraint Eq. 2.3 under some assumptions. First, let us consider a single straight and long vessel. Looking at a cross-sectional plane cutting through this vessel (see Fig. 2), we assume that the flow field in a small neighborhood around the vessel is strictly radial, and can thus be considered two-dimensional. Now, let us exemplarily consider the following radial kernel function

[TABLE]

The flux outside the kernel support region, $r>\varrho$, is identical to the flux of the formulation with line sources, as the same total mass is injected in both formulations if the flux scaling factor $\Xi$ is chosen as in Eq. 2.3. The pressure solution has to satisfy

[TABLE]

The pressure outside the kernel support region admits the following solution

[TABLE]

derived from the fundamental solution of the Laplace equation. Note that the mean pressure on the vessel surface, $p_{t,\mathbb{W}}$, is equal to the pressure evaluated at $r=R$, if $\varrho\leq R$. To derive the pressure solution inside the kernel support, we apply Leibniz’s rule to Eq. 2.5

[TABLE]

and integrating once more yields

[TABLE]

It is easy to verify that for $r=\varrho$, Eqs. 2.6 and 2.8 have equal derivatives which ensures flux continuity. The integration constant $C_{0}$ is determined such that the pressure is continuous at $r=\varrho$, yielding

[TABLE]

This function is shown qualitatively in Fig. 2. To derive a suitable $\Xi$, let us first evaluate Eq. 2.9 at $r=0$, so that $p_{t,\mathbb{W}}$ is expressed in terms of $p_{t,0}$,

[TABLE]

where $\hat{q}_{m}$ was replaced by inserting Eq. 2.1c. It directly follows from Eqs. 2.10 and 2.3 that the flux scaling factor can be expressed independently of the pressure,

[TABLE]

The flux factor has a very similar structure, when derived for other kernel functions. For instance, for the cubic kernel function

[TABLE]

we obtain

[TABLE]

and

[TABLE]

The dimensionless group $\Theta=RL_{p}\mu_{I}\kappa^{-1}$ is the ratio of the hydraulic conductivity of the vessel wall to the hydraulic conductivity of the extra-vascular space. If the filtration coefficient $L_{p}$ is low relative to $\kappa$, the regularized pressure profile is rather flat, so that the difference between $p_{t,R}$ and $p_{t,0}$ is low. This is reflected in a flux scaling factor close to $1$. If the filtration coefficient is elevated, e.g. when the BBB is impaired in a tumor, $\Theta$ can be larger than $1$, hence $\Xi$ may significantly defer from $1$. Comparing Eqs. 2.15 and 2.12, we observe that a kernel function with higher weights towards the vessel center-line tends to result in a $\Xi$ further away from $1$.

2.2 Kernel support radius $\varrho$

Considering that the physically meaningful part of the pressure solution is actually given for $r\geq R$, if $\varrho<R$, our formulation gives identical solutions (for $r\geq R$) to the methods discussed in [11, 34]. However, the smoothness of the solution for $r<R$ can be controlled by the choice of kernel functions. For $\varrho>R$, the physical pressure solution $p_{t}$ is altered. However, the source term as well as $p_{v}$, remain the same for the derived $\Xi$. For $\varrho>R$, the action of the coupling term on $p_{t}$ is distributed over an artificially enlarged vessel volume.

2.3 Multiple vessels

In the capillary bed, vessels form a dense network. The vessel diameters are small, such that average distances between vessels are about one order of magnitude larger than the radii, and the vessel volume fraction is in the range of $25\text{\,}\mathrm{\char 37\relax}$. In the following, we present our considerations for choosing $\Xi$ in such a network.

To this end, let us consider a setup of $N$ parallel long vessels with given boundary conditions, such that the resulting pressure solution in each plane, $\mathbb{P}\in\mathbb{R}^{3}$, perpendicular to the vessel is equal to the two-dimensional problem on that plane (the in-plane solutions are independent of the solution in any parallel plane). Due to the linearity of the Laplace operator the solution can be split into contributions by the individual vessels (superposition principle),

[TABLE]

where $r_{i}=||\boldsymbol{x}_{i}-\boldsymbol{x}||_{2}$ is the distance between a point $\boldsymbol{x}\in\mathbb{P}$ and the center $\boldsymbol{x}_{i}$ of the vessel $i$. The source term contributions of the individual vessels are denoted by $\hat{q}_{m,i}$. Note that it depends on $p_{t}$ rather than $p_{t,i}$,

[TABLE]

so that the partial solutions $p_{t,i}$ are not independent of each other. Each partial solution, $p_{t,i}$, assumes the same form as Eq. 2.9. Furthermore, $p_{t,i}$ is a harmonic function, satisfying the Laplace equation $\nabla\cdot\!\nabla\!p_{t,i}=0$, for $r_{i}>\varrho_{i}$. Assuming, that the kernel support regions of two neighboring vessels do not overlap, we observe that

[TABLE]

using the mean value property of harmonic functions [36, p.4f], where $p_{t,j,\mathbb{W}_{i}}$ is defined as the average of $p_{t,j}$ over the cross-sectional vessel surface of vessel $i$. From Eqs. 2.16 and 2.18 follows that

[TABLE]

where

[TABLE]

This allows us to derive the flux scaling factors $\Xi_{i}$ analogously to the single vessel case (cf. Eqs. 2.10 and 2.12),

[TABLE]

In Section 4.2, we construct an analytical solution with multiple parallel vessels using these observations.

Unfortunately, the derivation using the mean value property of harmonic functions is no longer valid for arbitrarily-oriented vessels. In a fully three-dimensional setup, the mean value property concerns the integral over a sphere, while in our model we are integrating over the boundary of a circular cross-section. However, numerical experiments have shown that using an unmodified flux scaling factor does not introduce a significant model error. A possible explanation is that at a large enough distance to the neighboring vessels the integral over the circular cross-section is expected to be very close to the integral over a sphere with the same radius and center point, given that the pressure gradient decays with $r_{i}^{-1}$. Moreover, realistic segmented vessel networks contain vessel bifurcations and may contain bends with sharp angles. In the vicinity of such features, the kernel support regions of two connected vessels may overlap. In such a case the correct flux scaling factor is not easily determined. However, we herein still employ the flux scaling factor shown in Eq. 2.12, assuming that possible errors are small, only occur in the close vicinity of such features and that the influence on the global solution is negligible in realistic applications. These assumptions and the corresponding modeling errors are investigated in a numerical experiment in Section 4.4.

3 Discretization

We discretize Eq. 2.1 using a cell-centered finite volume method with two-point flux approximation. To this end, $\Omega$ and $\Lambda$ are decomposed into two independent meshes $\Omega_{h}$ and $\Lambda_{h}$ consisting of control volumes $K_{\Omega}\in\Omega_{h}$ and $K_{\Lambda}\in\Lambda_{h}$. Herein, $K_{\Omega}$ are chosen to be hexahedra and $K_{\Lambda}$ are line segments. Furthermore, the control volume boundary, $\partial K$, can be split into a finite number of faces $\sigma\subset\partial K$, such that $\sigma=K\cap L$, given a neighboring control volume $L$. Integrating Eq. 2.1b over a control volume $K_{\Omega}$ and applying the Gauss divergence theorem to the left hand side yields

[TABLE]

where $\boldsymbol{n}_{K_{\Omega},\sigma}$ is the unit outward-pointing normal on face $\sigma\subset\partial K_{\Omega}$. The exact fluxes are approximated by numerical fluxes

[TABLE]

Let $\mathcal{I}$ be the set of intersections $\Gamma=K_{\Omega}\cap K_{\Lambda}$. Furthermore, let $\Phi_{\Gamma}\in\Phi_{\Lambda}$ denote a kernel function with the support $\mathcal{S}(\Phi_{\Gamma})$ associated with $\Gamma$. The discrete source term is computed as

[TABLE]

where $Q_{\Gamma}$ is the numerical approximation of the source term integrated over the intersection $\Gamma$,

[TABLE]

Hence, Eq. 3.1 can be reformulated as

[TABLE]

The discretization of the domains $\Omega$ and $\Lambda$ are visualized in Fig. 3, including an illustration of the kernel support region.

We mention that due to the approximation of the vessel segments as cylindrical tubes around the center-line, segments may overlap, in particular, at bifurcations and bends with sharp angles. To the knowledge of the authors, the resulting discretization errors are neglected throughout the present literature. This is a fair assumption, given that these overlaps are small and the corresponding discretization error is expected to be small in comparison with errors resulting from the vessel network segmentation.

3.1 Numerical fluxes $F_{K_{\Omega},\sigma}$

We compute the numerical fluxes using a two-point flux approximation,

[TABLE]

for two neighboring control volumes $K$ and $L$, with the pressure degrees of freedom $p_{K_{\Omega}}$ and $p_{L_{\Omega}}$ associated with these control volumes, and with the transmissibilities

[TABLE]

where $\boldsymbol{d}_{K_{\Omega},\sigma}=\boldsymbol{x}_{\sigma}-\boldsymbol{x}_{K_{\Omega}}$ is a vector from the center of the control volume $K_{\Omega}$ to the center of the face $\sigma$, and $|\sigma|$ denotes the area of face $\sigma$.

3.2 Numerical fluxes $F_{K_{\Lambda},\sigma}$

Analogously to the derivation above, a discrete representation of Eq. 2.1a is given by

[TABLE]

with

[TABLE]

where $p_{\sigma}$ denotes the pressure at $x_{\sigma}$, and the transmissibilities are defined analogously to Eq. 3.7,

[TABLE]

with $R_{K}$ denoting the vessel radius of control volume $K$. Since $\Lambda_{h}$ consists of a network of segments $K_{\Lambda}$, it occurs that a face $\sigma$ has more than two neighboring cells, i.e. a set of neighboring cells $\mathcal{K}_{\sigma}\subset\Lambda_{h}$. At such bifurcation faces, we enforce flux conservation, just like for faces with exactly two neighbors,

[TABLE]

Inserting Eq. 3.9 into Eq. 3.11 yields

[TABLE]

so that the face unknown in Eq. 3.9 can be replaced by an expression in terms of control volume unknowns $p_{K_{\Lambda}}$.

3.3 Numerical source term $Q_{\Gamma}$

In the following, we consider three different methods: the method proposed in [11] with line sources (ls), the method suggested in [34] with cylinder surface sources (css), and the above introduced method with distributed sources (ds). The methods can be distinguished by the choice of kernel functions (see Section 2), and the approximation of $Q_{\Gamma}$. For the ls and css method, we are looking for the discrete approximation

[TABLE]

where $\partial\mathcal{D}_{R}(s)$ is the boundary of a disc at position $s$, perpendicular to the intersection segment $\Gamma$ with radius $R$. We compute

[TABLE]

where $\mathcal{I}_{\partial\mathcal{D}_{R}}$ is the set of intersections $I_{K}=\partial\mathcal{D}_{R}(s)\cap K_{\Omega}$, $K_{\Omega}\in\Omega_{h}$, and the integral is approximated by a Gaussian quadrature rule. For the method ds we seek

[TABLE]

The center-line pressure $p_{t,0}(s)$ could be approximated by the pressure $p_{K_{\Omega}}$ in the control volume $K_{\Omega}$ containing $\Gamma$. However, this approximation is poor, if $K_{\Omega}$ is not significantly smaller than the vessel radius $R$. Assuming radial flux in a small neighborhood of $\Gamma$, we can reformulate Eq. 3.15 as

[TABLE]

where $p_{t,\delta}$ is evaluated at some distance $\delta\leq 0$ to the vessel center-line, and $\Xi_{\delta}$ is a modified flux scaling factor, ensuring equality. Now, $p_{K_{\Omega}}$ may better approximate $p_{t,\delta}$ than $p_{t,0}$, for a certain $\delta$. As $p_{K_{\Omega}}$ is commonly defined as the mean pressure in the control volume $K_{\Omega}$, we choose $\delta$ as the mean distance of $\Gamma$ to the control volume $K_{\Omega}$,

[TABLE]

The corresponding $\Xi_{\delta}$ is computed using the analytical derivations in Section 2, and is dependent on the kernel functions. Choosing $\Phi^{\text{const}}$, yields

[TABLE]

so that

[TABLE]

Denoting the maximum size of all control volumes $K_{\Omega}$ by $h_{\Omega}$, note that for $h_{\Omega}\rightarrow 0\implies\delta\rightarrow 0$, and $p_{t,\delta}\rightarrow p_{t,0}$, so that this approximation is consistent. For $\delta>\varrho$, a better approximation, independent of the chosen kernel function, can be derived from the analytical solution for $r>\varrho$, Eq. 2.9, yielding

[TABLE]

Note that this $\Xi_{\delta}$ is a sensible flux correction for the methods ls and css, too. However, herein, we do not modify $\hat{q}_{m}$ for these methods and retain the methods as described in [11] and [34].

3.4 Kernel integration

The kernel integral in Eq. 3.3

[TABLE]

is hard to approximate with standard quadrature rules, since $K_{\Omega}\cap\mathcal{S}(\Phi_{\Gamma})$ is difficult to compute. However, the integral over the entire support $\mathcal{S}(\Phi_{\Gamma})$ is known exactly. Hence, the problem can be reformulated as distributing the whole integral over the control volumes $K_{\Omega}$ weighted with the respective support volume fractions. To this end, we create $n_{\Gamma}$ points $\boldsymbol{x}_{i}\in\mathcal{S}(\Phi_{\Gamma})$ with associated volume elements $V_{i}$ of similar size and shape, so that

[TABLE]

4 Numerical experiments

The three considered methods, ls, css, ds, yield different solutions for $p_{t}$ and $r<R$ (ls, css), or $r<\varrho$ (ds). However, the solutions for $p_{v}$ and $\hat{q}_{m}$ are identical. We denote the different pressure solutions in $\Omega$ by $p_{t}^{\mathbb{M}}$, $\mathbb{M}\in\{\textsc{ls},\textsc{css},\textsc{ds}\}$. We analyze these methods with different vessel configurations in a series of numerical experiments. In Section 4.1, we consider a single straight vessel. The numerical methods are investigated in terms of the ratio of grid resolution to vessel radius, comparing with analytical solutions. In Section 4.2, we construct an analytical solution for three parallel vessels, that show that the optimal flux scaling factor $\Xi$ is independent of perturbations caused by neighboring vessels. For each numerical experiment, the setup is described and the results are presented and discussed.

4.1 Single vessel

Let us consider a slightly simplified problem, adapted from [11],

[TABLE]

with the domains $\Omega=[0,1]\times[0,1]\times[0,1]$ and $\Lambda=\{0.5\}\times\{0.5\}\times[0,1]$, i.e the vessel center-line coincides with the $z$ axis. By choosing the parameters as

[TABLE]

the pressure solutions,

[TABLE]

with $r=\sqrt{x^{2}+y^{2}}$, solve Eq. 4.1 given the boundary conditions

[TABLE]

where $\boldsymbol{n}$ is the outward-pointing unit normal vector on the boundary $\partial\Omega$ of the domain $\Omega$. From the analytical pressure solutions follows that $\hat{q}_{m,e}=1+z$ is the analytical source term.

The pressure discretization errors are computed in the normalized discrete norm

[TABLE]

where $p_{K_{\Omega}}$, $p_{K_{\Omega},e}$ denote numerical and exact pressure evaluated at the center of a control volume $K_{\Omega}$ and $|K_{\Omega}|$ its volume. The error for $p_{v}$ in $\Lambda_{h}$ is computed analogously. The error in the source term $\hat{q}_{m}$ is computed as

[TABLE]

where

[TABLE]

The maximum control volume size, $h$, is given by the maximum edge length in both domains. We choose the edge length so that $h=h_{\Omega}$ = $h_{\Lambda}$.

The numerical solutions $p^{\textsc{ds}}_{t}$, $p_{v}$, for $R=0.03$ and $\varrho=R$, are exemplarily shown in Fig. 4. The discretization error and the convergence rates are computed for $p_{t}$, $p_{v}$, and $q$. Fig. 5 shows the discretization errors for $R=0.1$, $\varrho=R$ and $\Phi^{\text{const}}$. It can be seen that ds is the only method with optimal convergence rate for $p_{t}$ in the given norm. The error in $p_{t}$ cannot directly be compared, since it is computed with respect to the respective analytical solution corresponding to the chosen method, which differ for $r<\varrho$. However, all methods are expected to converge to the same analytical solution for $p_{v}$ and $q$. It is evident from the error plots of $p_{v}$ and $q$ that ds shows the lowest discretization error of all three methods. Furthermore, for the presented numerical experiment ds achieves convergence rates in $q$ of approx. $2.5$, while $\textsc{ls},\textsc{css}$ show convergence rates of approx. $2$. We conclude that the increased discretization error for ls and css in $p_{t}$, influences the approximation of $p_{v}$, $q$, and is likely due to the insufficient approximation of $p_{t,\mathbb{W}}$.

Fig. 6 shows the discretization errors for $R=0.05$, $\varrho=R$ and $\Phi^{\text{const}}$. The first two error measurements show the situation $h>\varrho$. If the discretization length $h$ is larger than the kernel radius $\varrho$, the three methods do not differ in the representation of the source term. However, for the ds method, we introduced a source term correction by an adjusted flux scaling factor dependent on the discrete distance $\delta$ (see Eq. 3.18), which is proportional to $h$. It can be seen, that this adjustment significantly reduces the error for all quantities. Note that this scaling factor for $r>R$, can also be applied for the other methods (ls and css), however its motivation is directly drawn from the new formulation of the perfusion problem in Eq. 2.1.

4.2 Multiple parallel vessels

In the next experiment, we consider three parallel vessels, from which two are emitting fluid (arterial side), and one vessel is absorbing fluid (venous side). We assume constant pressures in the vessels, $\tilde{p}_{v,1}=\tilde{p}_{v,3}=$3400\text{,}\mathrm{Pa}, $\tilde{p}_{v,2}=$2300\text{\,}\mathrm{Pa}, which renders the problem effectively two-dimensional. Then, exploiting the linearity of the Laplace operator, we can construct a solution of Eq. 1.2 in $\Omega$ using the superposition principle, i.e.,

[TABLE]

where $\boldsymbol{x}_{i}$ is the position of the center-line of vessel $i$, and $\boldsymbol{x}=[x,y,z]^{T}$ a point in $\Omega$. Recall that $\hat{q}_{m,i}=\rho_{I}L_{p}2\pi R_{i}(p_{v,i}-p_{t,\mathbb{W}_{i}})$ is a linear function in the arguments $p_{v,i}$ and $p_{t,\mathbb{W}_{i}}$, and

[TABLE]

is the mean pressure on the surface of vessel $i$. Thus, using the mean value property of harmonic functions yields

[TABLE]

where $r_{ij}$ denotes the Euclidean distance of the center-lines of vessels $i$ and $j$ and $R_{i}$ the radius of vessel $i$. Eq. 4.9 constitutes a system of three linear equations with the unknowns $p_{t,\mathbb{W}_{1}}$, $p_{t,\mathbb{W}_{2}}$, $p_{t,\mathbb{W}_{3}}$. As the analytical solution of Eq. 4.9 results in a rather lengthy expression, we compute $p_{t,\mathbb{W}_{i}}$ numerically. Using the parameter values, $r_{12}=r_{23}=$40\text{,}\mathrm{\SIUnitSymbolMicro m}, $r_{12}=$80\text{\,}\mathrm{\SIUnitSymbolMicro m}, $R_{1}=R_{2}=R_{3}=$4\text{,}\mathrm{\SIUnitSymbolMicro m}, $L_{p}=$1\text{\times}{10}^{-9}\text{\,}\mathrm{m}\text{\,}{\mathrm{Pa}}^{-1}\text{\,}\mathrm{s}, $\kappa=$8.3\text{\times}{10}^{-18}\text{,}{\mathrm{m}}^{2}, $\mu_{I}=$1.339\text{\times}{10}^{-3}\text{\,}\mathrm{Pa}\text{\,}\mathrm{s}, $\rho_{I}=$1030\text{,}\mathrm{kg}\text{,}{\mathrm{m}}^{-3}, $\Pi_{t}=$666\text{\,}\mathrm{Pa}, $\Pi_{v}=$3300\text{,}\mathrm{Pa}, $\sigma=1.0$, cf. [[4](#bib.bib4)], we obtain $\tilde{p}_{t,\mathbb{W}_{1}}=\tilde{p}_{t,\mathbb{W}_{3}}\approx$587.29\text{\,}\mathrm{Pa}, and $\tilde{p}_{t,\mathbb{W}_{2}}\approx$-89.350\text{,}\mathrm{Pa}$$. Recall that $p_{t/v}=\tilde{p}_{t/v}-\Pi_{t/v}$.

Again, the solution with line sources is equal to the solution with distributed sources outside the kernel support radius. Furthermore, the superposition principle equally applies for the distributed source model so that, for instance, for the kernel function $\Phi^{\text{const}}$,

[TABLE]

is a solution to Eq. 2.1, given that the kernel support regions of the vessels do not overlap. For the css method the solution is given by a linear continuation of the pressure for $r_{i}\leq R_{i}$, cf. [34],

[TABLE]

The analytical solutions $p^{\textsc{ls}}_{t}(\boldsymbol{x})$, $p^{\textsc{css}}_{t}(\boldsymbol{x})$, $p^{\textsc{ds}}_{t}(\boldsymbol{x})$ for a domain $\Omega=[-100,100]\times[-50,50]$\mathrm{\SIUnitSymbolMicro m} with $\boldsymbol{x}_{1}=[-40,0]^{T}$\mathrm{\SIUnitSymbolMicro m}, $\boldsymbol{x}_{2}=[0,0]^{T}$\mathrm{\SIUnitSymbolMicro m}, $\boldsymbol{x}_{3}=[40,0]^{T}$\mathrm{\SIUnitSymbolMicro m} along the x-axis are shown in Fig. 8. The numerical solution $\tilde{p}^{\textsc{ds}}_{t}=p^{\textsc{ds}}_{t}-\pi_{t}$, for $\varrho=2R$, is exemplarily shown in Fig. 7.

The discretization errors with respect to the analytical solutions for $p^{\mathbb{M}}_{t}$ and $q$ are shown in Fig. 9. As for the numerical experiment with a single vessel, method ds shows the lowest error for the multi-vessel experiment. Additionally, it can be seen that the flux-scaling for $r>\varrho$ significantly improves the approximation of $q$, although the analytical solution is no longer strictly radial around the individual vessels.

In a second experiment, the grid resolution is fixed to $h=$1.25\text{,}\mathrm{\SIUnitSymbolMicro m}. Then, the kernel support radius is step-wise increased starting from $\varrho=R$, for all vessels. The discretization errors for $q$ are shown in [Fig. 10](#S4.F10). Recall that the distances $r_{12}=r_{23}=$40\text{\,}\mathrm{\SIUnitSymbolMicro m}, such that the kernel support region for two neighboring vessels start intersecting for $\varrho>5R$, and the kernel support region includes the location of the center-line of the neighboring vessel for $\varrho>10R$.

It can be seen that with increasing kernel radius the approximation of the source term improves significantly. Comparing Fig. 10 and Fig. 9 it seems that an increase of the kernel support radius has the same effect as increasing the grid resolution. An increase of the kernel support area by a factor of $2$ reduces the discretization error by a factor of $2$. However, note that the pressure solution is increasingly regularized and thus deviates from the physically sensible solution. If the kernel support radius is chosen too large, the regularization affects neighboring vessels, such that the error in the source term increases again. The results suggest that the variable kernel support region, to some extent, decouples the source term approximation error from the 3D grid resolution.

4.3 Multiple arbitrarily-oriented vessel

In this experiment, we consider multiple arbitrarily-oriented vessels embedded in a cubic extra-vascular domain, for which no analytical solution is given. Recall that for such a system, the methods ls, css and ds do not generally give the same solutions, but the differences are expected to be small. In order to isolate differences stemming from the different source models from other sources of error, the vessels do not bifurcate and are constructed in a way that they intersect the 3D domain boundary perpendicularly. The computational domain is shown in Fig. 11. Two arterial (fluid-emitting) and two venous (fluid absorbing) capillaries are embedded in a cubic extra-vascular domain with dimensions $100\text{\times}100\text{\times}100\text{\,}\mathrm{\SIUnitSymbolMicro m}$.

For each vessel, one end is chosen as inflow and the other as outflow boundaries. On the inflow boundary of vessel $i$, a constant inflow rate, $q_{v,i,\text{in}}=\rho_{B}\pi R_{i}^{2}v_{i}$ is enforced, where $R_{i}$ is the inflow segment radius and $v_{i}$ the inflow velocity. On outflow boundaries, we fix the pressure, $\tilde{p}_{v,i,\text{out}}$. The data for the geometry of the vessel and boundary conditions is given in Table 1 in the appendix. The remaining model parameters are chosen as in the previous multiple vessel experiment. The boundaries of the extra-vascular domain are considered symmetry boundaries, hence $\frac{\partial p_{t}}{\partial\boldsymbol{n}}=0$ on $\partial\Omega$.

We produce a reference solution using the css method, with $h_{\Omega}=$0.625\text{,}\mathrm{\SIUnitSymbolMicro m}, $h_{\Lambda}=$0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}. The source term $\hat{q}_{m}$ is computed for every cell $K\in\Lambda_{h}$, resulting in a source vector $q_{\text{ref}}$. Then, the corresponding source term, $q$, for $h_{\Lambda}=$0.5\text{,}\mathrm{\SIUnitSymbolMicro m}$$, is computed for different $h_{\Omega}$ using the methods css, ls, and ds with different kernel radii $\varrho/R_{i}=1,3$, and $5$. Furthermore, we compute the total mass flux, $q_{\text{out}}$, emitted by the arterial vessels, as the sum of all fluxes leaving the vessel domain into the extra-vascular domain. The results are shown in Fig. 12.

First, it can be seen that with a coarse grid, the difference in the source term computed by the different schemes is quite large ($\approx$10\text{,}\mathrm{\char 37\relax}$$ relative to the reference solution). Notably, the lowest difference with respect to the reference solution at coarse resolutions is achieved by the ds with the largest kernel support region. With grid refinement these differences decrease to less than $1\text{\,}\mathrm{\char 37\relax}$. However, it can be seen that the curve flattens for the ds methods with larger kernel for fine grid resolutions. The results suggest that the difference cannot be reduced to much less than $1\text{\,}\mathrm{\char 37\relax}$. At the given grid resolutions, such a behavior is not observed yet for the ds with the smallest kernel support.

We conclude that the flattening of the error curve for the ds methods with larger kernel is rather caused by overlapping kernels of neighboring vessels as well as at bends. This effect, which is further investigated in the next numerical example, is minimized for the smallest kernel support. From the fact, that the kernel curve still suggests convergence until a very fine grid resolution leads us to the conclusion that the error caused by the approximation of the flux scaling factor $\Xi$ (as discussed in Section 2.3) is much smaller. In perspective of the rather big uncertainties stemming from vessel segmentation and modeling error, the difference of $1\text{\,}\mathrm{\char 37\relax}$ between the different methods present in Fig. 12 will most certainly be acceptable in practical simulations. Furthermore, the observations in this experiment support our results from the previous experiment that for coarser mesh resolutions it is better to choose a larger kernel support, if a good approximation of the source term $\hat{q}_{m}$ is important.

4.4 Vessel network

In the last numerical experiment, we consider a network of capillaries extracted from the superficial cortex of a rat brain [37, 18]. Inlets and outlets are annotated in the data set. For the inlets, velocity estimates based on the vessel radius are given in [18], and herein enforced as Neumann boundary conditions. The vessel radii are in the range of $24.5\text{\,}\mathrm{\SIUnitSymbolMicro m}$. Dirichlet boundary conditions enforce $p_{v,\text{out}}=$1.025\text{\times}{10}^{5}\text{,}\mathrm{Pa}at the outlets. The extra-vascular domain $\Omega$ is given by a rectangular box,200\text{,}\mathrm{\SIUnitSymbolMicro m}$\times$210\text{,}\mathrm{\SIUnitSymbolMicro m}$\times$190\text{,}\mathrm{\SIUnitSymbolMicro m}$$. All boundaries $\partial\Omega$ are considered symmetry boundaries, $\frac{\partial p_{t}}{\partial\boldsymbol{n}}=0$ on $\partial\Omega$. The network boundaries are extended by $30\text{\,}\mathrm{\SIUnitSymbolMicro m}$ segments with perpendicular intersections on $\partial\Omega$. This adjustment to the network structure is necessary to better match the assumption of symmetry boundaries in the extra-vascular domain. Vessels intersecting the boundary at acute angles lead to non-physical, very non-radial flows around the vessel end at symmetry boundaries on $\partial\Omega$.

A reference solution is computed using the css method with $h_{\Omega}=$1.3125\text{,}\mathrm{\SIUnitSymbolMicro m} and $h_{\Lambda}=$0.5\text{\,}\mathrm{\SIUnitSymbolMicro m}. The network geometry and the $p_{v}$ reference solution are shown in Fig. 13. The source term, $q$, for $h_{\Lambda}=$0.5\text{,}\mathrm{\SIUnitSymbolMicro m}$$, is computed for different $h_{\Omega}$ using the methods css, ls, and ds with different kernel radii $\varrho/R_{i}=1,3$, and $5$. As in the previous example, we compute the total mass flux, $q_{\text{out}}$. The results are shown in Fig. 14.

Firstly, it can be seen that all methods agree well for the total mass flux on the finest grid with maximum differences of about $2\text{\,}\mathrm{\char 37\relax}$. Secondly, it is evident that the ds approximates $q_{\text{out}}$ much better for coarser grids, and the approximation gets better with larger kernel support. The ds method with a kernel support radius $\varrho_{i}=5R_{i}$ very closely approximates $q_{\text{out}}$ for the coarsest grid, where $h_{\Omega}$ is about five times the radius of the largest vessel. If the difference is measured in the relative $2$-norm, $\Delta q=||q-q_{\text{ref}}||_{2}/||q_{\text{ref}}||_{2}$, we observe slightly larger differences. On the one hand, this is due to differences along vessels oscillating around zero, i.e., due to terms that cancel when computing $q_{\text{out}}$, but not for $\Delta q$. On the other hand, $\Delta q$ emphasizes larger differences stronger.

The locality of the difference is visualized in Fig. 15, where we computed the absolute local differences of source terms $\hat{q}_{m,i}$ between the ls and ds methods ($\varrho/R_{i}=1$) and the reference solution (with css) for every vessel cell $K_{\Lambda}$.

The largest differences can be observed in the neighborhood of bifurcations and sharp bends. This difference can be explained by the fact that the vessels are discretely represented as cylindrical tubes and no explicit intersection geometry is assumed for such features. To the best of our knowledge, there is currently no mixed-dimension embedded method available that considers particular bifurcation geometries. For the ls method this means that at such features the mean surface pressure $p_{t,\mathbb{W}}$ contains contributions from the non-physical part of the pressure solution ($r<R$) of the neighboring vessels, leading to an overestimation of $p_{t,\mathbb{W}}$. This defect is improved in the css method, where, due to the linear continuation of the solution for $r<R$, the value of $p_{t,\mathbb{W}}$ is closer to the actual pressure on the vessel wall surface. For the ds method, the approximation of $p_{t,\mathbb{W}}$ is based on an analytical reconstruction from the center-line pressure $p_{t,0}$. However, at bifurcations the kernel support of the neighboring vessels overlap, which leads to imprecision in the approximation of $\Xi$. Fig. 14 shows that these effects lead to differences of $34\text{\,}\mathrm{\char 37\relax}$ in $\Delta q$ between the methods. Furthermore, it is shown that if the source contributions in the vicinity of bifurcations ($K_{\Lambda}$ closer than $10\text{\,}\mathrm{\SIUnitSymbolMicro m}$) are excluded from the norm, the differences decrease to $23\text{\,}\mathrm{\char 37\relax}$. This shows that the differences are rather local to the bifurcation neighborhood. Finally, while it cannot be concluded which method is best for the fine grid solutions, we again observe that the ds method, especially the variants with larger kernel support, better approximate the fine scale solution for large $h_{\Omega}$ (coarse grid). In fact, the difference to the fine scale solution (obtained with the css method) for the coarsest grid ($h_{\Omega}=$21\text{,}\mathrm{\SIUnitSymbolMicro m}$$) is $10\text{\,}\mathrm{\char 37\relax}$ for the ds method with $\varrho_{i}=5R_{i}$, while it is $30\text{\,}\mathrm{\char 37\relax}$ for the ls and css methods.

5 Summary and conclusion

We presented a new method for modeling tissue perfusion using a mixed-dimension embedded method with distributed sources. The most prominent difference to existing schemes is that the source term, coupling the vascular and the extra-vascular domains, is spatially distributed using kernel functions. The mean pressure on the vessel surface is not explicitly computed but locally reconstructed from the pressure at the vessel center-line using an analytically derived scaling factor. We showed in four numerical experiment that the result obtained with the new method match well with the results obtained with existing methods. It was consistently shown in the experiments that the new method converges with a higher rate for the source term $\hat{q}_{m}$ and the vessel pressure $p_{v}$ due to an increased regularity in the extra-vascular pressure $p_{t}$. Furthermore, in all experiments the new method provided better approximations of the source term $\hat{q}_{m}$ for coarse extra-vascular grids. In a test with a realistic vessel network from the rat cortex, using the new method resulted in a three-fold reduction of the error in the source term for coarse grids in comparison to state-of-the-art methods with respect to a reference solution computed on a fine grid.

In four numerical experiments the new method was compared to existing methods. To this end, all methods have been implemented in the open source software framework DuMuX [38, 21]. The implementation effort and computational costs (at the same grid resolution) of the new method are comparable with existing methods. However, our results suggest that the new methods can be considered computationally more efficient, since a good approximation of the source term is already achieved at lower grid resolutions.

The results in this paper suggest that the new method provides the best source term approximations, if the source is distributed over a volume larger than the vessel itself. In such a case the pressure solution in the extra-vascular space is regularized and thus may deviate from the physical solution. In turn, the regularized extra-vascular pressure results in better approximations of source term and vessel pressure. A good approximation of the fluid exchange between vascular and extra-vascular compartments is crucial in many applications involving transmural transport processes such as the estimation of contrast agent leakage from the brain microvasculature in multiple sclerosis [4]. A regularized extra-vascular pressure may be acceptable in cases where it is more interesting how much of a substance leaves the vascular system rather than its accurate distribution in the extra-vascular space. Furthermore, it is always acceptable to choose the distribution volume similar to the size of the coarsest neighboring extra-vascular discretization cell. In such a case, a possible error in the extra-vascular pressure is masked by other discretization errors. Nevertheless, the description of the source term and thus the vessel pressure may still be improved significantly by a better reconstruction of the mean vessel surface pressure.

The last numerical experiment revealed that the largest differences between the different mixed-dimension embedded methods occur in the vicinity of bifurcations. This is due to an imprecise description of bifurcation geometries in the discrete setting. In order to evaluate methods with respect to the error at bifurcations, and to develop improved descriptions of the flow around bifurcations, a comparison with methods with spatially resolved interfaces is necessary. However, our results suggest that discretization errors around bifurcations only affect small parts of the entire system, such that the current models might be sufficient approximations for most applications.

Finally, due to the improved accuracy at coarse grid resolutions, we consider the new method an important step towards simulations of larger vessel networks, where fine grid resolutions in the extra-vascular space may get computationally prohibitively expensive.

Acknowledgements

This work was financially supported by the German Research Foundation (DGF), within the Cluster of Excellence in Simulation Technology (EXC 310), and the Collaborative Research Center on Interface-Driven Multi-Field Processes in Porous Media (SFB 1313, Project Number 327154368).

Appendix

A1 Domain information for the numerical experiment in Section 4.3

The vessel network used in Section 4.3 is described in Table 1.