Branching exponents of synthetic vascular trees under different optimality principles

The branching behavior of vascular trees is often characterized using Murray's law. We investigate its validity using synthetic vascular trees generated under global optimization criteria. Our synthetic tree model does not incorporate Murray's law explicitly. Instead, we assume it holds implicitly and investigate the effects of different physical constraints and optimization goals on the branching exponent that is now allowed to vary locally. In particular, we include variable blood viscosity due to the F{\aa}hr{\ae}us--Lindqvist effect and enforce an equal pressure drop between inflow and the micro-circulation. Using our global optimization framework, we generate vascular trees with over one million terminal vessels and compare them against a detailed corrosion cast of the portal venous tree of a human liver. Murray's law is implicitly fulfilled when no additional constraints are enforced, indicating its validity in this setting. Variable blood viscosity or equal pressure drop leads to deviations from this optimum, but with the branching exponent inside the experimentally predicted range between 2.0 and 3.0. The validation against the corrosion cast shows good agreement from the portal vein down to the venules. Not enforcing Murray's law explicitly reduces the computational cost and increases the predictive capabilities of synthetic vascular trees. The ability to study optimal branching exponents across different scales can improve the functional assessment of organs.


Introduction
The cardiovascular system is responsible for transporting blood to and from all cells in the human body, leading to hierarchical networks of vessels, called vascular trees, inside each organ.According to Murray [1], this hierarchy obeys scaling relations based on the minimization of the total energy expenditure of the system.Many factors influence and constrain this minimization process, such as the type and shape of the organ supplied, the demand for the organ's cells, and the existence of vascular diseases.The goals and constraints guiding their structural development and influence on the vascular system have yet to be entirely understood, even though extensive work has been carried out for over a century.Thus the analysis of vascular diseases based on the anatomy and physiology of the vascular structure remains a challenge.
Murray first described a minimization problem for vascular segments in 1926 [1,2].Here, a tree is approximated as a bifurcating network consisting of rigid tubes, and the physical principles for fluid flow follow Poiseuille's law.The goal is to minimize the total power of the tree network with its minimum characterized by Murray's law.It describes the relationship of the radius of a parent vessel r 0 against the radii of its children's vessels (r 1 , r 2 ) as a power law with ( The branching exponent γ became an essential parameter for characterizing the branching behavior of vascular trees.In Murray's original formulation, γ = 3.0 is constant across the entire network.
An extensive number of studies have been conducted to investigate Murray's law experimentally [3][4][5].In general, exponents between 2.0 and 3.0 were measured.In [6], exponents were observed even going over the theoretical limit of 3.0 with γ = 3.2.Multiple theoretical studies have analyzed the possible factors contributing to these branching behaviors.An extension to Murray's law was proposed by Uylings [7], which incorporated the effects of turbulent flow into the minimization problem.Results show branching exponents as low as 2.33 for turbulent flow.In [8], a vascular model was investigated, which considered the role of elastic tubes.Compared to rigid tubes, the effect of pulsatile flow lowered the optimal value to 2.3.Zhou, Kassab, and Molloi [9,10] generalized Murray's law hypothesis to an entire coronary arterial tree by defining a vessel segment as a stem and the tree distal to the stem as a crown.They showed that γ deviates from 3.0 even for steady-state flow and depends on the ratio between metabolic demand and viscous power dissipation.
An alternative approach to investigate branching exponents is to construct vascular trees synthetically.The most well-known generation method here is constrained constructive optimization (CCO) [11].The local optimization approach is directly based on Murray's minimization principles and allows to investigate different, albeit constant, values for γ, like 2.55 [12] or 3.0 [13].Another approach, based on Simulated Annealing (SA), included the branching exponent as an optimization parameter [14].Results show that the vascular topology and the metabolic demand significantly influence the value of the branching exponent.
Recently, the authors extended the CCO approach to finding a synthetic tree optimal both in (global) geometry and topology [15].Finding the optimal geometry is cast into a nonlinear optimization problem (NLP), which allows the investigation of various possible goal functions and constraints.
In this paper, we utilize this flexibility and go beyond previous studies by allowing the branching exponent γ to vary locally.Furthermore, we include a blood viscosity law based on the Fåhraeus-Lindqvist effect [16] and enforce equal pressure drop to terminal vessels.The goal is to investigate the change in branching exponents under these influences.We start by introducing the relevant definitions and assumptions to generate synthetic trees.We then cast our goals and constraints into NLPs and introduce our optimization framework in more detail.Finally, we generate full portal venous trees of the human liver with up to one million terminal vessels and compare them against a vascular corrosion cast of a human liver [17,18].

