A multi-scale ﬂow model for blood regulation in a realistic vascular system

Background: In the last decade, numerical models have been an increasingly important tool in medical science both for the fundamental understanding of the physiology of the human body as well as for diagnostics and personalized medicine. One challenging problem in this modelling is that the vascular systems are made of blood vessels at diﬀerent scales. Methods: In this paper, a multi-scale model is developed for blood ﬂow and regulation in a full vascular structure of an organ. In this model, a 1D vascular graph model that represents blood ﬂow in larger vessels is coupled with a porous media model that describes ﬂow in smaller vessels and capillary bed. The vascular model is based on Poiseuille’s law, with pressure correction by elasticity and pressure drop estimation at vessels junctions. The porous capillary bed is modelled as a two compartments domain (arterial and venal) and Darcy’s law. The ﬂuid exchange between the arterial and venal capillary bed compartments is deﬁned as blood perfusion. Results: The numerical experiments show that the proposed model for blood circulation: 1) is closely dependent on the structure and parameters of both the vascular vessels and of the capillary bed, and 2) it provides a realistic blood circulation in the organ. Conclusions: The advantage of the proposed model is that it is complex enough to capture the underlying physiology reliably, yet highly ﬂexible as it oﬀers the possibility of incorporating various local eﬀects. Furthermore, the numerical implementation of the model is straightforward and allows for simulations on a regular desktop computer. While the resistance of venous structures in the artery-1 occlusion increased to 0.73 and the resistance of venous structures in the arteries-2 occlusion decreased to 0.47 compared to non-occluded structure. It shows that the ﬂow in the venal structure is aﬀected by the


Background
In the last decades, numerical models and computational approaches have become important tools to study structure, function, and blood regulation of vascular systems [1,2,3,4,5,6]. The fundamental purpose of developing a numerical model is to understand how changes in the vascular structure affect transport mechanisms in organs, giving important clinical information. However, one of the challenging modelling issues for these systems is the fact that vascular systems are made of vessels at different scales [4], ranging from large arteries with turbulent blood flow, to smaller arteries and arterioles with mostly laminar flow, and to capillaries, with perfusion of cellular particles between artery and vein capillaries. This makes it difficult to fully understand the mechanisms between each scale and the effect of changes on a localized part or the whole organ.
Multi-scale mathematical modeling for blood circulation offers a reasonable solution to these difficulties. By combining well-established flow models at various vascular scale levels, one may achieve a global model based that also takes into account various local properties. The coupling between each level is the key to create a complete model with a suitable balance between computational cost and detail, as required for pathological characterization [4,3,7]. Passerini et al. [2] propose a 3D/1D coupling for combining large and small arteries in cerebral vasculature. In [6], the authors consider a network representation which includes arteries, veins and capillaries. This model assumes that the whole vasculature is an interconnected tube network and results in a much larger system. However, the computational cost of this approach limits its applicability at a clinically relevant scale. In a recent study [8], Hodneland et al. treat vessels as a 1D network model and the capillary bed as a 3D continuum model, thus decreasing significantly the computational cost and allowing for full brain simulations. It must be mentioned, however, that 1D models have a reduction in accuracy when compared to full 3D models, in particular when applied to describe large vessels [9].
In the current work, a multi-scale model is proposed for blood circulation that incorporates nonlinearities induced by the vascular structure. The blood flow in the vascular network (arteries and veins) is described using a vascular graph model [6] in which vessels are represented as long cylindrical tubes with constant radius. To compensate for the accuracy loss of this model, a model for pressure drop is included at vessel bifurcations [9,10] and allow for a vessel radius dependence on pressure due to vessel elasticity [11]. Darcy flow is described the continuum representation of the micro-circulation in the smaller vessels and the capillary bed. Further, a continuous distribution is defined to model the unresolved structure between the vascular network terminals and the Darcy domain. The purpose of the fluid distribution function is to reflect the micro-vessel structure surrounding each terminal node. Differently from [8], we look at the impact of nonlinear interactions between pressure and flow, as well as the role of vessel bifurcations in the flow network. Our model was developed with a minimum number of boundary conditions and unknown parameter settings to obtain result purely from governing equations and the given anatomy.

