A 3D-1D coupled blood flow and oxygen transport model to generate microvascular networks
AA 3D-1D coupled blood flow and oxygentransport model to generate microvascularnetworks
T. K¨oppl , E. Vidotto , B. Wohlmuth Key words: blood flow simulations, oxygen transport,dimensionally reduced models, vascular growth, flows inporous media, 3D-1D coupled flow modelsMathematics Subject Classification:
Abstract
In this work, we introduce an algorithmic approach to generate mi-crovascular networks starting from larger vessels that can be reconstructedwithout noticeable segmentation errors. Contrary to larger vessels, thereconstruction of fine-scale components of microvascular networks showssignificant segmentation errors, and an accurate mapping is time and costintense. Thus there is a need for fast and reliable reconstruction algo-rithms yielding surrogate networks having similar stochastic propertiesas the original ones. The microvascular networks are constructed in amarching way by adding vessels to the outlets of the vascular tree fromthe previous step. To optimise the structure of the vascular trees, we useMurray’s law to determine the radii of the vessels and bifurcation angles.In each step, we compute the local gradient of the partial pressure of oxy-gen and adapt the orientation of the new vessels to this gradient. At thesame time, we use the partial pressure of oxygen to check whether theconsidered tissue block is supplied sufficiently with oxygen. Computingthe partial pressure of oxygen, we use a 3D-1D coupled model for bloodflow and oxygen transport. To decrease the complexity of a fully coupled3D model, we reduce the blood vessel network to a 1D graph structureand use a bi-directional coupling with the tissue which is described by a3D homogeneous porous medium. The resulting surrogate networks areanalysed with respect to morphological and physiological aspects.
The reconstruction of vascular networks from medical images plays an impor-tant role in different research areas and has several application fields. Modellingand simulation of blood flow and transport processes within vascular networksis of great interest, since it provides the possiblity to study the distribution of Department of Mathematics, University of Technology Munich, Boltzmannstr. 3,85748 Garching bei M¨unchen, Germany, [email protected], [email protected],[email protected] a r X i v : . [ c s . C E ] J u l njected medications [6, 9]. Moreover, a reliable and robust blood flow modelallows one to investigate the impact of diseases and therapeutic procedures in anon-invasive way. Numerical simulation techniques have become an importantpart in testing hypotheses when experimental investigations are not possible orare limited by accessibility, size and resolution [10, 11, 24, 22, 34, 40, 56]. To beable to perform conclusive blood flow simulations, it is crucial to have a precisedescription of the network. This means that accurate data on the radii, lengthsand connectivity of the vessels as well as the location of the vessels are required.In order to obtain high-quality data, remarkable progress using high-resolutionimaging techniques has been made. As a result, the acquisition of suitable im-age data with respect to complex vascular systems on different scales has beenfacilitated [47]. Nevertheless, an accurate reconstruction of vascular networks isstill challenging, in particular, if fine-scale microvascular networks are consid-ered. Most algorithms are concentrated on vessel detection and segmentation,while the determination of the vascular connectivity is neglected [49]. Furtherproblems are associated with the image recording (image noise and artifacts)and sample preparation (air bubbles and clotting) [27].Motivated by such problems, different approaches for generating microvas-cular networks have been developed [8]. One approach considers the theory ofporous media. Thereby, the fine scaled vessels are homogenised and consideredas a porous media. The pores and pore throats within these porous media aregiven by the blood vessels and its branchings [59]. In order to account for thevessel hierarchy multi-compartment flow models have been considered [20, 50].A challenge in this context is to estimate the permeability tensors and porosi-ties for each compartment and to derive suitable coupling conditions betweenthe different compartments. Other types of models directly simulate vasculargrowth taking optimality principles like the minimisaton of building materialor the minimisation of energy dissipation into account [29, 35, 44]. Insteadof simulating the growth of the vessels, a direct optimisation of the microvas-cular system can be considered. In this case, the optimisation processes areconcentrated on intravascular volume minimisation with constraints based onphysiological principles [17].In this work, we adapt concepts developed in [53, 54] to generate artificialarterial trees. Thereby the authors use ideas from modelling approaches ofangiogenesis processes [41, 42, 55]. One of these concepts is based on the as-sumption that the arterial tree generation is driven by the oxygen distributionin the tissue. After measuring the oxygen concentrations, the intensity of vascu-lar endothelial growth factors (VEGFs) triggering the sprouting of new vesselsis computed. Concerning the choice of the vessel radii and bifurcation angles,minimisation principles developed by Murray are used [32, 33]. Furthermore,other algorithms may be applied to optimise the vascular tree e.g. by removingtiny and redundant vessels. The arterial tree is constructed stepwise. In eachstep, firstly, the oxygen concentration in the tissue is computed. According tothe new distribution of oxygen and Murray’s princples, the directions of growthfor new vessels are determined. Redundant vessels are removed before obtainingthe final network structure.For our modelling approach, we adapt some of the ideas presented in [54] tosimulate the growth of an artificial network. However, we are not only focus-ing on the generation of an arterial tree. In a first step, we consider arteriolesand venoles that are well reconstructed from given image data. Based on these2essels, the missing part of the microvascular network linking the larger vesselsis generated. In this context, we do not restrict ourselves to the arterial com-ponent, but our starting network can comprise both venous and arterial com-ponents. A further enhancement is related to the computation of the oxygenconcentration in a tissue block that is supplied by the microvascular network.According to [54, Appendix A], the oxygen concentration is approximated by asuperposition of heuristically chosen functions. In the submitted work, we con-sider 3D-1D coupled models for simulating blood flow and transport of oxygenin microvascular networks. The term 3D-1D indicates that the microvascularnetwork is described by a one-dimensional (1D) graph-like structure while thesurrounding tissue is considered as a three-dimensional (3D) porous medium.This implies that blood flow and the transport of oxygen within the microvas-cular are governed by 1D PDEs, while flow and oxygen transport in the tissueare modelled by Darcy’s law [18] and a standard convection-diffusion equa-tion. The convection-diffusion equation is enhanced by the Michaelis-Mentenlaw modelling the consumption of oxygen by the tissue cells [6, Chapter 2]. Tocouple the PDE systems for flow and transport, we use a specific coupling con-cept presented in [5, 25, 26]. The key ingredient of this concept is to couple thesource terms of the PDEs for flow and transport. Thereby, Starling’s filtrationlaw or the Kedem-Katchalsky equation [4] are used to determine the exchangeof fluids and oxygen. In order to embed the 1D PDEs into the source termsof the 3D PDEs, we use Dirac measures concentrated on the vessel surfaces,instead of other approaches, such as [7]. This is motivated by the fact that theactual exchange processes between the vascular system and the tissue take placeat these locations.The remainder of this work is structured as follows: Section 2 is devoted tothe modelling assumptions, the mathematical equations for blood flow and oxy-gen transport, as well as the generation of the blood vessel network. Moreover,we discuss the numerical methods that have been used to solve the mathemati-cal model equations. In Section 3, we show numerical results and report on theinfluence of typical parameters as well as statistical characterisations. Finally,in Section 4, we summarise our main insights and the key ideas developed inthis paper. The main section of this paper is divided into six subsections. In the first sub-section, we present the data set for a microvascular network that is used forour numerical simulations in Section 3. Next, the basic modelling assumptionsthat are used in the following subsections are listed and motivated. To simulateblood flow and transport processes within a microvascular network, we developin the third and fourth subsection dimension-reduced models for flow and trans-port in the network and the surrounding tissue. Thereby, flow and transport inthe vascular system are governed by 1D models, while for flow and transport intissue 3D porous media models are considered, as noted earlier. Concerning thetransport processes, we restrict ourselves to the transport of oxygen. This ismotivated by the fact that a suitable microvascular network has to prevent anunbalanced distribution of the partial pressure of oxygen (
P O ) within the tis-sue block under consideration. Based on the P O distribution in the tissue, we3evelop in the next subsection a model for artificial vascular growth to producea network able to supply the tissue block with a sufficient amount of oxygen.In the final subsection, we briefly describe the numerical methods, used to solvethe equations governing the mathematical model.For our numerical simulations, we use the microvascular network shown inFigure 1 (left). This microvascular network is part of a larger vascular networktaken from the cortex of a rat brain. The whole vascular network is describedin [48]. According to this reference, its vessel midlines together with the 3Dcoordinates of branches, kinks as well as inlets and outlets are reconstructedfrom X-ray images. In the remainder of this paper, we denote these points asnetwork nodes. Curved midlines are approximated by straight edges and foreach vessel a mean radius is estimated. To each network node that is locatedwithin the sample, a fixed pressure value is assigned. The larger vessels of thisnetwork considered in this work are contained in a domain Ω roi ⊂ R given by(roi: region of interest):Ω roi = (cid:8) x ∈ R (cid:12)(cid:12) L r ,l < x < L r ,u , L r ,l < x < L r ,u , L r ,l < x < L r ,u (cid:9) . However, we consider a larger cuboid Ω , Ω roi ⊂ Ω ⊂ R as our tissue domain:Ω = (cid:110) x ∈ R | L ,l < x < L ,u , L ,l < x < L ,u , L ,l < x < L ,u (cid:111) . Thereby, we choose Ω such that Ω roi ⊂ Ω. This is motivated by the fact thatif we consider Ω roi as a tissue domain, the resulting artificial network might beaffected by the chosen boundary conditions. In order to determine the largerdomain Ω, each edge of Ω roi is enlarged on both ends by about 10 %. The valuesfor the lengths defining Ω roi are provided in Table 2. We assume further thatthe given microvascular network can be represented by a union of N cylinders V k with a radius R k and a straight center line or edge E k , k ∈ { , . . . , N } , seeStep 2 in Figure 2. By this, the domain Ω can be decomposed in two parts:Ω v = N (cid:91) i =1 V k and Ω t = Ω \ Ω v , where Ω v stands for the vascular system and Ω t the tissue. We denote the twoendpoints of Λ k by x k, , x k, ∈ Ω and the length of V k by l k = (cid:107) x k, − x k, (cid:107) ;the orientation λ k of the cylinder V k is defined by the following expression: λ k = x k, − x k, l k . Accordingly, the edge E k of V k and is given by: E k = { x ∈ Ω v | x = x k, + s · λ k , s ∈ (0 , l k ) } , where s denotes the arc length parameter. Γ k is the lateral surface of the cylin-der V k . The union Γ = (cid:83) Nk =1 Γ k provides an approximation of the surface of thevascular system, due to the fact that the cylinders V k may not perfectly matchto each other. Furthermore, the set Λ = (cid:83) Nk =1 E k denotes the 1D network com-posed of the different edges. In the next step, we extract the larger vessels fromthe microvascular network Ω v . For this purpose, we choose a fixed threshold R T . × . × . lv = (cid:83) Nk =1 { V k ⊂ Ω v | R k > R T } andΛ lv = (cid:83) Nk =1 { E k ⊂ Λ | R k > R T } . The motivation for defining this set is thatsuch vessels may be obtained by means of imaging techniques with acceptableaccuracy. Based on Ω lv and Λ lv , we construct the surrogate networks. Before presenting the modelling equations for blood flow, oxygen transport andvascular growth, we discuss the most important modelling assumptions and sim-plifications, which are fundamental for the choice of the corresponding models.In total, nine assumptions for our modelling concept are formulated:(A1)
Gravity is neglected.
Since the tissue block and the vascular system comprise a volume of about1 . , we can expect that only little fluid mass is contained in thissystem. As a result, the gravity force caused by the weight of the bloodvolume in the upper part of this system has almost no influence on theflow in the lower part.(A2) Blood is an incompressible fluid, i.e., its density is assumed tobe constant.
Assumption (A2) does not hold in general, due to the fact that the redblood cells in a microvascular network are not distributed in a homoge-neous way [36, 43, 52]. Besides the red cells, blood is composed of plasma,white cells and platelets [14, Chap. 1], which implies that on the level ofmicrocirculation blood can not be considered as a homogeneous fluid witha constant density. Despite of this fact, in many publications that are5 . Outline of a blood vessel network 2. Approximation of the blood vessels by cylinders with a constant radius network nodes centerlines3. Decomposition of the blood vessel network, modelling the exchange with tissue and derivation of 1D fl ow models for each vesselexchange with tissue inner network nodes4. Coupling of the 1D fl ow models at the inner network nodes and formulation of boundary conditions inner boundary nodesouter boundary nodes Figure 2: This figure illustrates the four basic steps for modelling of blood flowand oxygen transport in microvascular networks. In a first step, we consider theoutline of a microvascular network Ω v . Next the blood vessels are approximatedby straight cylinders. After that, the network is decomposed, and 1D models forblood flow and oxygen transport for each vessel are derived. In the last step, the1D models are coupled and boundary conditions are derived, such that bloodflow and oxygen transport can be simulated within the whole 1D network.concerned with modelling of blood flow in microvascular networks, bloodis considered as an incompressible fluid [4, 12, 38, 40, 55].(A3) The non-Newtonian flow behavior is accounted for by an alge-braic relationship.
The red blood cells govern the viscosity of blood, significantly. As the redblood cells have to deform such that they can move through capillaries,the viscosity varies within the microvascular network. A quantitative re-lationship between the vessel diameter d is given by the following formulafor the in vivio viscosity µ bl [Pa · s] [45]: µ bl ( d ) = µ p (cid:32) µ . −
1) (1 − H ) C − − . C − · (cid:18) dd − . (cid:19) (cid:33) · (cid:18) dd − . (cid:19) . (1)In (1), the diameter d is dimensionless. The physical diameter [ µ m] has6o be divided by 1 . µ m to obtain d . µ p [Pa · s] denotes the viscosity ofblood plasma, and H stands for the discharge hematocrit which is definedby the ratio between the volume of the red blood cells and the total bloodvolume. The apparent viscosity µ . is given by: µ . = 6 . − . · d ) + 3 . − .
44 exp (cid:0) − . · d . (cid:1) , and C is a coefficient determining the influence of H on µ bl : C = (0 . − . · d )) (cid:18) − − d (cid:19) + 11 + 10 − d . In this context, one should be aware of the fact that the constitutiverelationship (1) is known to hold for human blood. For rat blood, we donot have a suitable function at hand. However, we apply this functionalso to rat blood.(A4)
Inertial effects concerning flows in the vascular system and tissueare not considered.
According to [14, Tab. 1.7], blood velocity is about 0 . / s in thearterioles and venules and about 0 .
01 mm / s in the capillary bed in ahuman system. Therefore, it can be concluded that the Reynolds numbersare significantly lower than 1 . Pulsatility of blood flow is neglected.
The pulsatility of blood flow can be neglected, since the Womersley num-bers in the microcirculation are much smaller than 0 .
1. The Womersleynumber is given by a dimensionless number comparing the frequency of apulse and the viscosity of a fluid to each other [14, Tab. 1.7].(A6)
Tissue surrounding the microvascular network is considered asan isotropic, homogeneous porous medium.
Considering the tissue, it can be observed that it is mainly composed ofcells, fibers and the interstital space filled with a fluid similar to bloodplasma. The interstital space exhibits several pores that are connectedby pore throats. Therefore, it is reasonable to consider tissue as a porousmedium. In order to model flow within the tissue, we do not use porenetwork models [61] that attempt to resolve each pore and pore throat, butwe apply homogenised or REV-based flow models such as Darcy’s law [21,40, 58]. In this context, it is assumed that the parameters characterisingthe porous structure like porosity and permeability are constant. Since wedo not have any information on the tissue surrounding the given network,we assume that the porosity is constant and that the permeability tensoris diagonal and isotropic with constant values. However, this assumptiondoes not hold for every tissue type. In muscle tissues, for example, themuscle fibers form small channels in which the interstitial fluid can flowfaster than in other space directions. This means that the permeabilitytensor is not isotropic in this case. If larger domains of an organ like theheart or the liver have to be considered the permeability and porosity ofthe different tissue layers my vary, due to the fact that different partsof an organ have different functions. A possible way to incorporate thisfeature into a computational model is to consider a multi-compartment7odel to simulate flow in tissue composed of multiple layers having adifferent permeability or porosity. In [20, 50] this modelling approach hasbeen used to simulate vascular systems that are partially homogenised andrepresented as a porous medium. However, this technique could also beadapted to a heterogeneous tissue.(A7)
No lymphatic drainage is considered.
This is due to the fact that we consider a portion of brain tissue. Ac-cording to medical literature, brain tissue exhibits no lymphatic systemas other organs. Over time, several hypotheses have been made on theway interstitial fluid is drained from the brain tissue. Most researchersthink that cerebrospinal flow is mainly responsible for the drainage tissue,but the exact exit route of cerebrospinal fluid is still subject of controver-sies [1, 28]. Therefore, we do not incorporate lymphatic drainage into themodel that is presented in this work. Some numerical models and furtherinformation on lymphatic systems in different organs can be found e.g. in[30, 46, 57].(A8)
Distribution of oxygen in both the vascular system and tissue isrepresented by a continuous unit.
In our case, we use the partial pressure of oxygen
P O . The propagationof oxygen is governed by diffusion and advection in the vascular systemand tissue. In order to model the transport of oxygen in blood, one hasto take into account that oxygen molecules are attached to the red bloodcells. Therefore, a particle tracer model could be used to determine theoxygen concentration within the vascular network. In this work, we useboth for the blood vessel network and the tissue a continuous unit. Acommonly used unit to describe the oxygen distribution is the partialpressure ( P O ), which is related to the oxygen concentration by Henry’slaw [19, 60]. The consumption of oxygen by the tissue cells is modelledby means of the Michaelis-Menten-law [6, Chapter 2][24] containing anaverage consumption rate and a mean value of the partial pressure inbrain tissue.(A9) New vessels are only added at the boundary nodes of the termi-nal vessels that are not adjacent to the boundary of the tissuedomain.
Finally, we assume that the creation of new vessels occurs only at the out-lets of vessels that are contained in the inner part of the tissue domain.Since our network is formed out of the arterioles and venules which can bedetected in the data set, one can imagine that the missing capillary bedcan be created by branches or the growing of single vessels at the outletsof the terminal vessels. In medical terms, the formation of new branchesat the outlets of blood vessels is referred to as apical growth, while theemergence of new vessels at other places is denoted as vascular growth bysprouting [54]. 8 .2 3D-1D modelling of blood flow in microvascular net-works and tissue
Let us consider the network Ω v . As a first step towards a mathematical modelfor flow in the vascular network, we decompose Ω v into its single vessels V k .Then for each vessel V k , a 1D flow model is derived (see Step 3 in Figure 2).In order to obtain a flow model for the whole network flow, coupling conditionsconnecting the single vessels as well as boundary conditions are specified (seeStep 4 in Figure 2). Finally, a model for tissue flow is determined.Taking the modelling assumptions (A4) and (A6) from the previous sub-section into account, the inertia terms in the standard fluid equations may bedropped. Due to (A5), there is no need to consider complex fluid structureinteraction models relating blood flow and deformations of the vessel walls toeach other. All in all, it is admissible to simulate both tissue and vascularflow by means of Darcy type equations combined with conservation equationsfor incompressible fluids. According to [62], the permeability K V k for the vessel V k ⊂ Ω v is given by: K V k = ( R k ) / . Using (1), the blood viscosity is computedas: µ bl ( x ) = µ bl (cid:18) R k . µm (cid:19) =: µ ( k )bl for x ∈ V k ⊂ Ω v . Denoting the pressure in the vascular system by P v , the 3D flow model for asingle vessel V k reads as follows: u v = − K V k µ ( k )bl ∇ P v in V k ⊂ Ω v , (2)where u v is the velocity field within V k . This relationship is of Darcy or Hagen-Poiseuille type and can be derived from the Navier-Stokes equations assuminga stationary and laminar flow with a parabolic profile (see [3, Section 2.4] or[31]). It remains to specify the boundary conditions for both flow problems. Atfirst, we turn our attention to the boundary conditions on Γ coupling flow inthe microvascular network and the tissue. For this purpose, Starling’s filtrationlaw is considered to determine the flux across the vessel surface Γ k : u v · n = L p (( P v − P t ) − σ ( π v − π t )) =: J p ( P v , P t ) on Γ k , (3)where n is the outward unit vector normal to the vessel surface Γ k and P t represents the fluid pressure in the tissue. In (3), besides the pressure differencebetween fluid pressures, the oncotic pressures π v and π t are also taken intoaccount, where the oncotic pressure difference is weighted by a constant oncoticreflection coefficient σ . Finally, both pressure differences are multiplied by thehydraulic conductivity L p of the vessel wall. For simplicity, we assume that L p is equal and constant for all the vessels contained in the network Ω v .Next, we apply a topological model reduction transforming the 3D flow prob-lem (2) into a 1D flow problem. Since the radius of V k , k ∈ { , . . . , N } is smallcompared to the diameter of Ω, we assume that the function u v has a uni-form profile with respect to a cross section D k ( s ) i.e. P v is constant and equalto p v ( s ) on D k ( s ) for all s ∈ (0 , l k ). Integrating over an arbitrary cylinder Z k ⊂ V k having the radius R k and ranging from s to s , where s < s and9 , s ∈ (0 , l k ), we obtain by (2) and the divergence theorem:0 = (cid:90) Z k ∇ · (cid:32) K V k µ ( k )bl ∇ P v (cid:33) dV = (cid:90) ∂Z k K V k µ ( k )bl ∇ P v · n dS == − (cid:90) D k ( s ) K V k µ ( k )bl ∂P v ∂s dS + (cid:90) D k ( s ) K V k µ ( k )bl ∂P v ∂s dS + (cid:90) Γ Zk K V k µ ( k )bl ∇ P v · n dS. Γ Z k is the lateral surface of Z k . By means of the interface condition (3), thelast surface integral can be further simplified to: − (cid:90) Γ Zk K V k µ ( k )bl ∇ P v · n dS = (cid:90) Γ Zk L p (( P v − P t ) − σ ( π v − π t )) dS == L p (cid:90) Γ Zk P v dS − L p (cid:90) Γ Zk P t + σ ( π v − π t ) dS == 2 L p πR k (cid:90) s s p v ( s ) ds − L p (cid:90) Γ Zk P t + σ ( π v − π t ) dS. Using the fundamental theorem of integral calculus and K V k = R k /
8, we obtain: − (cid:90) D k ( s ) K V k µ ( k )bl ∂P v ∂s dS + (cid:90) D k ( s ) K V k µ ( k )bl ∂P v ∂s dS = R k π µ ( k )bl (cid:90) s s ∂ p v ∂s ds. Summing up the previous equations, a 1D integral equation for the unknown p v results:0 = (cid:90) s s R k π µ ( k )bl ∂ p v ∂s − L p πR k p v ds + L p (cid:90) Γ Zk P t + σ ( π v − π t ) dS = (cid:90) s s R k π µ ( k )bl ∂ p v ∂s − L p πR k p v + L p (cid:90) ∂ D k ( s ) P t + σ ( π v − π t ) dγ ds. Since s and s are chosen arbitrarily, we conclude by the fundamental theoremof variational calculus that the following 1D differential equation holds for eachvessel: 0 = R k π µ ( k )bl ∂ p v ∂s − πL p R k p v + L p (cid:90) ∂ D k ( s ) P t + σ ( π v − π t ) dγ. (4)Using (3), we get: − R k π µ ( k )bl ∂ p v ∂s = − π R k J p (cid:0) p v , P t (cid:1) , where P t ( s ) = 12 πR k (cid:90) ∂ D k ( s ) P t dγ. This means that the flow equation for V k is now defined on the 1D centerlineΛ k parameterised by s ∈ (0 , l k ). As a consequence, the 3D network Ω v is fromnow on considered as a 1D network Λ, which is the union of the 1D edges E k . Itremains to specify the boundary conditions for x ∈ ∂E k . Thereby, we considerthe following three cases (see Figure 2, Step 4): • Boundary node: x ∈ ∂E k ∩ ∂ ΛAt boundary nodes, we set a Dirichlet boundary value: p v ( x ) = p boundary .10 Inner network node: x ∈ ∂E k ∧ x / ∈ ∂ ΛAt a boundary node located in the inner part of the network Λ, we enforcethe continuity of pressure and the conservation of mass. Let N ( x ) be theindex set numbering the edges E l whose boundaries are intersecting at x i.e.: x = ∂E k ∩ ∂E l . Using this notation, the coupling conditions read asfollows: – Continuity of p v : p v | x ∈ ∂E k = p v | x ∈ ∂E l , ∀ l ∈ N ( x ) . – Mass conservation: (cid:88) l ∈ N ( x ) ∪{ k } (cid:32) − R l π µ ( l )bl ∂p v ∂s (cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x (cid:33) = 0 . A derivation of these coupling conditions (mass conservation, continuity ofpressure) can be found in [13]. In this publication, the coupling conditionsfor a 1D hyperbolic blood flow problem have been derived by means of anenergy inequality. Similar concepts can be applied to derive the couplingconditions for an elliptic flow problem like (7).After deriving a simplified flow model for the blood vessel network, we turn ourattention to flow in the porous matrix. Since we reduced our 3D network Ω v to a1D network given by Λ, we identify the porous tissue matrix Ω t with the wholedomain Ω. According to (A6), it is assumed that the tissue is an isotropic,homogeneous porous medium. Consequently, the corresponding permeabilitytensor K t is given by a constant scalar K t : K t ( x ) = K t · I for x ∈ Ω . Similarly, the viscosity of the fluid in the interstitial space is taken to be constantand is denoted by µ int . This can be justified by the fact that this fluid consistsmostly of water which is a Newtonian fluid. At the boundary ∂ Ω, we set forsimplicity homogeneous Neumann boundary conditions: n · (cid:18) K t µ t ∇ P t (cid:19) = 0 on ∂ Ω . Posing a no-flow boundary condition implies that the tissue block is supplied bythe circulatory system via the arterial blood vessels and drained via the venouspart of the network.To model the fluid exchange with the vascular system, we adapt a couplingconcept discussed in [5, 25, 26]. Thereby the influence of the vascular system isincorporated by means of a source term in the corresponding Darcy equation.The source term contains a Dirac measure concentrated on the lateral surfaceΓ of the network. This is motivated by the observation that the fluid exchangebetween the vascular system and tissue occurs across the vessel surface. An-other coupling concept consists in concentrating the Dirac source term on thecenterline Λ of the blood vessel network. However, when using this strategy,line singularities along the centerline [6, 7] arise, which have to be handled withcare. To determine the complete source term for the tissue matrix, we multi-ply the Dirac measure by Starling’s law (3). Combining Darcy’s equation and11 conservation equation, we have the following flow model for the 3D tissuematrix: −∇ · (cid:18) K t µ t ∇ P t (cid:19) = J p (Π ( p v ) , P t ) δ Γ in Ω , n · (cid:18) K t µ t ∇ P t (cid:19) = 0 on ∂ Ω . (5)Since p v is defined on a 1D manifold, we have to project it onto the vesselsurface such that the pressure difference between the vascular pressure and thetissue pressure can be computed. For an edge E k ⊂ Λ the projection operatorΠ occurring in (5) is defined by:Π ( p v ) | ∂ D k ( s ) ≡ p v ( s ) , ∀ s ∈ (0 , l k ) . Contrary to existing 3D-1D coupled flow models [4, 5, 6, 7, 25], we project bymeans of the operator Π the 1D pressure onto the vessel walls, to compute thepressure differences on Γ k . In the aforementioned references, the 3D pressuresare projected on the 1D vessels E k , by means of an average operator [7] (seeFigure 3). Summarising (4) and (5), we finally have the following 3D-1D coupledFigure 3: In this figure, we illustrate two possible coupling concepts for 3Dtissue and 1D vascular flow. On the left hand side, the tissue pressure p t isaveraged and projected on the main axis E k of the vessel V k . The right handside shows the projection of the 1D vascular pressure p v ( s ) onto the vessel wallof V k .PDE-system governing flow in the vascular system and the tissue: −∇ · (cid:18) K t µ t ∇ p t (cid:19) = J p (Π ( p v ) , p t ) δ Γ in Ω , (6) − R k π µ ( k )bl ∂ p v ∂s = − πR k J p ( p v , p t ) , ∀ E k ⊂ Λ , (7)where J p ( P v , P t ) = L p (( P v − P t ) − σ ( π v − π t )) . Henceforth, we adopt a moreuniform notation and denote the 3D tissue pressure as well as the 1D pressurein the vascular system by lower case letters, i.e. P t is replaced by p t in (6). p v in (7) still indicates the 1D blood pressure. A similar 3D-1D coupled flow modelcan be found in [23]. One major difference between the PDE-system in (6) and127) and the model in [23] is the type of coupling vascular and tissue flow. Ourmodel is based on Dirac source terms in 3D which are concentrated on the vesselsurfaces. In [23] this is not the case. Here, the 3D source term is concentrated onthe center lines of the microvascular network. However, modelling the exchangeprocesses between tissue and vascular system based on the vessel surfaces as itis done in (6) and (7) is more realistic, since the exchange processes occur atthe permeable vessel surfaces. To determine the oxygen distribution in both the blood vessel network and thesurrounding tissue matrix, we assign to each system a variable quantifying thepartial pressures. In the following, we denote them by P ( v ) O for the vascularsystem and by P ( t ) O for the tissue block. As the fluid pressures, P ( v ) O is definedon the 1D network Λ, and P ( t ) O is defined on Ω. Concerning the distribution ofthe partial pressures, we assume that it is governed by convection and diffusion,where the diffusion coefficients for the vascular system and the tissue are givenby two constants D v and D t .As in Section 2.2, we use a 3D-1D coupled model for these processes. Thederivation is similar to the one developed for the blood flow problem. For con-venience, we briefly describe here the most important steps: First, the networkis decomposed and a 1D PDE for each vessel is derived. Next, we assume thatthe vascular partial pressure is assumed to be constant across the section area.Based on this assumption a 3D convection-diffusion equation for each vesselis reduced to a 1D convection-diffusion equation using similar arguments as inSubsection 2.2. After that suitable coupling and boundary conditions are chosenfor the boundary nodes of the 1D network.The 3D model for the tissue is again defined on Ω. Similarly for the bloodflow process, the exchange between tissue and vascular system is accountedfor by means of a Dirac measure concentrated on the lateral surfaces Γ k ⊂ Γ.Modelling the flux of oxygen from the vessels into the tissue matrix, we considerthe vessel walls as semipermeable membranes allowing a selective filtration ofmolecules. According to [4, Section 5.2] the Kedem-Katchalsky equation is asuitable model to determine the flux J P O of oxygen across the vessel wall: J P O (cid:16) Π ( p v ) , p t , Π (cid:16) P ( v ) O (cid:17) , P ( t ) O (cid:17) = (1 − σ ) J p (Π ( p v ) , p t ) P ( t/v ) O (8)+ L P O (cid:16) Π (cid:16) P ( v ) O (cid:17) − P ( t ) O (cid:17) , ∀ Γ k ⊂ Γ . The symbol P ( t/v ) O denotes the average partial pressure in vessel walls. In thiswork, we define it as the arithmetic mean of the partial pressures on the twosides of the walls [4, Chapter 5.2] [37]: P ( t/v ) O = 12 (cid:16) P ( t ) O + Π (cid:16) P ( v ) O (cid:17)(cid:17) , The parameter L P O in (8) indicates the permeability of the vessel wall withrespect to oxygen. As mentioned in Section 2.1, we use the Michaelis-Menten13ormula [55], m ( ox ) (cid:16) P ( t ) O (cid:17) = m ( ox )0 P ( t ) O + P ( t ) O , P ( t ) O (9)to model the oxygen consumption of the tissue cells. m ( ox )0 represents the max-imal oxygen demand. For this work, we assume that it is constant. P ( t ) O , is thepartial pressure at half maximal consumption. In this work, we assume thatboth parameters are constant. Summarising the considerations of this subsec-tion and using (8) and (9), we obtain the following 3D-1D coupled model foroxygen transport: ∇ · (cid:16) u t P ( t ) O − D t ∇ P ( t ) O (cid:17) = − m ( ox ) (cid:16) P ( t ) O (cid:17) (10)+ J P O (cid:16) Π ( p v ) , p t , Π (cid:16) P ( v ) O (cid:17) , P ( t ) O (cid:17) δ Γ in Ω ,∂∂s (cid:32) u v P ( v ) O − D v ∂P ( v ) O ∂s (cid:33) = − πR k J P O (cid:16) p v , p t , P ( v ) O , P ( t ) O (cid:17) , ∀ E k ⊂ Λ , (11) n · (cid:16) u t P ( t ) O − D t ∇ P ( t ) O (cid:17) = 0 on ∂ Ω . The velocities u v and u t are computed by means of Darcy’s law: u v | Λ k = − R k µ ( k )bl ∂p v ∂s and u t = − K t µ t ∇ p t . (12)As for the fluid problem, we and use Dirichlet conditions for the boundary nodes x ∈ ∂E k ∩ ∂ Λ. Thereby, we have to check first, if E k is part of a vein or anartery (see also Subsection 3.1). Depending on the decision a typical value fora venous or an arterial partial pressure is chosen. In this work, we compute thefluid pressures and velocities and consider a vessel as vein, if its fluid pressure orvelocity is below the corresponding average. This is motivated by the fact thatveins belong to the low pressure and velocity region in the vascular system. Inthe other case, the vessel is considered as an artery. At the inner networks nodes( x ∈ ∂E k ∧ x / ∈ ∂ Λ), the continuity of partial pressures and a mass conservationequation is considered. In order to formulate the coupling equations, we usethe same notation as for the fluid problem and the coupling conditions for thepartial pressure of oxygen at inner network nodes read as follows: • Continuity of P ( v ) O : P ( v ) O (cid:12)(cid:12)(cid:12) x ∈ ∂E k = P ( v ) O (cid:12)(cid:12)(cid:12) x ∈ ∂E l , ∀ l ∈ N ( x ) , • Mass conservation: (cid:88) l ∈ N ( x ) ∪{ k } (cid:32) u v P ( v ) O − D v ∂P ( v ) O ∂s (cid:33)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) x ∈ ∂E l = 0 . .4 Generation of surrogate microvascular networks The process of generating a surrogate microvascular network out of the well-segmented arterioles and venoles is divided into three different phases:
Phase 1 (P1)
Generation of the larger blood vessels i.e. blood vesselswith a radius ranging from about 4 . µm to 16 . µm . Phase 2 (P2)
Generation of the fine scaled vessels i.e. blood vessels witha radius ranging from about 2 . µm to 4 . µm . Phase 3 (P3)
Removal of terminal vessels ending in the inner part of thetissue block.This is motivated as follows: Since even a microvascular network like the onedepicted in Figure 1, exhibits a multiscale structure, it is reasonable to firstgenerate the large scaled vessels in a first step. Otherwise, the growth of thelarger vessels might be blocked by some smaller vessels and the main vesselsof the network can not be developed. In a second step, we add the fine scaledvessels sprouting out of the large scaled vessels. The lower bound of 2 . µm Algorithm 1:
Collision detection algorithm
Preliminary work:
The computational domain Ω is decomposed intoeight disjunctive sub-control volumes in form of cuboids SV l , l ∈ { , . . . , } . Next the edges of the 1D network Λ are assigned tothe SV l in which they are contained. If a new edge E new is contained in SV l only the edges in SV l have to be checked for intersection and a loopover all edges in Λ is avoided; INPUT:
New edge E new to be added to Λ, 1D network Λ, sub-controlvolumes SV l , l ∈ { , . . . , } ;Determine the sub-control volumes SV l with E new ∩ SV l (cid:54) = ∅ ;Add E new to the detected sub-control volumes SV l ;Get radius R new of the edge E new ;isIntersecting = false ; for all edges E k ∩ SV l (cid:54) = ∅ do Get radius R k of the edge E k ;Compute the eucledian distance dist nk between E new and E k ; if E new and E k do not share a common network node ANDdist nk < R new + R k thenOUTPUT: E new is intersecting an existing edge;isIntersecting = true ;Remove E new from SV l ; BREAK ; endendif !isIntersecting thenOUTPUT: E new is intersecting no existing edge;Add E new to Λ ; end for the vessel radii in (P2) is motivated by considering the distribution of radii15n Figure 6. It can be observed that most of the radii are larger than 2 . µm .Therefore, we choose this value as a lower bound in (P2) . In order to preventan excessive runtime of our iteration procedure, we restrict the total number ofgrowth steps for (P1) and (P2) by N (1 , it, max = 35 and for (P3) by N (3) it, max = 15.Both the generation of the larger vessels and the generation of the smaller vesselsis carried out in an stepwise manner (see Subsection ?? ), where new vessels areonly added at the outlets of the terminal vessels. Thereby, both arterial andvenous vessels are updated.However, after finishing Phase (P2) , it can be observed that in the innerpart of the tissue there are several vessels, which are not fully integrated in thenetwork. In a healthy microvascular network this is usually not the case, sinceit exhibits fine scaled vessels or a capillary bed forming continuous connectionsbetween the arterial and venous part of the microvascular network. An impor-tant tool that is used several times in the algorithm below, is a method to detectcolliding or intersecting vessels (see Algorithm 1). Avoiding the intersection ofblood vessels is crucial for generating a well structured blood vessel network. Inthe remainder of this subsection, we provide the key ingredients for the threephases listed above: • Phase 1 (P1):
Given a 1D network Λ j from the j -th step in (P1), we compute the fluidpressures in Ω as described in Subsection 2.2. Please note that for j = 0,Λ j is given by the large scaled vessels Λ lv presented in Figure 1 (right).Based on the velocity fields, the distribution of the partial pressures inboth tissue and vascular system is determined using the PDEs providedin Subsection 2.3. Next, we add new blood vessels at the boundary nodesof the network Λ j . This results in a new 1D vascular network Λ j +1 forthe next step. Let us consider a boundary node x contained in an edge E ( j ) k ⊂ Λ j : x ∈ ∂ Λ j ∩ ∂E ( j ) k . For each boundary node x the normalisedgradient of P ( t ) O is determined: d p ( x ) = ∇ P ( t ) O ( x ) (cid:13)(cid:13)(cid:13) ∇ P ( t ) O ( x ) (cid:13)(cid:13)(cid:13) . (13)Choosing d p as the new direction of growth reflects the property of mi-crovascular networks to grow towards the tissue region with an increasedmetabolic demand. However, our numerical tests showed that followingstrictly d p can result into sharp bendings posing a high resistance on bloodflow. Therefore, we follow an idea presented in [54]. In this publication,the orientation d k of E ( j ) k is multiplied with a regularisation parameter λ g and combined with d p to avoid such sharp bendings: d g ( x ) = d p ( x ) + λ g · d k ( x ) , (14)where d g ( x ) denotes the new direction of growth at the node x . Normal-ising the vector d g ( x ), we obtain the new growth direction. It remainsto specify the radius and length of the new segment. If a single vesselis generated out of E ( j ) k , the vessel radius is taken over from the fathersegment E ( j ) k . To determine the length l of the new vessel, we assume a16og-normal distribution of the ratio r = l/R between the vessel length andthe vessel radius:log N ( r, µ r , σ r ) = 1 r (cid:112) πσ r exp (cid:32) − (ln r − µ r ) σ r (cid:33) , (15)where µ r and σ r denote the mean value and standard deviation of thelog-normal distribution, respectively. In accordance with the distributionin (15), we compute r . Using the radius R k of the father segment E ( j ) k ,the length l of the new vessel is given by: l = R k · r .Besides the formation of a single vessel at the inner boundary of E ( j ) k , weconsider also the possibility of generating a bifurcation. As in [54], thecumulative distribution function P bif of (15) is used to decide whether abifurcation or a single vessel is created. The function P bif is given by thefollowing expression: P bif ( r ) = Φ (cid:18) ln r − µ r σ r (cid:19) = 12 + 12 erf (cid:32) ln r − µ r (cid:112) σ r (cid:33) . (16)In this context, Φ represents the standard normal cumulative distributionfunction and erf the Gaussian error function. If P bif ( r ) exceeds for E ( j ) k a certain threshold P th = 0 .
6, a bifurcation is attached to the boundarynode. We observed that this threshold leads to a reasonable number ofbifurcations in the surrogate network. If P bif ( r ) ≤ P th a single vesselis generated as described before. For our simulations, we used µ r = 2 . σ r = 0 .
3. The choice of µ r and σ r results from several numericalsimulations. Using these values our vascular trees exhibit radii and lengthsin a reasonable range, as we will discuss in Section 3. At this point, wedeviate from [54]. The authors of this publication use the following values: µ r = 14 . σ r = 6 .
0. Now the issue arises of, how to choose the radii,orientations and lengths of the new branches b and b . The radii of thenew branches are computed based on a Murray type law, which relatesthe radius R k of the father vessel to the radius R b of branch b and theradius R b of branch b [33]: R γk = R γb + R γb , (17)where γ denotes the bifurcation exponent. However, it is not quite clear,how to choose this parameter. In [15, 33, 54, 63] one can find valuesfor γ ranging from 2 . .
5. Our numerical simulations show that thisparameter influences the shape of the microvascular network significantly.Therefore, we study in Subsection 3.2 how a variation of this parameteraffects the morphology of the network. Besides (17), a further equationis required to determine the radii of the branches. According to [54], weassume that R b follows a Gaussian normal distribution: R c = 2 − γ R k , R b i ∼ N ( R, µ = R c , σ = R c / , i ∈ { , } . (18)The heuristic choice of R b i has the following interpretation: Based on theradius R k of the father vessel, we compute the expected radius R c result-ing from Murray’s law for a fully symmetric bifurcation ( R b = R b ). R c
17s used as a mean value for a Gaussian normal distribution, with a smallstandard derivation. By this, we obtain bifurcations that are slightly de-viating from a symmetric bifurcation whose radii are chosen in accordancewith Murray’s law. Having R b and R b at hand, we compute the corre-sponding lengths l b and l b as in the case of a single vessel. The creationof a bifurcation is accomplished by specifying the orientations of the twobranches. At first, we define the plane, in which the bifurcation is con-tained. The normal vector n p of this plane is given by the cross product ofthe vessel orientation d k and the growth direction d g of the non-bifurcatingcase: n p = d k × d g (cid:107) d k × d g (cid:107) . The exact location of the plane is determined such that the vessel E ( j ) k iscontained in this plane. Further constraints for the bifurcation configura-tion are related to the bifurcation angles. In [32] and [51], it is shown howoptimality principles of minimum work and minimum energy dissipationcan be used to derive formulas relating the radii of the branches to thebranching angles φ (1) k and φ (2) k :cos (cid:16) φ (1) k (cid:17) = R k + R b − R b · R k R b and cos (cid:16) φ (2) k (cid:17) = R k + R b − R b · R k R b . (19) φ ( i ) k denotes the bifurcation angle between branch i, i ∈ { , } and thefather vessel (see Figure 4, right). Rotating the vector d k at x aroundthe axis defined by n p counterclockwise by φ (1) k and clockwise by φ (2) k , weobtain two new growth directions d b and d b . These vectors can be usedto define the main axes of the two cylinders representing the two branches.However, this approach does not allow us to incorporate the direction ofthe local oxygen gradient into the growth direction. For this reason theoptimality conditions (19) for the bifurcation angles are relaxed. This isdone by choosing the vector d b i , i ∈ { , } minimising the difference to d g : d ( k )min = argmin d bi (cid:107) d b i − d g (cid:107) . As orientations for the bifurcation, wekeep the direction d b j (cid:54) = d ( k )min , j ∈ { , } , while the other orientation isgiven by a vector d b j splitting the angle between d ( k )min and d g (see Figure,4, right) in two halves. This choice for the second growth direction can beconsidered as a compromise between the optimality principles providedby (19) and the tendency of the network to adapt its growth directionto the oxygen demand of the surrounding tissue. Before the new vesselsare attached to the inner boundaries of Λ j , Algorithm 1 is used to testwhether the newly created vessels are intersecting an existing vessel. Ifa vessel is attached to a boundary node x , we obtain a new boundarynode y for Λ j +1 . The new boundary values for p v and P O ( v )2 are takenover from the old boundary node x . In the final step of (P1) , we check,whether there is a terminal vessel left at which a large scaled vessel can beattached. A terminal vessel E ( j +1) k ⊂ Λ j +1 is considered as a large scaledvessel, if its radius R k is larger than 4 . µm . In addition, it is observedwhether the P O ( t )2 distribution is changing significantly. In order to testthis, we consider the region of interest Ω roi ⊂ Ω and compute an averaged18alue
P O ( t )2 ,roi : P O ( t )2 ,roi = 1 | Ω roi | (cid:90) Ω roi P O ( t )2 dV. (20)If the relative change of P O ( t )2 ,roi with respect to the previous iterationis below 1 . P O ( t )2 distribution as stationary. Phase(P1) is summarised in Algorithm and the result of (P1) is denoted by Λ P .Figure 4: The picture on the left hand side shows the plane defined by thenormal n p and the boundary node x . Furthermore the orientation vector d k of the vessel E ( j ) k is shown. d p represents a vector pointing into the directionof the local gradient of the partial pressure at x . It can be seen that d p formsa sharp bending at x with respect to E ( j ) k . A linear combination of d p with λ g d k results in a smoothing of this sharp corner. Thereby, λ g plays a role of aregularisation parameter. On the right hand side the two bifurcation angles φ (1) k and φ (2) k resulting from Murray’s law (19) are shown. The vector d b j representsthe direction of growth resulting from the regularised vector d g and the direction d b obtained from Murray’s law. It can be considered as a compromise betweenthe metabolic demand and the optimality principles given by (19). • Phase (P2):
Having the 1D network Λ P from the first growth phase at hand, we com-pute as in (P1) the fluid pressures as described in Subsection 2.2. Basedon the resulting velocity fields, the distribution of the partial pressures inboth tissue and vascular system is determined using the PDEs providedin Subsection 2.3. The third step of the second phase uses essentially thesame concepts presented in (P1). An exception is the choice of the radii.If the computed radius is below a threshold of 3 . µm , we deviate fromthe strategy discussed in (P1). This is motivated by the observation thata huge number of vessels becomes too tiny, if we follow strictly the stepsof (P1). Instead, it is assumed that a new radius R k < . µm followsanother normal distribution: R k ∼ N ( R, µ = 2 . µm, σ = 0 . µm ) . (21)By this choice, the distribution of the radii is mainly concentrated on therange of the fine scale vessels i.e. [1 . µm, . µm ]. In addition to that, wedemand that 2 . µm ≤ R k ≤ R p , where R p denotes the radius of the fathervessel. As it has already been mentioned, this is motivated by the fact that19 lgorithm 2: Algorithm for Phase 1
INPUT:
1D network Λ lv consisting of the larger arteries and veins,domain Ω ;Set j = 0, Λ j ← Λ lv , P O ( t )2 ,roi,old = 0, StopPhase1 = false ; while !StopPhase1 do Determine the pressures p v and p t by solving (6) in Ω and (7) in Λ j ;Compute the partial pressures P ( v ) O and P ( t ) O by solving (6) in Ω and(7) in Λ j ; for all boundary nodes x ∈ Ω ∩ ∂ Λ j do Let x ∈ ∂E ( j ) k , get radius R k of edge E ( j ) k ⊂ Λ j ; if R k > . µm then Compute the growth direction d g = d p + λ g d k according to(14);Compute the length and radius of the new vessel according to(15);Evaluate P bif according to (16); if P bif > P th then Create a second branch by means of Murray’s laws (17),(18) and (19);Intersection tests (Alg. 1) for the new vessels, adapt theboundary conditions; else
Intersection test (Alg. 1) for the new vessel, adapt theboundary conditions; endendend
Add the new vessels and update the network Λ j +1 ← Λ j , compute P O ( t )2 ,roi,j according to (20), determine the number of the largeterminal vessels ( N term,lv ) ; if (cid:12)(cid:12)(cid:12) P O ( t )2 ,roi,j − P O ( t )2 ,roi,old (cid:12)(cid:12)(cid:12) /P O ( t )2 ,roi,j < . · − OR j > N (1) it, max OR N term,lv = 0 then StopPhase1 = true ; end P O ( t )2 ,roi,old = P O ( t )2 ,roi,j , j = j + 1; endOUTPUT: Network containing the main vessels: Λ P ← Λ j ;20lmost all the vessels of the segmented network are larger than 2 . µm (seeFigure 6). The upper bound R k ≤ R p ensures that the radius R k of thenew vessel does not exceed the radius R p of the father vessel. Applyingthis strategy, the main part of the vessel radii are within a reasonablerange and the distribution of the radii has a smooth shape. Before thenew vessels are attached to the boundaries of Λ j , we test as in (P1.4) ifa new vessel intersects an existing vessel or the boundary ∂ Ω. Dependingon the result it is decided whether a certain vessel can be integrated intothe new network. By this, we obtain a new network ˜Λ j . If a vessel isattached to a boundary node x , we obtain a new boundary node y for˜Λ j . The new boundary values for the blood pressure p v and the partialpressure of oxygen P O ( v )2 are taken over from the old boundary node x . In new vessel connecting two terminal vessels Figure 5: This figure shows a terminal vessel ˜ E ( j ) k with a boundary node x .For convenience, we show a 2D sketch of a part of the network ˜Λ j . Connectingthe boundary node x with other network nodes that are contained in ˜Λ j , weconsider only those network nodes that are within a cone of opening angle π .From these network nodes, we choose those that have a distance less than d x from x (red nodes) and connect one of them with x . Thereby, we considerpreferably nodes of minimal distance to x and of maximal pressure differencebetween x and the other node. The node connected with x is denoted by y ,such that x and y form a new blood vessel that can be added to the networkΛ j +1 (red segment).a next step, we check if the new inner boundaries of ˜Λ j can be connectedwith another network node. For this purpose the network nodes in thevicinity of a boundary node x ∈ ∂ ˜Λ j are taken into account. Thereby,only network nodes that are contained in a cone with an opening angle of2 π/ x . By this convention, the formation ofsharp bendings can be circumvented. From the network nodes within thecone, we choose only those network nodes as candidates that have at mosta distance d x to x . In order to create a smooth distribution of the vessel21engths, it is assumed that d x can be chosen according to the followingdistribution: d x ∼ N ( d x , µ = 60 . µm, σ = 10 . µm ) . This is motivated by the fact that blood vessels in a microvascular networkgrow at least about 50 . µm per day [44]. Among the network nodesfullfilling these conditions, we search for those with minimal distance andthe highest pressure difference. The pressure difference is computed bymeans of the pressure value at x and the pressure value at the networknode that could be connected with x . As a consequence, we link preferablyvessels from the venous and the arterial side of the microvascular network,in order to enable a circulation of blood within the microvascular system.The radius of the new segment or blood vessel is given as the arithmeticmean of the terminal vessel containing x and the segment containing theother network node (see Figure 5). Next, we check by means of Algorithm1, whether they intersect other vessels and adapt the boundary conditions.In the final step of (P2), Ω roi is decomposed in N cv control volumes CV i having the shape of a cuboid: Ω roi = (cid:83) N cv i =1 CV i . For our computations,we used in each space direction 4 control volumes i.e. N cv = 64. Withrespect to each CV i and Ω roi , we compute averaged P O ( t )2 values: P O ( t )2 ,i = 1 | CV i | (cid:90) CV i P O ( t )2 dV and P O ( t )2 ,roi = 1 | Ω roi | (cid:90) Ω roi P O ( t )2 dV. (22)Having these values at hand, the following decisions are made: If P O ( t )2 ,i > . CV i . By this an adaptivegrowth of the network is achieved. The whole growth process in thisphase is stopped, if P O ( t )2 ,roi > . . P O ( t )2 isaround 32 . − . . P O ( t )2 = 34 . • Phase (P3):
In the third phase of the generation process, the network structure Λ P obtained from (P2) is further optimised by removing dead ends i.e. ter-minal vessels contained in Ω roi . In addition to that, we try to connectterminal vessels with other vessels in the network. This is again executedin an iterative way by means of the following four steps: In the first step,we remove all the terminal vessels in Ω roi .After removing the terminal vessels, we try to connect the new terminalvessels with other vessels in the network as it is done in (P2). As in (P2),it is tested whether the new vessels are intersecting other vessels (see Al-gorithm 1). If less than 10 terminal vessels can be found in Ω roi , we stopthe growth process of the third phase. Finally we extract the networkcontained in Ω roi . A pseudo code for (P3) can be found in Algorithm 4.22 lgorithm 3: Algorithm for Phase 2
INPUT:
1D network Λ P obtained from (P1), domain Ω ;Set j = 0, Λ j ← Λ P , P O ( t )2 ,roi,old = 0, StopPhase2 = false ; while !StopPhase2 do Determine the pressures p v and p t by solving (6) in Ω and (7) in Λ j ;Compute the partial pressures P ( v ) O and P ( t ) O by solving (6) in Ω and(7) on Λ j ;Compute the local partial pressures P O ( t )2 ,k , k ∈ { , . . . , N cv } according to (22); for all boundary nodes x ∈ Ω ∩ ∂ Λ j do Determine CV k with x ∈ CV k ; if P O ( t )2 ,k > . then No new vessel is attached to x else Construct new vessels as in Algorithm 2, if the radii of thenew vessels are below 3 . µm choose the new radii accordingto (21) ; endend Add the new vessels, try to link terminal vessels within the resultingnetwork ˜Λ j as shown in Figure 5 and update the network Λ j +1 ← ˜Λ j ;Compute P O ( t )2 ,roi,j according to (20) ; if (cid:12)(cid:12)(cid:12) P O ( t )2 ,roi,j − P O ( t )2 ,roi,old (cid:12)(cid:12)(cid:12) < . · − OR j > N (2) it, max then StopPhase2 = true ; end P O ( t )2 ,roi,old = P O ( t )2 ,roi,j , j = j + 1 ; endOUTPUT: Network containing the main vessels and fine scaled vessels:Λ P ← Λ j ; 23 lgorithm 4: Algorithm for Phase 3
INPUT:
1D network Λ P obtained from (P2), domain Ω;Set j = 0, Λ j ← Λ P , P O ( t )2 ,roi,old = 0, StopPhase3 = false ;Determine the number of terminal vessels of Λ j in Ω roi : N term ; while !StopPhase3 dofor all boundary nodes x ∈ Ω ∩ ∂ Λ j doif x ∈ Ω roi then Remove vessel E ( j ) k containing x from Λ j ; endend Obtain an intermediate network ˜Λ j ;Try to link terminal vessels with other vessels in the network ˜Λ j asshown in Figure 5;Update the network Λ j +1 ← ˜Λ j ;Determine the number of terminal vessels of Λ j +1 in Ω roi : N term ; if j > N (3) it, max OR N term < then StopPhase3 = true ; end j = j + 1; end Extract network Λ roi = Λ j ∩ Ω roi ; OUTPUT:
Final network: Λ P ← Λ roi ; Concluding the main section of this work, we briefly describe the numericalmethods that have been used to determine a solution of the mathematical modelsdeveloped in the previous subsections. To solve the 3D PDEs (6) and (10),we use a standard cell-centered finite volume method [18]. The fluxes acrossthe boundaries of a finite volume cell are approximated by the two-point fluxmethod. Using a simple two-point approximation of the fluxes yields a consistentdiscretisation, since a simple uniform hexahedral mesh decomposing Ω is used.Moreover, we have assumed that the tissue is an isotropic porous medium andthat the diffusion coefficient for oxygen in tissue is constant. To incorporate theDirac source term in the 3D equations (6) and (10) using finite volumes, we firstdetermine which finite volume cells are intersected by the cylindrical vessels. Ina next step, the areas of vessel surfaces contained in the affected finite volumecells are estimated. Thereby, it is crucial that the area of a cylindrical vesselis consistently distributed among the respective finite volume cells. This meansthat the sum of the single surface areas assigned to the 3D finite volume cells isequal to the surface area of the respective cylinder. By this feature, we obtain amass conservative discretisation of the flow and transport models. An advantageof this numerical solution method is that the 1D graph structure governing thelocation of the vascular network does not have to be aligned with the edges ofthe 3D finite volume mesh as e.g. in [23] the authors require that the edges ofthe 3D mesh are adapted to the edges of the vascular graph.For the numerical solutions of the 1D PDEs (7) and (11), we employ the24ascular graph model (VGM) [12, 38, 48, 59], which corresponds in principle toa vertex centered Finite Volume method with a two-point flux approximation.The discretised versions of the 3D and 1D PDEs are linear in the fluid pressuresand partial pressures of oxygen. However, there is one exception caused by thesource term of (10). Since this source term involves the Michaelis-Menten law(9), a nonlinearity is introduced and a nonlinear solver is required. For oursimulations, we choose a simple fixed point iteration with a damping factor toenforce the convergence of the scheme. As an initial guess for the fixed pointiteration, we use the distribution of partial pressure of oxygen from the previousgrowth step.
After describing a numerical model for generating surrogate networks, we testthe performance of this numerical model. This is done by comparing someimportant features of the segmented network shown in Figure 1 (left) with thecorresponding features of the created network. In particular, the total length L = (cid:80) k ∈ I (cid:107) E k (cid:107) , the total vessel surface area A = (cid:80) k ∈ I (cid:107) E k (cid:107) R k π , thetotal volume V = (cid:80) k ∈ I (cid:107) E k (cid:107) R k π and the number of segments N seg = | I | of a network Λ = (cid:83) k ∈ I E k are considered. I is an index set numbering thesegments of Λ. The reason L is reported, is due to the fact that the surrogatenetwork should exhibit a similar flow path as the segmented network. If thesurface area A is in both networks in the same range, it can be expected thatunder similar boundary conditions the fluid exchange between the tissue and thevascular system is in both cases in the same range. The volume V is a quantityof interest, since networks having a similar volume contain a similar amount ofblood volume. In addition to that P O ( t )2 ,roi and an averaged fluid pressure intissue p t,roi : p t,roi = 1 | Ω roi | (cid:90) Ω roi p t dV are presented. Furthermore, the flux F tv [ µ g / s ] from the vascular system into thetissue and vice versa is calculated. Therefore, we multiply the volumetric fluxby the density of water ( ρ w = 1000 . kg / m ). Finally, we take the total numberof growth steps N it for all the three growth phases into account. All thesedata are reported varying two model parameters that affect the shape of thegenerated networks in a significant way. These parameters are the consumptionrate of oxygen m ( ox )0 and the Murray parameter γ . The rest of this section isdivided into two parts. In the first part, we list all the model parameters thatare used for the simulations, while in the second part, the simulation results arepresented and discussed. Figure 6 shows the distribution of radii and lengths of the segments that arepart of the given network Λ s shown in Figure 1 (left). By means of these distri-butions, we compute the corresponding length L s , surface area A s , volume V s and the total edge number N s . The different values can be found in Table 1.The next table contains all the remaining parameters occurring in the 3D-1Dcoupled PDE based models (7), (6) and (11), (10). It can be seen at a first25igure 6: This figure shows two histogramms with respect to the segmentednetwork shown in Figure 1. On the left hand side, the distribution of theradii is shown, while on the right hand side the distribution of the lengthsis presented. The terms avg and std stand for the mean value and standarddeviation, respectively.Table 1: Characteristic features of the segmented network (see Figure 1, left).Name Sign Value Unittotal length L s .
595 mtotal surface area A s . · − m total volume V s . · − m total number of segments N s,seg − glance that most of the parameters are fixed apart from m ( ox )0 and γ . We havechosen γ ∈ { . , . } , since these are usual values for the Murray parameter γ one can find in literature [15, 33, 54]. In [55, Table 1] the authors provide forthe oxygen consumption rate c ( ox )0 a range of 0 . O (cid:0)
100 cm (cid:1) − min − to 2 . O (cid:0)
100 cm (cid:1) − min − . By means of Henry’s law m ( ox )0 = c ( ox )0 /α H and the Henry constant α H = 2 . · − cm O cm − mmHg − [24, Table4], we obtain: m ( ox )0 ∈ [2 . mmHg / s , . mmHg / s ]. Due to the fact that a por-tion of rat brain is considered, it can be assumed that low oxygen consumptionrates should be used for the simulations. Therefore, we have chosen for m ( ox )0 the values 3 . mmHg / s and 4 . mmHg / s . The information on the networkΛ lv containing the larger and well segmented vessels can be found in the file extracted network.dgf , which is uploaded together with the supplementary ma-terial. In this file the nodes of the graph-like structure are listed. Next to thecoordinates of the boundary nodes a blood pressure in [Pa] is listed. These val-ues are used as boundary conditions for the 1D PDEs (7). In addition to that thesegments of Λ lv are defined i.e. there is a list consisting of three columns, wherethe first column contains the index of the first node and the second column theindex of the second node. Finally, the third column presents the radius of thecorresponding segment in [m]. It remains to specify the boundary conditions forthe partial pressure of oxygen P O ( v )2 . According to [55], the P O ( v )2 in arterial26able 2: Model parameters of the 3D-1D coupled models for blood flow andoxygen transport.Name Sign Value Unit SourcePermeability of tissue K t . · − m [4, Sec. 4.5]Permeability (plasma) L p . · −
12 m / Pa · s [4, Sec. 4.5]Viscosity blood plasma µ p . · − Pa · s [45]Reflection parameter σ . − − Onc. pressure blood plasma π v . π t . µ t . · − Pa · s [4, Tab. 5.1]Diff. of oxygen (vascular) D v . · − / s [4, Tab. 5.2]Diff. of oxygen(tissue) D t . · − / s [4, Tab. 5.2]Metabolic rate m ( ox )0 .
00, 4 . mmHg / s [55, Tab. 1] P O (half consump.) P ( t ) O , .
00 mmHg [55, Tab. 1]Permeability (oxygen) L P O . · − / s [4, Table 5.2]Murray parameter γ . , . − [54, Tab. 1]Regularisation parameter λ g . − [54, Tab. 1]Lower bound Ω roi , x L r ,l .
038 mm − Upper bound Ω roi , x L r ,u .
13 mm − Lower bound Ω roi , x L r ,l . · − mm − Upper bound Ω roi , x L r ,u .
05 mm − Lower bound Ω roi , x L r ,l . · − mm − Upper bound Ω roi , x L r ,u .
50 mm − vessels is around 75 . . lv and the average velocity for each edge. If an edge exhibits a velocitybelow the average velocity, we regard it as a vein. Otherwise, it is consideredas an artery. At first, we investigate, how many simulations are required to obtain, for eachparameter choice, a representative data set. This is necessary, since the growthprocess depends on several stochastic processes e.g. the formation of bifurca-tions. For this purpose, we perform N m (cid:16) γ, m ( ox )0 (cid:17) simulations and report allthe quantities of interest that are listed at the beginning of this section. Then,for each quantity q ∈ (cid:110) L, A, V, N seg , P O ( t )2 ,roi , p t,roi , F tv , N it (cid:111) , the followingmean values q m i are computed: q m i = 1 i i (cid:88) n =1 q n (cid:16) γ, m ( ox )0 (cid:17) , ≤ i ≤ N m (cid:16) γ, m ( ox )0 (cid:17) , q n (cid:16) γ, m ( ox )0 (cid:17) denotes the quantity q for a certain data set and the n -th simulation. As it can be observed numerically (see Figures 7 and 8), asmall number of samples is sufficient to find a quite good approximation ofthe quantities of interest. More precisely with only N m (cid:16) γ, m ( ox )0 (cid:17) = 20, weget for almost all cases satisfying results. Accordingly, we can assume that inthis case that we have a representative data set. In Table 3, the mean valuesand standard deviations for the different parameter pairs are listed. Comparingthese data with the ones in Table 1, we see that the quantities which determinethe morphology of the network are closest to the segmented ones for γ = 3 and m ( ox )0 = 3 . mmHg / s . All the other data sets deviate to a larger extend from thesegmented data.Next we compare for this parameter set two distributions of radii and lengthswith their segmented counterparts (see Figure 6 and Figure 9). It can be seenthat the distributions of the radii have a similar shape as the one of the seg-mented data and that the mean value as well as the standard deviation arerelatively close to each other. The main difference in this context is that morevessels are distributed in a small vicinity of the mean value. Contrary to that thedistributions of the edge lengths exhibit different behaviours. While in the sur-rogate case, almost all segments have lengths between 1 . µ m and 100 . µ m, thelength distribution in the segmented case is decaying much slower, which meansthat a significant amount of the edges has a length of more than 100 . µ m. Thisdifference might be explained by the fact that an edge length in the segmentednetwork does not necessarily correspond to a biological blood vessel, whereas inthe other case we try to adapt the edge length to a typical blood vessel length.However, the total length of the networks is quite similar due to the fact thatthe surrogate network contains more vessels than the segmented one (see Table1 and 3).Finally, we illustrate in Figures 10 and 11, how the different parameters affectthe shape of the networks and how the networks contained in Ω typically looklike after the different phases. The networks become denser, if the rate of oxygenconsumption is higher. This is due to the fact that for a higher consumption ratethe propagation front extends across a much shorter distance than for a lowerconsumption rate. In order to compensate the reduced distance and maintainan average P O ( t )2 of about 33 mmHg −
35 mmHg, the network has to becomedenser. As a consequence, the total surface area A of the vessels and the totalvolume V increases, which implies that the flux into tissue as well as the bloodvolume in the vascular system is increased. However, constructing the densernetworks requires more growth steps N it compared to the low consumption case(see Table 3). The blood pressure in the network ranges between 61 . . p t,roi , we have mean values between 33 . . roi . n n -5 n -11 n Figure 7: Asymptotic behaviour of the total length L , total surface area A ,total volume V and the total number of segments N seg .Table 3: Mean values and standard deviations of the quantities of interest. γ [ − ] 3 . . . . m ( ox )0 [ mmHg / s ] 3 . . . . L [m] 0 . ± .
145 0 . ± .
121 0 . ± .
159 0 . ± . A (cid:2) µ m (cid:3) · − . ± .
49 9 . ± .
10 1 . ± .
27 1 . ± . V (cid:2) m (cid:3) · − . ± .
33 1 . ± .
29 2 . ± .
38 2 . ± . N seg [ − ] 13845 ±
850 11689 ±
685 18459 ±
845 15327 ± P O ( t )2 ,roi [mmHg] 34 . ± . . ± . . ± . . ± . p t,roi [mmHg] 34 . ± . . ± . . ± . . ± . F tv [ µ g / s ] · − . ± .
28 7 . ± .
08 9 . ± .
00 8 . ± . N it [ − ] 34 . ± . . ± . . ± . . ± . In this paper, we have presented a model for simulating flow and oxygen trans-port processes in vascular systems embedded in a homogeneous tissue block.29 n n n -3 n Figure 8: Asymptotic behaviour of the pressures
P O ( t )2 ,roi , p t,roi , the flux intotissue F tv and the total number of growth steps N it .Using the distribution of the partial pressure of oxygen in tissue, a stochasticmodel for generating surrogate microvascular networks out of given arteriolesand venoles has been constructed. Thereby, the gradient in the oxygen distri-bution governs the growth direction while by means of log-normal distributionsthe length of the vessels and the formation of bifurcations is determined. Thechoice of the vessel radii at a bifurcation as well as the branching angles arebased on Murray’s law. Combining characteristic distributions for lengths andradii with Murray’s law and the oxygen gradient as the main growth direction,the resulting microvascular network incorporates several well-known optimalityprinciples. In addition, the algorithm avoids an intersection of newly createdvessels. Finally, it is guaranteed that in the interior part of the tissue domainno terminal vessels occur, i.e., the flow paths between the arterioles and venolesare not interrupted.Our simulation results reveal that about 20 simulations are needed to ob-tain a representative data set. Furthermore, we have observed that a Murrayparameter of γ = 3 and a low oxygen consumption rate yields a satisfactorymatch with respect to some characteristic quantities like the total surface areaor the volume of the vascular network.Future work in this field could be concerned with applying the presentednumerical model to further similar data sets i.e. microvascular networks thatare contained in a volume of a few cubic millimeters. For example it would beof great interest, how the values for the consumption rate of oxygen and theMurray parameter γ have to be chosen to obtain satisfactory results. Moreover,one could try to find a relationship for the radii and the branching angles at a30igure 9: Two distributions of radii and segment lengths for γ = 3 and m ( ox )0 =3 . mmHg / s .bifurcation that is more precise than Murray’s law, which is only an approxima-tion and does not hold in general. In addition to that, it could be investigatedwhether the remaining parameters like the mean values and standard deviationsof the statistical distributions as well as the probability threshold for bifurca-tions P th have to be recalibrated. Finally, the performance of the numericalmodel could be further tested by considering CT scans showing the saturationdistribution of a contrast fluid [50]. Besides the transport of oxygen furtherconvection-diffusion equations modelling the injection of a contrast fluid couldbe added and the resulting simulation data could be compared with the imagedata. Acknowledgements
This work was partially supported by the DFG grant (WO/671 11-1). Wegratefully acknowledge the authors of [48] and, in particular, Prof. Dr. BrunoWeber and Prof. Dr. Patrick Jenny, for providing the data set of the vascularnetwork considered in this paper.
References [1] J. Ahn, H. Cho, J. Kim, S. Kim, J. Ham, I. Park, S. Suh, S. Hong, J. Song,Y. Hong, et al. Meningeal lymphatic vessels at the skull base drain cere-31igure 10: This figure shows examples of four networks contained in Ω roi andobtained after the third phase. All the depicted networks are simulated by meansof different parameters and represent the typical results for the respective dataset. The colouring of the networks is given by the blood flow pressure rangingfrom 26 mmHg to 61 mmHg.brospinal fluid.
Nature , page 1, 2019.[2] P. Bastian, M. Blatt, C. Engwer, R. Kl¨ofkorn S. Kuttanikkad, andT. Neubauer M. Ohlberger O. Sander. Dune, distributed and unified nu-merics environment. , 2011.[3] S. ˇCani´c and E. Kim. Mathematical analysis of the quasilinear effects ina hyperbolic model blood flow through compliant axi-symmetric vessels.
Mathematical Methods in the Applied Sciences , 26(14):1161–1186, 2003.[4] L. Cattaneo.
FEM for PDEs with unfitted interfaces: application to flowthrough heterogeneous media and microcirculation . PhD thesis, Italy, 2014.32igure 11: In the top left picture of this figure the starting network for vesselgrowth is shown, while the top right picture shows the network after the firstgrowth phase. At the bottom the networks after the second phase (left) andthird phase (right) is presented. Each picture shows a tissue slice. The colourbar is adapted to the
P O within the tissue slice.[5] D. Cerroni, F. Laurino, and P. Zunino. Mathematical analysis, finiteelement approximation and numerical solvers for the interaction of 3dreservoirs with 1d wells. GEM-International Journal on Geomathematics ,10(1):4, 2019.[6] C. D’Angelo. Multiscale modelling of metabolism and transport phenomenain living tissues. Technical report, EPFL, 2007.[7] C. D’Angelo and A. Quarteroni. On the coupling of 1d and 3d diffusion-reaction equations: application to tissue perfusion problems.
MathematicalModels and Methods in Applied Sciences , 18(08):1481–1504, 2008.[8] J. Dankelman, A. Cornelissen, J. Lagro, E. VanBavel, and J. Spaan. Re-lation between branching patterns and perfusion in stochastic generatedcoronary arterial trees.
Medical & biological engineering & computing ,45(1):25–34, 2007. 339] M. Dewhirst and T. Secomb. Transport of drugs from blood vessels totumour tissue.
Nature Reviews Cancer , 17(12):738, 2017.[10] D. Drzisga, T. K¨oppl, U. Pohl, R. Helmig, and B. Wohlmuth. Numeri-cal modeling of compensation mechanisms for peripheral arterial stenoses.
Computers in biology and medicine , 70:190–201, 2016.[11] W. El-Bouri and S. Payne. Investigating the effects of a penetrating vesselocclusion with a multi-scale microvasculature model of the human cerebralcortex.
NeuroImage , 172:94–106, 2018.[12] K. Erbertseder, J. Reichold, B. Flemisch, P. Jenny, and R. Helmig. A cou-pled discrete/continuum model for describing cancer-therapeutic transportin the lung.
PLoS One , 7(3):e31966, 2012.[13] L. Formaggia, D. Lamponi, and A. Quarteroni. One-dimensional modelsfor blood flow in arteries.
Journal of engineering mathematics , 47(3-4):251–276, 2003.[14] L. Formaggia, A. Quarteroni, and A. Veneziani.
Cardiovascular Mathemat-ics: Modeling and simulation of the circulatory system , volume 1. SpringerScience & Business Media, 2010.[15] Y. Fung. Biomechanics: Circulation 2nd, 1997.[16] Y. Fung and B. Zweifach. Microcirculation: mechanics of blood flow incapillaries.
Annual Review of Fluid Mechanics , 3(1):189–210, 1971.[17] M. Georg, T. Preusser, and H. Hahn. Global constructive optimization ofvascular systems. 2010.[18] R. Helmig.
Multiphase flow and transport processes in the subsurface: acontribution to the modeling of hydrosystems.
Springer-Verlag, 1997.[19] W. Henry. Iii. experiments on the quantity of gases absorbed by water, atdifferent temperatures, and under different pressures.
Philosophical Trans-actions of the Royal Society of London , (93):29–274, 1803.[20] E. Hyde, C. Michler, J. Lee, A. Cookson, R. Chabiniok, D. Nordsletten, andN. Smith. Parameterisation of multi-scale continuum perfusion models fromdiscrete vascular networks.
Medical & biological engineering & computing ,51(5):557–570, 2013.[21] A. Khaled and K. Vafai. The role of porous media in modeling flow andheat transfer in biological tissues.
International Journal of Heat and MassTransfer , 46(26):4989–5003, 2003.[22] T. Koeppl, G. Santin, B. Haasdonk, and R. Helmig. Numerical modellingof a peripheral arterial stenosis using dimensionally reduced models andkernel methods.
International journal for numerical methods in biomedicalengineering , 34(8):e3095, 2018.[23] M. Kojic, M. Milosevic, V. Simic, E. Koay, J. Fleming, S. Nizzero, N. Kojic,A. Ziemys, and M. Ferrari. A composite smeared finite element for masstransport in capillary systems and biological tissue.
Computer methods inapplied mechanics and engineering , 324:413–437, 2017.3424] T. K¨oppl, M. Schneider, U. Pohl, and B. Wohlmuth. The influence of anunilateral carotid artery stenosis on brain oxygenation.
Medical engineering& physics , 36(7):905–914, 2014.[25] T. K¨oppl, E. Vidotto, B. Wohlmuth, and P. Zunino. Mathematical model-ing, analysis and numerical approximation of second-order elliptic problemswith inclusions.
Mathematical Models and Methods in Applied Sciences ,28(05):953–978, 2018.[26] J. Kremheller, A. Vuong, B. Schrefler, and W. A. Wall. An approach for vas-cular tumor growth based on a hybrid embedded/homogenized treatmentof the vasculature within a multiphase porous medium model.
Internationaljournal for numerical methods in biomedical engineering , 35(11):e3253,2019.[27] D. Lesage, E. Angelini, I. Bloch, and G. Funka-Lea. A review of 3d vessellumen segmentation techniques: Models, features and extraction schemes.
Medical image analysis , 13(6):819–845, 2009.[28] E. Lindstrøm, J. Schreiner, G. Ringstad, V. Haughton, P. Eide, andK. Mardal. Comparison of phase-contrast mr and flow simulations for thestudy of csf dynamics in the cervical spine.
The neuroradiology journal ,31(3):292–298, 2018.[29] B. Lloyd, S. Hirsch, and G. Sz´ekely. Optimization of case-specific vasculartree models based on vessel size imaging. In
International Symposium onBiomedical Simulation , pages 38–48. Springer, 2010.[30] K. Margaris and R. Black. Modelling the lymphatic system: challenges andopportunities.
Journal of the Royal Society Interface , 9(69):601–612, 2012.[31] S. Mayer. On the pressure and flow-rate distributions in tree-like andarterial-venous networks.
Bulletin of mathematical biology , 58(4):753–785,1996.[32] C. Murray. The physiological principle of minimum work applied to theangle of branching of arteries.
The Journal of general physiology , 9(6):835,1926.[33] C. Murray. The physiological principle of minimum work: I. the vascularsystem and the cost of blood volume.
Proceedings of the National Academyof Sciences of the United States of America , 12(3):207, 1926.[34] M. Nabil and P. Zunino. A computational study of cancer hyperthermiabased on vascular magnetic nanoconstructs.
Royal Society open science ,3(9):160287, 2016.[35] F. Nekka, S. Kyriacos, C. Kerrigan, and L. Cartilier. A model of growingvascular structures.
Bulletin of mathematical biology , 58(3):409–424, 1996.[36] D. Obrist, B. Weber, A. Buck, and P. Jenny. Red blood cell distribution insimplified capillary networks.
Philosophical Transactions of the Royal Soci-ety A: Mathematical, Physical and Engineering Sciences , 368(1921):2897–2918, 2010. 3537] K. Perktold, M. Prosi, and P. Zunino. Mathematical models of mass trans-fer in the vascular walls. In
Cardiovascular mathematics , pages 243–278.Springer, 2009.[38] M. Peyrounette, Y. Davit, M. Quintard, and S. Lorthois. Multiscale mod-elling of blood flow in cerebral microcirculation: Details at capillary scalecontrol accuracy at the level of the cortex.
PloS one , 13(1):e0189474, 2018.[39] C. Poelma. Exploring the potential of blood flow network data.
Meccanica ,52(3):489–502, 2017.[40] L. Possenti, S. di Gregorio, F. Gerosa, G. Raimondi, G. Casagrande,M. Costantino, and P. Zunino. A computational model for microcircu-lation including fahraeus-lindqvist effect, plasma skimming and fluid ex-change with the tissue interstitium.
International Journal for NumericalMethods in Biomedical Engineering , 35(3):e3165, 2019.[41] A. Pries, B. Reglin, and T. Secomb. Structural adaptation of vascularnetworks: role of the pressure response.
Hypertension , 38(6):1476–1479,2001.[42] A. Pries, B. Reglin, and T. Secomb. Remodeling of blood vessels: responsesof diameter and wall thickness to hemodynamic and metabolic stimuli.
Hypertension , 46(4):725–731, 2005.[43] A. Pries and T. Secomb. Microvascular blood viscosity in vivo and theendothelial surface layer.
American Journal of Physiology-Heart and Cir-culatory Physiology , 289(6):H2657–H2664, 2005.[44] A. Pries and T. Secomb. Making microvascular networks work: angiogen-esis, remodeling, and pruning.
Physiology , 29(6):446–455, 2014.[45] A. Pries, T. Secomb, and P. Gaehtgens. Biophysical aspects of blood flowin the microvasculature.
Cardiovascular research , 32(4):654–667, 1996.[46] NP Reddy and K Patel. A mathematical model of flow through the terminallymphatics.
Medical engineering & physics , 17(2):134–140, 1995.[47] J. Reichold.
Cerebral blood flow modeling in realistic cortical microvascularnetworks . PhD thesis, ETH Zurich, 2011.[48] J. Reichold, M. Stampanoni, A. L. Keller, A. Buck, P. Jenny, and B. Weber.Vascular graph model to simulate the cerebral blood flow in realistic vas-cular networks.
Journal of Cerebral Blood Flow & Metabolism , 29(8):1429–1443, 2009.[49] M. Rempfler, M. Schneider, G. Ielacqua, X. Xiao, S. Stock, J. Klohs,G. Sz´ekely, B. Andres, and B. Menze. Extracting vascular networks underphysiological constraints via integer programming. In
International Con-ference on Medical Image Computing and Computer-Assisted Intervention ,pages 505–512. Springer, 2014.[50] E. Rohan, V. Lukeˇs, and A. Jon´aˇsov´a. Modeling of the contrast-enhancedperfusion test in liver based on the multi-compartment flow in porous me-dia.
Journal of mathematical biology , 77(2):421–454, 2018.3651] R. Rosen.
Optimality principles in biology . Springer, 2013.[52] F. Schmid, J. Reichold, B. Weber, and P. Jenny. The impact of capillarydilation on the distribution of red blood cells in artificial networks.
Ameri-can Journal of Physiology-Heart and Circulatory Physiology , 308(7):H733–H742, 2015.[53] M. Schneider, S. Hirsch, B. Weber, G. Sz´ekely, and B. Menze. Tgif: Topo-logical gap in-fill for vascular networks. In
International Conference onMedical Image Computing and Computer-Assisted Intervention , pages 89–96. Springer, 2014.[54] M. Schneider, J. Reichold, B. Weber, G. Sz´ekely, and S. Hirsch. Tis-sue metabolism driven arterial tree generation.
Medical image analysis ,16(7):1397–1414, 2012.[55] T. Secomb, J. Alberding, R. Hsu, M. Dewhirst, and A. Pries. Angiogenesis:an adaptive dynamic biological patterning problem.
PLoS computationalbiology , 9(3):e1002983, 2013.[56] R. Shipley and S. Chapman. Multiscale modelling of fluid and drug trans-port in vascular tumours.
Bulletin of mathematical biology , 72(6):1464–1491, 2010.[57] I. Shvab and V. Nimaev. Mathematical modeling of microcirculatory pro-cesses. In , pages 531–533. IEEE, 2017.[58] K. Støverud, M. Darcis, R. Helmig, and M. Hassanizadeh. Modeling con-centration distribution and deformation during convection-enhanced drugdelivery into brain tissue.
Transport in porous media , 92(1):119–143, 2012.[59] E. Vidotto, T. Koch, T. K¨oppl, R. Helmig, and B. Wohlmuth. Hybridmodels for simulating blood flow in microvascular networks.
MultiscaleModeling & Simulation , 17(3):1076–1102, 2019.[60] J. Wagner. Properties of the michaelis-menten equation and its integratedform which are useful in pharmacokinetics.
Journal of Pharmacokineticsand Biopharmaceutics , 1(2):103–121, 1973.[61] K. Weishaupt, V. Joekar-Niasar, and R. Helmig. An efficient coupling offree flow and porous media flow using the pore-network modeling approach.
Journal of Computational Physics: X , 1:100011, 2019.[62] S. Whitaker. Flow in porous media i: A theoretical derivation of darcy’slaw.
Transport in porous media , 1(1):3–25, 1986.[63] M. Zamir and R. Budwig. Physics of pulsatile flow.