Definitions and assumptions
We represent a vascular tree as a directed branching network T = (V, A) with nodes u ∈ V and segments a ∈ A. Each segment a = uv connects a proximal node x u with a distal node x v .It approximates a vessel as a rigid and straight cylindrical tube and is defined by its radius r a , length a = x u − x v , volumetric flow Q a and apparent blood viscosity η a .The distal nodes of terminal segments are terminal nodes (leaves) v ∈ L, and the proximal node of the (single) root segment is the root node x 0 .A synthetic vascular tree perfuses blood at a steady state from the root segment down to the terminal segments inside a given (non-convex) perfusion volume Ω ⊂ R 3 , schematically shown in Fig. 1 Figure 1: Schematic of a vascular tree and its relation to nodes and segments.Red circles denote a node, and white rectangles a segment.This tree has a given inflow Q 1 = Q perf and equal terminal outflow As in Murray's original paper [1], we assume laminar flow and approximate blood as an incompressible homogeneous Newtonian fluid.We express the hydrodynamic resistance R a of segment a by Poiseuille's law with The pressure drop over a segment can now be computed as and the pressure at a node v follows with We further assume that the (known) perfusion flow Q perf is homogeneously distributed among all N terminal segments, leading to a terminal flow value Q term = Q perf /N.All remaining flow values can then be computed using Kirchhoff's law with We aim at generating vessels down to the smallest arterioles/venules with typical radii in the range of 0.015 mm to 0.1 mm.The Fåhraeus-Lindqvist effect [19] should be accounted for at this scale.It describes how the blood viscosity decreases as the vessel diameter decreases.The tendency of red blood cells to migrate toward the vessel center is largely responsible for this effect.In turn, this forces plasma toward the walls and decreases peripheral friction.At the smallest vessels with radii approaching the radii of red blood cells, the viscosity sharply rises again.Pries et al. [16] derived an empirical relationship for this behavior with 2.44 exp(−8.08(ra /mm) 0.645 ) + 3.2, (7) where η p is the viscosity of the plasma, which we set to η p = 1.125 cP.η 45 is the relative apparent blood viscosity for a discharge hematocrit of 0.45.This relationship is depicted in Fig. 2 for the relevant radii between 0.015 mm and 10 mm.2.2 Design goals and constraints

Murray's minimization problem
The original minimization formulation by Murray states that the total power of a vascular tree consists of the metabolic power required to sustain blood P vol and the viscous power P vis required to pump blood from the root down to the micro-circulation.The cost function of a tree is then defined as where m b is the metabolic demand of blood in µW mm −3 .As described in [15], we can now include the nodal positions x, the length and radii r as well as the blood viscosity η in the vector of optimization variables y, leading to y 1 = (x, , r, η).We add physical lower bounds − , r − and η − and, for numerical efficiency, upper bounds + , r + and η + .The best geometry is then found in the rectangle defined as Our NLP "Power minimization" finally reads: min Eq. ( 12) fixes the position of terminal nodes and Eq. ( 13) ensures consistency between nodal positions and segment length.The third constraint in Eq. ( 14) enforces the Fåhraeus-Lindqvist effect as defined in Eq. ( 5).

Enforcing equal pressure drop
In Murray's original formulation, no consideration of the resulting pressure values at terminal segments was given.This terminal pressure is a crucial parameter for the regulation of blood flow and blood velocity at the microcirculatory domains.Since we assume these domains are roughly homogeneous across the organ and have equal demand, the pressure should not differ significantly.Therefore we enforce equal pressure at each terminal segment by adding the pressure p v at each node as a new unknown in our (second) NLP with variables y 2 = (x, , r, η, p), leading to Secondly, we constrain the pressure drop between the root and the terminal nodes as a prescribed constant value ∆ p .Since the viscous power at each segment a is directly proportional to the pressure drop by a factor of Q a , the total viscous power P vis becomes a constant.Thus, we can remove it from the cost function.Finally, we can drop the constant factor m b , leading to a minimization goal proportional to the tree volume V T .This formulation is used in most synthetic tree studies, e.g., [11,12,14,15,20].Our NLP "Volume minimization" then reads: Remark 1: Murray's law, as stated in Eq. (1), is not incorporated explicitly.Instead, we assume it holds implicitly and compute the associated branching exponent γ using the Newton-Raphson method after the optimization is finished.