Materials and Methods
Graph structure model for vascular networks A system of equations was constructed based on a network structure for both arteries and veins. Assuming laminar flow and non-slip conditions on the vessel walls, each vessel n was described as a long cylindrical tube of length L n with constant radius r n L n . The vessel radius was computed as an average radius value from the segmentation data. The pressure drop ∆P n in a single vessel segment n was computed using Hagen-Poiseuille's law [6] ∆P h n = 8µL n q n πr 4 where the upper index h stands for hidrodynamic, µ is the blood viscosity and q n is a volumetric blood flow. At a junction node, an additional pressure drop estimation was defined from [10,9], where the upper index b stands for bifurcation and the index dat refer to the datum vessel, i.e. the vessel from which the flow approaches the junction. Further, θ (dat,n) is the angle between the datum vessel and vessel n. The pressure drop equation was constructed based on Mynard et. al.'s work for converging and diverging flow at a junction [10]. Hence, the total pressure drop after a bifurcation node was computed as the sum of both equation (1) and (2), The other governing equation was the conservation of mass at a node with q in representing the blood that flows into the node and q out is a flow out of the node.

Elasticity on vessel walls
Vessel wall elasticity is a physiological regulatory factor that forces blood flow in a particular direction. On the other hand, the vessel wall elasticity also allows for a vessel radius dependence on the pressure gradient, to favour blood supply in case of an alteration of the vascular system [11]. The change in radius due to pressure was described by where r 0 is the initial vessel radius without pressure gradient, h is the vessel wall thickness, E is the Young modulus, λ is Poisson ratio, and P in , P ext are the pressures inside and outside the vessel, respectively. The vessel thickness was assumed to be proportional to the vessel radius, with h = 0.2r for the arteries and h = 0.1r for the veins [11]. For the elasticity parameters, the Young modulus in the arteries E a is two times bigger than E v in the veins, E a = 2E v = 0.5 MPa, and λ = 0.5. P in is defined as an average pressure in the segment, that is the average pressures at both segment endpoints, and P ext as the capillary pressure at the midpoint outside the vessel. The elasticity equation (5) had the effect of introducing nonlinearity to the system. Thus, the pressure difference between two adjacent nodes in a vessel was determined by substituting equation (5) into equations (1) and (2) Capillary model Under resolved, smaller vessels and the capillary bed were discretized with a uniform grid and described by Darcy's single-phase flow equation. Darcy's law, which is describing the filtration of a fluid in a porous medium, states that a fluid flows from regions of higher pressure to regions of lower pressure. Thus where v is the Darcy flux (volumetric flow rate per unit area [m 3 s −1 m −2 ]), K is the permeability tensor of the porous medium and µ is the viscosity [12]. In addition, we assume conservation of mass (continuity equation), where Φ is the porosity of the capillary bed [-] and Q is the source term [s −1 ]. In this model, Q is a source or sink, for instance describing the flow in or out a terminal node arteries or veins, or a describing the passage of flow from one compartment to the other. We assumed blood to be an incompressible fluid, and by incorporating Darcy flow into the continuity equation, we obtain The capillary bed was described as a two-compartment system, using one compartment for the arterial part and one for the venous part. Blood perfusion is the interchange between the two compartments, i.e. the exchange of oxygenated blood with deoxygenated blood in the capillary bed. The driving force for perfusion is the pressure difference between the two compartments [8].
where α is the perfusion parameter. The Darcy model for the two compartments becomes where the index β = {a, v} stands for artery and vein respectively.

Coupling vascular structure and capillary models
To complete the system, both the arterial and venal Darcy systems were combined with the terminal nodes of the arterial and venal networks as point sources/sinks. The analytic solution of the Darcy equation around these points is a Dirac delta function, with a singularity at the source/sink point, creating a problem at the 3D-1D interface. To avoid this problem, a fluid distribution function was introduced [8] and use this function with a finite radius , f (x) = f ( x ) on the Darcy domain. Thus, at a terminal node l (source/sink), we take Beside mass continuity, we had to describe pressure continuity between the 1D vascular graph and the 3D Darcy model. The pressure drop between the terminal node i and the surrounding tissue was thus given as where κ represents the resistance estimation for capillary system around terminal i in the Darcy model. Both flow and pressure couplings between the vascular graph and Darcy model close the loop from the arterial roots to the venal roots.

Time of flight
Time-of-flight (TOF) was used as a parameter to measure the distribution of blood flow or nutrients in a domain. Time of flight is a total time (T ) needed for a passive particle to move along a streamline Ψ from a point source to a certain location (x) in the domain. It is usually defined as where s is the path-length and Φ is the porosity of the medium (in our case, the capillary bed). In this paper, we assume Φ = 1. If we apply a directional derivative of T (x) with respect to the distance along a streamline in the velocity direction, we obtain The equation (16) is solved with Dirichlet boundary conditions at the point source: where T is the TOF in the computational domain Ω, T N is the TOF in the arterial terminal nodes, and ∂Ω − is a point source. From the Darcy flux, we obtain v as the total flux velocity in each computational cell. T N is obtained from solving the TOF in the arterial network. a similar system of equations was constructed with the capillary bed based on the flow rate from the graph vascular model, where q is a 1D vector volumetric flow rate, V /L is volume per length of the related segment vessel, and T = 0 on ∂Ω − is the boundary condition on the arterial root terminals.
The solution of this system is obtained efficiently by arranging equations according to their one-sided directional dependency. Both the analytic and the numerical solutions at any point in the domain are obtained from the solution of the point upstream. A pointwise solution was constructed, starting from the source (inflow boundary condition) and continuing downstream to the next points to cover the entire domain (Fast Marching algorithm [13]).
To simplify further, it focus on the volumetric flow rate, rather than the flux velocity. The equation (17) is modified by multiplying both sides with the cell face area, q ij = v ij · A ij and using the the cartesian grid of the Darcy discretization, we obtain where N i is the the set of upstream cells for cell i, q ij is the volumetric flux at the face between cells i and j, V i is the volume of cell i, and d ij is the distance between cells i and j. An upstream cell was defined to be any cell j adjacent to cell i that gives influx to cell i.

Numerical implementation of the flow
The vascular network The vascular network, consisting of segments and nodes, is governed by (1), (4) and (11). For each graph segment i, the pressure drop (1) (linear part) and the pressure drop at junctions (2) (nonlinear part) were considered. When considering elasticity, equations (5) was used.
The equations at the nodes are three different types: at the root nodes, it is control equations/boundary conditions (given with the problem); at the internal nodes (connecting two or more segments) it is conservation of volume (4), and at the terminal nodes we have coupling equation with the Darcy model. In our four roots terminal (two arterial and two venal roots), there were several options to assign the boundary conditions.
Further five control equations must be considered for the system to have a unique solution. These were constructed by assigning constant pressure on the arterial and the venal root nodes (Dirichlet BC) and conservation of mass for the whole system.
where root, β = a, v is a root node on the arterial or the venal structure, P ext represent the constant pressure on two arterial roots (input) and two venal roots (output), and N root are root segments, with q i > 0 for the artery, q i < 0 for the vein.