Additional optimization variants
To better isolate the individual influence of different factors, we define additional variants of our two minimization problems.Firstly, we simplify both problems to a constant apparent viscosity η const = 3.6 cP, removing η from the vector of optimization variables and dropping the corresponding constraints Eq. ( 14) and Eq. ( 22).Secondly, we investigate the influence of the metabolic demand m b and the total pressure drop ∆ p .We consider values between 0.1 µW mm −3 and 1.0 µW mm −3 for the metabolic demand, but note that estimates of this parameter vary significantly [21].For the total pressure drop, we set the terminal pressure to p term = 6 mmHg and vary the root pressure between 10 mmHg and 14 mmHg.Finally, for each variant, a separate tree is generated, where Murray's law (see Eq. ( 1)) is enforced directly with a single exponent γ opt , included in the optimization variables.All variants include the same root flow Q perf = 1.1 l/min.

Generation framework
We generate our synthetic trees using the framework introduced in [15], which we summarize in the following.First, we generate N topo terminal nodes on a regular cubic grid inside our organ's volume.The root position is manually set and connected to the geometric center of the volume, which in turn is connected to all terminal nodes.We swap segments to explore new topologies from this initial (fan-shaped) tree.A swap detaches a segment from its parent and connects it with another segment.After each swap, the geometry is optimized by solving the corresponding NLP.The newly created topology is accepted on the basis of an SA approach.
After topology optimization, we grow the tree using a modified CCO approach.Here, we optimize the global geometry each time after adding N geo new terminal nodes and then increase N geo heuristically based on the current density of the tree.Notably, we drop the local optimization of branching positions and set them to their flow-weighted mean, similar to [22].Due to our repeated global geometry optimization, this simplification had no significant impact on the final tree structure.In the last step of the optimization, we delete all segments that reached the lower bound − (degenerate segments), possibly creating n-furcations (n ≥ 3).
We then classify the hierarchy throughout the finished tree by assigning each segment an order number corresponding to the Strahler ordering method [23].Continuous segments of the same Strahler order correspond to one vessel.Additionally, we employ an ordering scheme based on [4] to allow direct comparison (in reverse order) to the generation notation used for the vascular corrosion cast in [18].

Remark 2:
The complete optimization framework is only applied once to obtain a common topology.We solve each NLP variant with this topology to get the corresponding global geometry.We choose this method to focus on the geometry changes and to allow a direct comparison of branching exponents and radii at the same branch types.

Overall structure of the synthetic portal vein
The full synthetic portal vein tree for the case of volume minimization with variable viscosity is depicted in Fig. 3.The left side shows the complete tree inside the non-convex liver domain with two zoom levels.
The tree splits into four major branches, which further split into 8 main branches.These results align with previous results of a sparser tree in [15].
For a detailed comparison between the different optimization variants discussed in Section 2.2.3, we summarized the results in Table 1.Here, a column represents the results of a single variant, with the first four columns corresponding to the NLP "Power minimization" and the last four columns corresponding to the NLP "Volume minimization".For "Power minimization", we included the results for metabolic demands m b of 0.1 µW mm −3 and 1.0 µW mm −3 .Similarly, for "Volume minimization", we included the results for root pressures p root of 10 mmHg and 14 mmHg.The last row indicates the results of each variant after enforcing a single (optimal) branching exponent γ opt .
For "Power minimization", an increase in metabolic demand m b from 0.1 µW mm −3 to 1.0 µW mm −3 leads to an increase in viscous power P vis , shown in the first row of Table 1, by around 364 % and a reduction in volume V T , shown in the second row Table 1, by around 53 %.Similarly, for the "Volume minimization", an increase in the root pressure from 10 mmHg to 14 mmHg leads to a 320 % increase in viscous power and a 51 % reduction in volume.The Fåhraeus-Lindqvist effect had a minor influence.It decreased the viscous power by around 1 % and the total volume by around 4 % in all cases.

Vessel radii
The root r root is shown in the third row of Table 1.It decreased by 31 % after the metabolic demand m b was increased for "Power minimization".Similarly, in the case of "Volume minimization", it decreased by 30 % after the root pressure p root was increased to 14 mmHg.In all cases, the Fåhraeus-Lindqvist effect had a limited influence on the root radius with changes less than 0.1 %.However, a significant decrease in radius can be observed for vessels between Strahler order 1 and 6, shown in Fig. 4. In both NLP cases, this decrease was highest at the terminal vessels with 2 % for "Power minimization" (Fig. 4a) and 1.5 % for "Volume minimization" (Fig. 4b).A notable difference between both cases is the variance of radii at each Strahler order.In the case of "Power minimization", the highest variance is observed at Strahler order 6, whereas the terminal radii are constant.In the case of "Volume minimization", the highest variance is at the terminal segments and decreases as the Strahler order increases.