The Darcy domains (capillary bed)
The Darcy equations (10) in the capillary bed were solved using a two-point flux approximation (TPFA). First, we computed the fluid transmissibility between adjacent cells i and j, where S ij is the face area between cells i and j, and, as above β ∈ {a, v}. Because of the homogeneity of the system, constant parameters and uniform discretization, equation (22) simplifies to Then, applying TPFA to (6) in single cell i and adjacent cells j ∈ N i , neighbour cells around i, we obtain Applying this procedure for all cells in the domain, we obtain a system of equations We note that while the TPFA method is not a consistent discretization for general grids, it is consistent and convergent on the quadrilateral grids used in this study [14].
The full coupling Finally, the vascular networks and the Darcy domains were combined in a (nonlinear) system of equations Ax = b, with unknown x consisting of the pressure and the flow rate in the model, The indexes N and T refer to the internal and terminal nodes in the vascular network, and index D stands for the Darcy equation discretization. The unknown solution x N is the pressure at the internal nodes and flow rate in the corresponding segment, x T is the pressure in the terminal nodes, and x D is the pressure on the Darcy domain. The blocks matrices represent systems of equations in the following domains: • Matrix A N −N ∈ R N ×N is the system of nonlinear equations that describe the vascular graph model. It is based on equation (3), (4), and (5) and consists of as many equations as the number of vessels and nodes in the internal nodes network. It also has no more than four non-zeroes entries per row. • Matrix A T −T ∈ R T ×T is the system of nonlinear equations on terminal nodes at the coupling between the vascular graph and the capillary model. It is based on equation (3) and (5). The matrix is a diagonal matrix.
The above nonlinear system (27) is solved using the solution of the linearized system as initial approximation. Frog tongue image data The model described above is generic and applicable to any 2D or 3D domain.
As an illustrative example, we used the model to simulate blood flow in a frog tongue from a 2D anatomical image from a classical biology textbook [15]. This example resembles a realistic vascular system while still being more easy to visualize, interpret, and analyze compared to a full 3D system. A drawback of this choice is that the original sample is not available, thus we cannot validate the correctness of the vessel structure, as well as there is a deformation in geometry. Incorrect vessel structure and incorrect geometry may have consequences for the blood flow in the entire domain, but this is outside the scope of the present paper.  Figure 2 shows a frog tongue, with a physiologist's segmentation of the main vessels. The frog tongue was stretched flat and nailed to a canvas. Thereafter, a trained physiologist manually traced the arterial and venal network structures on it. Only a few vessels serve the middle part of the tongue.
Both the arterial and venal network structures consist of two main vessels (left and right). Our network structure was based on an image segmentation of these using the method in [16]. The small arteries and veins that are near the roots but not visibly connected to the main vessels were discarded in our segmentation, as we need control over inlets and outlets, so that only the principal ones was kept. We observed that both the arterial and venal networks have some level of anastomoses. The two arterial networks are connected in uppermost top of the tongue (point D in figure 2) by a small vessel and the venal networks are connected in the middle of tongue by a large vessel. These connected structures have a regulation effect, maintaining blood circulation for the whole organ when some vessel is severed or occluded.
The frog tongue model parameters for our simulations are found in table 1.

Result and discussion
Using the model proposed above model, we have simulated blood flow in several numerical experiments: a baseline model simulation, linearized simulations, and simulations with vessel occlusions. In particular, the occlusion experiments demonstrate possible applications to the study of blood circulation in case of alteration or pathology of the vascular system.

Baseline model
The baseline model is the fully nonlinear model with the nonlinear equation due to pressure drop at junctions and vessel elasticity. It was based on equations (1), (2), (4), (10), (5), (12), and (13). Figure 3 shows the computed pressure distribution in the whole domain. We observed that the pressure is higher around the top edges and is lower in the middle of the domain. This result is related to the vascular structure having a small number of terminals in the central region.
The total pressure drop in the whole system was 6 kPa (see table 1), the arterial network being the major contributor (60.96 % of total), arterial compartment 6.45 %, venal compartment 13.29 %, and the vein vascular network 1.16 kPa pressure drop (19.30 %). These values indicated the importance of the arterial vascular structure to provide blood circulation through the whole organ. If there is an alteration on the arterial vessel, its impact for blood circulation will be greater than for a similar vein alteration. Figure 4 shows TOF distribution in the Darcy compartments for the arterial and venal compartment. These value represents the time needed for blood to fill a certain location. The TOF in both compartments shows a stagnant point in the mid-bottom region of the domain (yellow area, see Fig. 4). This artefact is likely due to the model using a wrong geometry (mis-segmentation, ignoring some small vessels in the region) rather than model errors. In turn, observing such anomalies can also be used as a feedback on the quality of the segmentation/geometry.