Pressure drop
After "Power minimization", the terminal pressures are not constant across the tree but exhibit a wide range of values, see Table 1 row 4 column 1-4.This range widens further for higher metabolic demands and also increases the mean total pressure drop from the root to the terminal segments, shown in Fig. 5.In Fig. 6a,  Including the Fåhraeus-Lindqvist effect leads to slightly higher pressure values for the Strahler orders 1 to 6.The terminal pressures after "Volume minimization" are fixed to p term = 6 mmHg as enforced by Eq. ( 19) -Eq.( 21).The effect of these constraints is highlighted in Fig. 6b for root pressure p root = 10 mmHg.In contrast to "Power minimization", variances are significantly higher at the intermediate Strahler orders 3 to 11.Furthermore, the influence of the Fåhraeus-Lindqvist effect is more pronounced, decreasing the pressure values between Strahler order 2 to 10.

Branching behavior
The resulting branching exponents of all variants are summarized in row 5 of Table 1.For "Power minimization" with constant viscosity η const , the exponents are constant with γ = 3.0 across all branches regardless of metabolic demand m b .In contrast, the inclusion of the Fåhraeus-Lindqvist effect leads to deviations from 3.0, with branching exponents reaching a minimum of 2.9.For "Volume minimization," exponents are not For a more detailed comparison, the probability density function of the branching exponents for both optimization cases is compared in Fig. 7.The influence of the Fåhraeus-Lindqvist effect shifts most branching exponents from a constant 3.0 to 2.9 during "Power minimization" (Fig. 7a).During "Volume minimization" with constant blood viscosity, most exponents are at 3.0 (Fig. 7b in red) and are shifted to 2.9 when including the Fåhraeus-Lindqvist effect (Fig. 7b in green).
Fig. 8 highlights the distribution of mean branching exponents across different branch types.Each cell (i, j) corresponds to a branch with child segments of Strahler order i and j.In Fig. 8(a) the effect of variable blood viscosity on "Power minimization" is depicted.The branching exponents decrease if the Strahler order of either child decreases, leading to the smallest branching exponent of 2.9 at branches with two terminal segments.If both child segments have a Strahler order over 8, the mean branching exponent is at its maximum of 3.0.The effect of enforcing equal terminal pressure is shown in Fig. 8(b).Here, the higher the difference between the Strahler orders of both children is, the smaller the branching exponent is.Again, the smallest mean branching exponents are observed at branches connecting two terminal segments with a value of 2.76.Fig. 8(c) shows the accumulated effect of both constraints with a minimum mean exponent of 2.7, again at terminal branches.

Comparison to vascular corrosion cast
We now directly compare our synthetic trees against a corrosion cast of the portal vein of the human liver [18].In Fig. 9a, the radii per generation (radius-adjusted Strahler order in reverse) are compared between the "Power minimization" with m b = 0.1 µW mm −3 , the volume minimization with p root = 10 mmHg and the corrosion cast data, including measurements and a best-fit trend line, based on the least sum of square errors.
The synthetic trees' radii fit the data and trend line well for the generations 5 to 15. Notably, they significantly deviate for the first four generations, especially against the measurements with errors of around 25%.In contrast, the number of vessels in Fig. 9b of both synthetic trees fit the data of the corrosion cast well for all generations.