Linearized simulation
The linearized simulation is the linearized baseline model by omitting nonlinearity due to vessel elasticity or pressure drop at junctions, or both from the system.
The vessel elasticity has a similar effect as a capacitor in an electronic system. Elasticity can maintain pressure constant by increasing or decreasing vessel radii. However, this effect also gives additional resistance to blood flow, equivalent to impedance in electricity. The arterial vessel tends to expand if the pressure in the surrounding Darcy domain is lower than vessel inner pressure. On the other hand, the vein vessel tends to tighten and increase resistance. This model gives a total resistance of 2.10 (see table 2), which is almost similar to the baseline model (The difference was not significant and unseen after rounding up the value). The arterial structure gave 60.94% of the total resistance (0.02% lower than the baseline) and venal structure 19.31% (0.01% higher than the baseline) with Darcy compartments representing only 6.45% and 13.30% of the total value. The elasticity model provided The pressure drop at the junctions increased the total pressure drop in the vascular structure, both in arteries and veins. It is also causing the total resistance in the whole system to increase: as the volumetric blood flow Q changes according to the relation ∆P = RQ, the pressure drop must increase proportionally to the resistance, in order to maintain a constant total flow.
The resistance in the arterial network still had a dominant value 61.07% in the junction model (0.11% higher than the baseline model). However, the vein resistance percentage was lower compared to the baseline model (0.06% lower). The compartments had 6.44% and 13.25% of total resistance. This results were more significantly changed compared to the elasticity model (see table 2).
The junction and elasticity model only incorporates one of the non-linear equations into the system. Both result shows a value almost similar to the simulation in baseline model. These results indicate that the effect of the junction pressure drop and the vessels elasticity do not stack up when combined. In the current experimental setup, the junction pressure drop hence becomes insignificant to the system. Its difference compared to the baseline model was insignificant (less than 0.2 %) when it was omitted from the model.
The linear model is the fully linear model based on equations (1), (4), (10), (12), and (13). The total resistance in this model was significantly lower compared to the previous models which allowed a bigger volumetric blood flow (see table 2). The arterial structure provided 73.4% of the total resistance and venal structure 19.95% with Darcy compartments representing only 3.21% and 3.39% of the total value which was a notable decreasing compared to the other linearized models. These The areas close to terminal nodes show a lower value(well-perfused areas) and the yellow pattern indicate the less oxygenated areas which are farther from terminal nodes. The pattern in the venal compartment is similar to the pattern in the arterial, but with more blurred areas of high TOF. The yellow pattern is likely due to mis-segmentation (leaving out some small input vessels in the region) rather than model error.
results indicated the removal of all nonlinearities significantly changed both the total resistance and resistance distribution.

Occlusion in root vessels
To understand the blood flow regulation due to pathology or alteration, a series of simulations were performed by occluding one root vessel in each network at a time. The occlusions factors for the vessel radii range between 5% and 95% of the original value. Occlusions of artery-1 and artery-2 gives an expected result. The pressure in the occluded part of the structure becomes significantly lower than in other part. Because of an asymmetric structure of left and right of arterial network, the pressure distributions are not mirroring each other. Furthermore, the artery-1 occlusion cause a steeper pressure gradient (left to right) compared to the artery-2 occlusion (see Figure 6). Table 2 reports a summary of the occlusion experiments. Note that occlusions in the arterial roots increase resistance in the arterial structure and decrease total flow into the system. On the other hand, we have to increase the driving pressure around 1.5 times compared to non-occlusion model. Even if the total resistance in the artery-1 occlusion experiment is larger than for the artery-2 occlusion, the arterial structure resistance in the artery-1 occlusion take smaller part compare to artery-2 (the actual value are almost similar, 2.34 and 2.39). While the resistance of venous structures in the artery-1 occlusion increased to 0.73 and the resistance of venous structures in the arteries-2 occlusion decreased to 0.47 compared to nonoccluded structure. It shows that the flow in the venal structure is affected by the Figure 5: Average system resistance R as a function of the occlusion rate in roots. The occlusion rate (x axis) in root radius varies from 5% to 95% of the original value for each simulation. Artery-1 and 2 refer to the right and left roots of the arterial structure. Vein-1 and 2 refer to the right and left roots of the venal structure. Occlusions of venal roots result in small changes of the macroscopic flow parameter in the whole system. In particular for vein-2, the total resistance is close to being constant(∆R < 1%). arterial structure. However, the compartments resistance for both occlusions tends to be constant. The alteration on artery change the resistance in the whole system.
The numerical results for occluded venal networks show that the resistance change is less than 12 % of the non-occluded resistance. In particular, vein-2 occlusion doesn't give any significant change to the total resistance (∆R < 1%R non-ocluded ), even with the occlusion 95 % of the original value. The structure of the venal network has a good collateral circulation and maintain blood circulation in the whole system. Even if the macroscopic flow remain similar, the occlusion in the venal structure changes microscopic pressure inside the system. The venal compartment resistance in the vein-1 occlusion experiment has twice the value of the baseline simulation (see Table 2). The resulting pressure distribution change in the venal compartment allows for blood going through the veins more than once and thereby also cause a reduced blood flow in capillary bed. Figure 6c shows an example of blood flow in reverse direction due to high pressure difference in the venal Darcy compartment. Higher pressure gradients between several terminal nodes in the bottom left corner makes blood flow through venal vessels rather than in capillary bed. These phenomena and the other pathological effects are partly provoked by the severe occlusion in the experimental setup and partly caused by a mismatch between the observed macro-anatomical structures and the assumptions about the micro-structures in the Darcy domain. With further prior knowledge about local variability i.e. the capillarity density of the tissue, the current constant permeability assumptions are easily exchanged by a heterogeneous permeability field.

Conclusion
A multi-scale flow model has been presented to simulate blood circulation in a vascular system. The model is a coupling between a Poiseuille flow in the macroscopic vascular networks and a Darcy flow for micro-vessels and the capillary bed. It is coupled by introducing a mollifier around terminal nodes to model unresolved micro-vessels and the coupling with the capillary bed. The model is highly flexible and allows for easy incorporation of other local effects. As an example, it is able to incorporate a pressure drop at the vessel junction nodes and vessel elasticity, allowing radius dependence on the pressure. In the numerical experiments, we simulate several instances of blood flow in a frog tongue anatomical image. In all the experiments, our model gives overall realistic simulation results. Anomalies in the results (like stagnation points) can be traced back to incomplete or erroneous data or geometry, thus also providing a feedback to quality of the geometry of the problem and the initial data. The linearized simulations indicate the nonlinearities have significant impacts to the system and the effects do not stack up when combined. The nonlinearity due to pressure drop at junctions is insignificant when the elasticity is applied to the system hence it can be ignored.
Vessel networks alteration simulations confirm that both local and global blood circulations in the vascular system depend on the vascular structure. In our simulations the existence of an anastomosis in form of a big connected vessel in venal network system has an important role in maintaining blood circulation when a main root is occluded. The occlusion in the venal structure does not give significant change in the macroscopic flow. At the microscopic level, however, the pressure drop in the venal compartment increases twice as much compared to the baseline model. This indicates that it is important to simulate the complete vascular system including venal system, and not just arterial system.