Discussion
The exclusion of Murray's law simplifies our two optimization problems and automatically allows the branching exponents to vary locally.Both problems can generate synthetic trees with similar geometry and topology using appropriate parameters.Power minimization with a metabolic factor of m b = 0.1 µW mm −3 leads to more realistic radii than a value of 1.0 µW mm −3 , which is in line with estimated values around 0.3 − 0.4 µW mm −3 for venous trees [21,24].Similarly, a root pressure of p root = 14 mmHg results in a pressure drop of ∆ p = 8 mmHg, which is related to portal hypertension.Such a high pressure drop leads to unrealistic small radii in comparison to a more realistic root pressure of p root = 10 mmHg.
Under power minimization, the radii are at their individual local optima, i.e., the radius of each segment can be solved independently as the minimum of the metabolic demand and the viscous power dissipation.This observation is in line with the findings of Murray and explains the constant branching exponent of 3 in Table 1.The variations in radii for Strahler orders 2 to 12 in Fig. 4a are based entirely on the branching type, which is completely defined by the flow values of the child segments.
For volume minimization, no such simplification can be made, as the constraint of equal pressure creates dependencies between segments on the same path to the root.This constraint also forces radii to deviate from their local optima, which means a deviation from the branching exponent γ = 3.0.The highest deviations are found at branches between two terminal segments, because they can be adjusted to a given pressure drop without significantly increasing the tree's volume.Given the same length, they also constitute a higher pressure drop than segments with bigger radius.
The inclusion of the Fåhraeus-Lindqvist effect reduces the viscosity of the smaller vessels and, in turn, increases their pressure drop.The corresponding branching exponents also deviate from γ = 3.0, as can be observed in Fig. 8.In contrast, the effect on bigger vessels is negligible and results in constant exponents γ = 3.0, as seen in Fig. 8(a), for branches where both children have Strahler orders over 8.
While the generated trees generally fit the corrosion cast data well, the underestimation of the largest radii (generation 1 to 4) is significant.Since the vessels of the portal venous tree are more elliptical than circular, the radii of the vascular corrosion cast were estimated based on their cross-sectional area.Perhaps the synthetic trees' radii may correspond better with estimations based on other criteria such as a maximal inscribed circle.Another explanation is that the total energy dissipation is likely higher than our models predict.This underestimation could be due to ignoring the effect of turbulent flow [8] and simplified geometric modelling of the branching of vessels.The effect of pulsatile flow, however, can be neglected because the blood comes not directly from the heart but first flows through the digestive tract, leaving only a limited amount of pulsatility inside the portal venous tree.

Conclusion
Our optimization framework can handle complex constraints and goal functions while generating synthetic trees up to but not including the capillary level of the microcirculation.We used our framework to investigate the local branching behavior for different constraints and goal functions.Branching exponents automatically lie in the experimentally predicted range between 2.0 and 3.0.Even small changes to Murray's original optimization problem, like the inclusion of variable blood viscosity, significantly affect the optimal branching exponents of vessels.The topology and geometry of our synthetic trees closely follow the vascular corrosion cast of a human portal venous tree, with significant deviations only in the largest vessels.
Without enforcing any pressure constraint, terminal pressures vary by up to 1.0 mmHg, leading to highly heterogeneous boundary conditions for the microcirculation.In contrast, enforcing equal terminal pressure in its current form leads to pressure variations up to 2.0 mmHg in the intermediate vessels of the mesocirculation.Both these results need to be critically evaluated and compared against measurements of real vascular trees.One approach to improve the pressure constraint could be to prescribe a certain physical range using inequality constraints rather than a fixed value.In the future, we plan to include pulsatile and turbulent flow effects in a closed form, improving our framework's computational ability.Furthermore, shear stress, a critical parameter for vascular growth, must also be incorporated into the model.
A more mature version of this model, specifically the ability to predict optimal branching exponents under different constraints, could have many potential applications in the medical field.An example would be to predict and relate the branching behavior across the scales to vascular diseases.These predictions could improve the interpretation of medical images by giving valuable input to the functional assessment of organs.

Figure 2 :
Figure2: Change in apparent blood viscosity due to the Fåhraeus-Lindqvist effect as approximated by Pries et al.[16]

Figure 3 :
Figure 3: Complete synthetic vascular tree of the portal vein of a human liver with 1, 000, 000 terminal vessels (volume minimization with variable viscosity and p root = 10 mmHg).Two zoom levels highlight the hierarchical structure at different scales.The radii are between 5.1 mm (root vessel) and 0.017 mm (smallest terminal vessel)

Figure 4 :
Figure 4: Influence of Fåhraeus-Lindqvist effect on vessel radii for different Strahler orders

Figure 5 :Figure 6 :
Figure 5: Density plot of the total pressure drop from root vessel to terminal vessels for the power minimization under different metabolic demand

3 )Figure 7 :Figure 8 :
Figure 7: Effect of enforcing variable viscosity and equal pressure onto the branching exponents

Figure 9 :
Figure 9: Comparison of synthetic trees (with variable viscosity) against corrosion cast data [18]

Table 1 :
Comparison of the results between the different variants of our two minimization problems as introduced in Section 2.2.3.The branching exponents γ opt in the last row (marked with *) are the results of separate runs for each variant, where a single constant branching exponent was enforced.