Asymmetries arising from the space-filling nature of vascular networks

Cardiovascular networks span the body by branching across many generations of vessels. The resulting structure delivers blood over long distances to supply all cells with oxygen via the relatively short-range process of diffusion at the capillary level. The structural features of the network that accomplish this density and ubiquity of capillaries are often called space-filling. There are multiple strategies to fill a space, but some strategies do not lead to biologically adaptive structures by requiring too much construction material or space, delivering resources too slowly, or using too much power to move blood through the system. We empirically measure the structure of real networks (18 humans and 1 mouse) and compare these observations with predictions of model networks that are space-filling and constrained by a few guiding biological principles. We devise a numerical method that enables the investigation of space-filling strategies and determination of which biological principles influence network structure. Optimization for only a single principle creates unrealistic networks that represent an extreme limit of the possible structures that could be observed in nature. We first study these extreme limits for two competing principles, minimal total material and minimal path lengths. We combine these two principles and enforce various thresholds for balance in the network hierarchy, which provides a novel approach that highlights the trade-offs faced by biological networks and yields predictions that better match our empirical data.


I. INTRODUCTION
The vital functions of the cardiovascular system are the distribution of oxygen and nutrient resources throughout the body, as well as the collection and filtration of waste by circulating blood.Transfer of resources and waste occurs primarily at the capillary level via diffusion through nearby tissue.Consequently, this smallest level of the network must reach all living cells in order to maintain them, filling the entire space of the body.In models developed by Krogh for effective diffusion of oxygen [1], cells cannot survive beyond a maximum distance from a capillary.This defines a service volume of cells that is associated with each capillary, which has a typical size that has been observed to vary across species based on cellular metabolic rate [2,3].The constraint on maximum distance from capillaries necessitates that the final levels of the cardiovascular network are also space-filling throughout the body.In this paper we investigate the relation between this space-filling property and basic optimization principles such as the minimization of costs from construction material and pumping power.Specifically, we highlight how this relation influences the asymmetries in sizes and flow rates of sibling segments as measured in empirical data.
A central focus of our investigation of cardiovascular systems is the space-filling properties of networks, but these properties are also of great interest in many other contexts.General space-covering hexagonal patterns appear in nature in the cell structure of beehives as well as in economic theories for market areas [4].Trees (the woody, perennial plants) have been studied for both how forests fill an area [5], as well as how the vascular structure within an individual plant determines the hydraulics of resource delivery [6,7].Apollonian networks [8] describe the space-filling packing of spheres of various sizes, which are similar in the cardiovascular system to considering the volumes of different subtrees of the network.Efficiently filling space in two dimensions is important for information visualization [9].In addition to these applications, Kuffner and LaValle study space-filling tree networks (i.e., networks that branch and have no closed loops) to determine a route from one location to another [10], but without the biological constraints that we impose here.For cardiovascular networks, this motion planning is analogous to how the vascular structure routes blood.Efficient routing of blood through a hierarchy is central to models that investigate allometric scaling of metabolic rate with body mass [2,3,[11][12][13][14], which build on metabolic scaling by Kleiber [15].Determining specific space-filling strategies will inform these models to better describe the cardiovascular system.
Developmental processes (i.e., growth) as well as evolutionary pressures, such as efficiency in material and energy use, shape the structure of cardiovascular networks.Filling a volume or surface efficiently with the terminal nodes of a branching network is nontrivial, especially when the distribution system must reliably deliver blood at each stage of development.The system must also be robust to changes in vessel lengths and the number of hierarchical levels that result from growth and obstructions from damage [16,17].We propound that these structural challenges lead to the asymmetric features in both the local relations of segments, as well as in the global paths from the heart to each service volume that we observe in empirical data.
At the local level, the ratio of lengths between parent and child segments may vary across the network, deviating from strict self-similarity.Additionally, sibling segments may vary in length, which we call asymmetric branching.Within our numerical method, these features arise as a result of optimizing branch point positions relative to adjacent branch points across the network.Variation in the length and number of levels in paths means that the network is also not symmetric or balanced in these global properties.By allowing these asymmetries and explicitly ensuring space-filling structure, we expand other models that are strictly balanced in network hierarchy and perfectly symmetrical in local quantities.
Asymmetries observed in real systems motivate our investigation of the space-filling properties and asymmetries in cardiovascular networks.These observations show that such networks have a tendency for the lengths of sibling segments to be distributed less symmetrically than is the case for radii.The empirical data in Sec.V come from two types of sources.Images collected through microtomography of mouse lung comprise one data set.The mice were part of a study with different manipulations of matrix GLA protein (which causes the vasculature to be under-or over-branched [18]), but we focus on the wild-type specimen for our analysis.The other data set, collected through MRI, excludes the lungs and instead focuses on the central vasculature in the human head and torso [19].We utilize the custom software Angicart [19] to analyze these two distinct vascular data sets.Because of the noninvasive nature of these data acquisition and analysis techniques, future studies have the opportunity to track the development of cardiovascular systems as individual organisms grow and age, including repair after the system incurs damage (i.e., wound healing).Such studies can elucidate the effects that patterns of growth and changes from damage have on the final, mature state of the network.
In this paper, we study the optimization principles that correspond to evolutionary pressures for efficiency in material cost and power loss.Our focus is the influence of space-filling patterns on length asymmetry distributions without the explicit inclusion of radius information.The list of candidate networks includes all possible hierarchical (topological) connections between the heart and all capillaries.For each hierarchy and unique permutation of pairings between terminal vessels and service volumes (see Fig. 1d), we must determine the positions of the branch points.The combination of the hierarchy, service-volume pairings and branch point positions defines the configuration of a candidate network.For these reasons, we must search through many candidate configurations to determine the most efficient structure.We quantify the fitness of each candidate network using individual segment lengths between branch points as well as full path lengths between each capillary and the heart.
To perform a reliable comparison between candidate configurations, it is crucial to determine branch point positions in a consistent way.We determine these positions iteratively for the entire network in order to identify the global optimum.While the local process of choosing branch points that minimize total vessel lengths (or sim-ilar features) is relatively straightforward to iterate over the network, any single branch point and its relation to its neighbors relies indirectly on updates that are applied elsewhere.This dependence emerges from the fact that each end of a vessel is connected to a branch point, which upon repositioning affects the lengths of all vessels that it joins.The uniqueness of the Fermat point -the branch point position that minimizes the sum of the lengths of vessels at a single junction -is already well established (for example, see [20]).This allows us to carefully construct an algorithm (described in Sec.II C) that reliably relaxes all branch point positions into the global optimum.
After determining the positions of branch points for a given hierarchy, we compare distinct configurations to find the optimum network.The search through configurations is also a central problem in phylogenetics, where the goal is to construct phylogenies to identify similar groups of species and trace the development of genes through speciation.Even in the case of genes that control biological traits, a loosely analogous space-filling phenomenon emerges in the form of species filling the niches in the environment.With our specific goal of complete spatial covering of network tips, we develop strategies in Sec.IV for exploring the space of hierarchies that are similar to those used on phylogenetic trees.
The organization of the subsequent sections is as follows.In Sec.II, we describe the basic assumptions for our space-filling network model, including the details of the local optimization of branch point locations and the global paths through the network.In Sec.III, we introduce the specific quantitative network properties that we use to compare the fitness of candidate networks.We introduce the properties of the space of tree hierarchies and our implemented exploration strategies in Sec.IV.In Sec.V, we detail the results from the several layers of optimization that we implement, and discuss the insights that they offer in Sec.VI.

II. CONSTRUCTION OF ARTIFICIAL VASCULAR NETWORKS
To better understand the connection between the local asymmetries of individual vessel lengths and the global constraint on space-filling capillaries, we optimize candidate networks that are embedded in two spatial dimensions (2-D) in silico with respect to specific optimization principles.We explore these optimized artificial networks and quantify their branching length asymmetries to compare with our empirical data.Our model's simplification of the cardiovascular network focuses on the lengths of segments as defined by the straight line between adjacent branch points.Reticulated structures occur within leaves to mitigate damage [16,17] and within animals as anastomoses (or pathologically as fistulas).However, we focus on the vast majority of the cardiovascular system that distributes resources through a hierarchical, tree-like structure, in which no segments or subtrees rejoin downstream to form closed loops before reaching the capillaries.This is sufficient for the focus of our investigation of the asymmetric, space-filling structure that distributes the resource-rich blood from the heart throughout the body.
The space-filling property of the cardiovascular network constrains the hierarchical structure of the network and the positions of branch points.Here we describe our process for the construction of individual networks and the space of possible networks through the following steps: defining a distribution of space-filling service volumes in the space of the body, identifying all unique hierarchies and pairings between tips of hierarchies and distinct service volumes, and determining the positions of branch points for each hierarchy and pairing.

A. Space-filling service volumes
Because real systems do not organize or grow on a regular (symmetric, isotropic) grid, we position service volumes randomly within the space they fill.Construction of service volumes begins by choosing a random point within the body volume that represents the location of a capillary.We then randomly choose other points (capillary locations) so that none lie within a predefined constant distance from another capillary location.After determining a set of capillary locations that span the 2-D area, the entire body is partitioned into Voronoi cells fed by the closest capillary.In this way, each capillary becomes associated with a specific service volume, and the sum of the service volumes fills the entire space (see Fig. 3 or 5).

B. Space of hierarchically distinct trees and pairings with service volumes
Because multiple branching levels connect the service volumes to the heart, there are many possible hierarchical orderings of branching junctions across these different levels.For example, there are two unique hierarchies when there are four service volumes: the top three configurations in Fig. 1d are the same perfectly balanced hierarchy, while the remaining trees have the same unbalanced hierarchy.The distinguishing feature is the pairings of tips in the hierarchy to specific service volumes (1-4).
There are many distinct pairings of terminal tips with the set of fixed service volumes for each bifurcating tree.For four service volumes, there are a manageable 15 unique bifurcating trees (shown in Fig. 1d).Before determining branch point positions for small networks, we exhaustively search through all possible hierarchies and pairings with service volumes.For large networks, the number of distinct hierarchies and pairings ((2n − 3)!! for n distinct service volumes) is prohibitively large, so we sample the space as described in Sec.V B.
We do not disqualify configurations if one vessel path crosses with another (these would likely separate in three dimensions), and there is no exchange of resources or interaction in blood flow at such locations.Crossings are not observed for networks that minimize only total network length without a constraint on hierarchical balance, but they often occur for optimal configurations with a strong constraint on hierarchical balance.

C. Optimization of branch point positions for a fixed hierarchy and pairing
We now detail our algorithm for the optimization of the positions of branch points that connect a distribution of service volumes to the heart through a hierarchical network.Within our algorithm, the position of each branch point depends solely on the location of the adjacent branch points in the network.Distant vessels affect each other indirectly, but not through any direct longrange process.Using the limited, local information given by the neighborhood of a branch point, each junction is assigned a uniquely-defined position that minimizes the sum of the Euclidean distances to each neighboring junction.This is equivalent to the Fermat point of the triangle formed by the two downstream ends of the child vessels and the one upstream end of the parent vessel (see Fig. 1b).The Fermat point of a triangle is a special case of the more general geometric median, the unique point that minimizes the sum of distances to an arbitrary number of other fixed points.We follow the algorithm presented in [21] to avoid errors in determining the geometric median.Assigning branch point positions as geometric medians effectively minimizes the construction costs for the local network structure.
We construct our networks from simple bifurcations, but using the Fermat point to assign branch point positions can lead to coincident (degenerate) bifurcations, as shown in Fig. 1c.Degenerate branch points are consolidated at the geometric median of the upstream endpoint of the parent and three or more downstream endpoints of the associated children segments.In this way, two degenerate bifurcations become a trifurcation and, more generally, n degenerate bifurcations become a single (n + 1)-furcation.Networks that are hierarchically distinct in their bifurcating structure can become identical networks by collapsing bifurcations.Through exhaustive explorations (described in Sec.V A), we find that this marginally reduces the number of possible configurations (see results in Fig. 12a).Because we have no a priori filter to identify which bifurcating trees are redundant, we must consider the entirety of the space of labeled, rooted, bifurcating trees throughout the algorithm to identify sufficiently optimal configurations.With positions defined in a consistent way, we can now compare the properties of distinct hierarchies to determine which is the best for a particular space-filling strategy.

III. SELECTION CRITERIA FOR BIOLOGICAL NETWORKS
All characteristics of an organism that affect fitness and are heritable are under selection.A key question is which features of the vascular network are under selection.Here we define specific fitness measures that are tied to the structure of the network configuration that allows us to rank candidate networks and determine the optimal configuration.

A. Global length properties of space-filling configurations
Here we introduce a standardized fitness measure that allows us to compare candidate networks for their suitability to transport blood and resources.Each indepen-dent measure for a network's fitness relates to a physical quantity that likely guides the evolution of the cardiovascular system toward a more efficient network.Specifically, the system's cost in construction material and the maintained volume of the blood relates to the total network length -the sum of the lengths of all segments.Competing with this minimization of materials, the dissipation of the heart's pumping power relates to the path lengths between each capillary and the heart.The power dissipated by smooth (Poiseuille) flow through a segment is directly proportional to the length of the segment [22].
In the absence of radius information, reducing the cost of pumping blood is equivalent to reducing the total path lengths that blood travels.
We define these two fitness measures -one dealing with total network length and the other with individual path lengths -as where N tips is the number of distinct service volumes, corresponding to the number of tips and distinct paths.The generalized total cost function C linearly combines these two measures by their respective weights C L and C H : This cost function connects minimization of material and power dissipation to study optimal networks that are space-filling.Because an increase in cost corresponds to a decrease in fitness, we place this approach in an evolutionary framework by inverting and scaling our cost function C to be a relative fitness function F : where C min is the most optimal network under consideration.By this definition, the optimal configuration has a fitness F = 1 and less optimal configurations have a fitness F < 1.

B. Equal Distribution of Resources through Hierarchical Balance
Because the network tends to exhibit nearly symmetric branching in radius and must distribute resources equally to each capillary in the body, the network hierarchy cannot be overly unbalanced, with one segment having many more tips to supply than its sibling.In accordance with this argument, empirical data do not show major arteries branching directly into capillaries.We address this constraint by selecting for networks with more balanced hierarchies.
A hierarchy is better balanced if there are roughly equal numbers of tips supplied downstream by each sibling segment.Conversely, a hierarchy becomes more poorly balanced as the disparity grows between the number of tips.In this sense, we define the degree that a hierarchy is unbalanced U as where tips is the number of distinct downstream service volumes supplied through segment i and segment j is always the sibling with the most downstream tips.In our algorithm, we select against hierarchies for which the degree of unbalance U is greater than some threshold U 0 .We eliminate configurations above this threshold before optimizing branch points and calculating fitness.

IV. GLOBAL OPTIMIZATION IN THE SPACE OF HIERARCHIES
To determine the optimal hierarchy and its connectivity, we search the space of rooted, labeled, bifurcating trees.The positions of the branch points are fixed by the process in Sec.II C. The globally optimal network of all configurations maximizes the fitness F while satisfying the space-filling constraint on service volumes.As an example, the optimal configuration in Fig. 1d corresponds to the hierarchy in the bottom right, where the fitness F (1, 0) = 1 includes only total network length L [Eq. ( 1)], resulting in a Steiner tree [23]).
Our exploration of configuration space has many similarities to phylogenetic trees, for which software is available to search through the space of hierarchies [24,25].Since the available software is not tailored to our specific goals of optimizing space-filling networks, we implement our own algorithms.Because of the large number of distinct bifurcating rooted trees (that grows factorially with size), efficient search strategies generally focus on regions with greater fitness.We develop strategies to search through possible configurations and find spacefilling networks that best satisfy the general biological constraints from Sec. III.

A. Navigating in the space of hierarchies
Our numerical technique guides the search by selecting changes that increase configuration fitness.Making small changes in the branching structure, such as a single swap of two segments in the hierarchy or a regraft of one segment to a spatially near part of the tree, yields new configurations.Because the change to the hierarchy is small, using the positions of branch points from the previous configuration saves time in optimizing the global positions in the new configuration.
For local swaps in the hierarchy, we exchange a segment with one child (including the associated downstream subtrees) of the segment's sibling (i.e., its nibling).There are 2(n − 2) possible nibling swaps for n (≥ 2) discrete service volumes.However, nibling swaps do not address changes for segments that are distant in the hierarchy but have small spatial separation.To account for these changes, we regraft single segments to spatially near branches of the hierarchy.We limit the search of spatially proximal branch points to those within twice the minimum service volume separation of each other.This restriction maintains a linear increase in the number of explored regrafts with the number of service volumes, in contrast to the factorial increase that would result from including all possible regrafts.

B. Seed for trajectories: Balanced Hierarchy Construction
We accelerate the identification of near-optimal networks by choosing an initial configuration that avoids many suboptimal structures (e.g.configurations with many repeated crossings or non-contiguous subtrees).To improve overall computation time, some approaches explore permutations of pairing tips with service volumes under a constant hierarchy [26].Fortunately, the dimensionality of the space for each branch point position never exceeds three in our study, which allows us to construct a favorable configuration directly through spatial partitioning.Such a favorable configuration avoids less-fit configurations and satisfies the intuitive guidelines that branch points connect nearby subtrees (efficiency by proximity) and that sibling subtrees have similar numbers of service volumes in accordance with symmetric branching in radius.
To ensure the maximal hierarchical balance for the seed, we begin with a single set that contains all terminal service volumes.This set is then partitioned into two subsets of equal size (or within 1 service volume if the number is odd, which guarantees U ≤ 0.5), using a straight line to define the boundary between the two sets.When appropriate, this line passes through the geometric center (i.e., the unweighted average position) of the previous set of points and the geometric center of the new subset.Resulting sets further split into smaller subsets to yield a complete, bifurcating hierarchy.We refer to this process and the resulting seed as the Balanced Hierarchy Construction (BHC).

C. Efficient search trajectories
We further accelerate our search by limiting the number of nearby configurations considered at each step.We accomplish this through a carefully guided greedy search through the space of hierarchies (effectively simulated annealing [27,28] at zero temperature), which often finds a near-optimal configuration.A greedy strategy offers expedited elimination of configurations that are far from optimal; the algorithm abandons configurations that saturate at a fitness lower than the current most optimal configuration.Our implementation allows five iterations of the process in Sec.II C, then excludes configurations that fail to reduce the cost C [Eq.(3)] by at least 5% of the remaining difference from the current optimal configuration.The algorithm with this exclusion scheme successfully identifies near-optimal configurations.
Because the sampling process is not exhaustive, the search through the space of possible hierarchies is not guaranteed to yield a globally optimal configuration.However, performing a reasonably thorough search as we outline here and conducting several runs from the BHC seed (in our simulations, at least 10 runs) increase the likelihood of identifying a configuration that is near-optimal and shares many of the branching length asymmetries that an optimal configuration exhibits.With dependable algorithms for determining branch point positions and for exploring the space of possible hierarchies, we can now investigate the length properties of spacefilling networks under several basic space-filling strategies.

V. RESULTS AND ANALYSIS
We now present the results of optimized networks and of the analysis on real vascular networks, including the properties of the most optimal networks.To build intuition about the space of hierarchies, we first explore the space exhaustively for small networks and establish the distinct patterns that the two optimizations L and H produce.In comparing optimal configurations with observations of real systems, we find better agreement by enforcing a constraint on the degree of unbalance U in the hierarchies of candidate configurations.

A. Exhaustive search for small networks
To become more familiar with the landscape of possible configurations, we exhaustively explore the space of hierarchies and pairings for networks that are small enough to quickly yield comprehensive results for a single realization of fixed service volumes.We collapse and reorganize the higher-dimensional space of branch point swaps into a single dimension by ranking each configuration based on the fitness F .This reorganization involves a normalization of rank so that the fittest configuration occurs at 0 and the least fit occurs at 1. Rescaling the rank is necessary even for networks of the same size and shape because different realizations may have different numbers of unique configurations after consolidating degenerate bifurcations (mentioned in Sec.II C), despite having the same number of service volumes.Larger networks tend to have greater range with respect to both costs in L and H.
The minimum distance between each service volume is constant for each of the networks that constitute the ensemble of realizations for the curves in Fig. 2a.Each curve represents the average fitness (relative to the fittest configuration for the particular set of capillary positions) over an ensemble of networks with a fixed number of tips and constant total area.In generating the ensemble, we exclude those that arise with a different number of tips than what is desired until we accumulate 1000 configurations of the target size.Across curves, the total area increases to produce networks with more service volumes more frequently.
One might expect a large set of similarly fit, nearoptimal networks, which would be represented by a plateau near the optimum.However, the sharp descent away from the optimal configuration in Fig. 2a indicates that there are few configurations that are near-optimal.From an evolutionary perspective, this implies that the vascular networks of organisms are under strong selection.Furthermore, the slope near the optimum becomes steeper as more service volumes are introduced, so that the best configurations become more distinct from other possibilities as the number of service volumes grows.
Considering the very large number of service volumes in real organisms, this again indicates that real vascular networks are under strong selection pressures for spacefilling and efficiency.
Optimal networks that have no constraint on hierarchical balance fall into two general classifications depending on the relative weights of total network length L and average path length H in the fitness measure F .As shown in the simple examples of Fig. 3, network fitness measures that are weighted to minimize L yield bifurcating trees, while measures that are weighted to minimize H yield hubs. Bifurcating trees better correspond to real networks, suggesting that total network length L plays a larger role than average path length H. Since a single hub is not observed (and not expected from material costs) in real systems, we do not consider configurations that ignore total network length L. However, optimizing only for L leads to meandering, bifurcating paths, which become shorter and more direct by including both costs (L and H).Furthermore, additional global information is necessary to directly minimize path lengths than the local environment that we consider in Sec.II C -specifically, the context of the entire path.This means that our analysis is best suited for optimality that always includes a significant contribution from total network length L and a weaker contribution from average path length H.

B. Trajectories for sampling larger networks
With better intuition about the space of hierarchies from small networks, we now explore the space for larger networks with more service volumes.The branching properties of larger networks give more applicable results to connect particular space-filling strategies with the observations of real cardiovascular systems.We first summarize the properties of optimized networks without any constraint on hierarchical balance (U 0 = 1).
Because the search through the space of hierarchies is not exhaustive for large networks, we cannot show ranked landscapes averaged over ensembles with different service volume positions as we did for small networks.Instead, we show landscapes from a single realization of service volume positions that come from an ensemble of trajectories that start with the BHC configuration and end at a local optimum (Fig. 11a in App.B).The greedy algorithm samples fewer less-fit configurations, yielding a shallower slope near the optimum than the exhaustively explored landscapes in Fig. 2a.Since the starting point of the search (the BHC configuration) is already favorable, we expect that the worst-ranked configuration of the partial search is already very near optimal.Searches through the space of hierarchies and the properties of optimal configurations do not vary with different convex body shapes.Fig. 5 shows example optimized networks for a maximally symmetric body shape (see App. B for other shapes).The general trends of long, meandering paths for solely minimizing total network length L and of more direct paths when including H are consistent across both isotropic, circular areas and elongated, rectangular areas.
To characterize branching features of these large configurations, we quantify the asymmetric branching attributes with the two ratios choosing ℓ c1 ≤ ℓ c2 for the lengths of child 1 and 2 and r c1 ≤ r c2 for the radii (shown in Fig. 1a).Note that perfect symmetry corresponds to λ L = λ R = 1 and smaller values of λ L and λ R correspond to more asymmetric branching.Distributions for the branching asymmetry ratio in length λ L for various sizes and shapes are shown in Fig. 6.Branch points for which one of the child segments does not exist because of degeneracy with a service volume center do not contribute to  the distribution for λ L .There is little change in the features of the distribution of λ L across different sizes of networks with U 0 = 1.This trend persists for both isotropic and anisotropic enclosing shapes (as Fig. 6 in App.B shows).A summary for the cross-generational length ratio is given in App. C.
Many branch points in these networks coincide with a service volume, predicting large trunks that feed capillaries directly.Similar results appear in the study of flow through a dynamic, adaptive network [29].However, such a trend does not agree with the empirical data.Although there is asymmetry in adjacent segments at branch points and a lack of strict balance in the hierarchy along different paths, we observe that large arteries do not branch directly to capillaries and arrive at the same expectation from the dynamics of blood flow.The major qualitative distinction between the BHC and the optimized configurations with U 0 = 1 is that the BHC is a network with a balanced hierarchy.Upon inspection of empirical data in Sec.V C, we find that the branching length asymmetries for the BHC configuration (given in App.B) motivate an additional constraint on hierarchical balance during the search through the space of hierarchies.

C. Comparison of optimized networks with empirical data
The results in Sec.V B show that optimization for total network length or average path length with no constraint on hierarchical balance leads to distributions of asymmetry in sibling vessel length that skew toward symmetry (λ L ≈ 1).We now present the analysis of λ L that characterizes the local length asymmetries at branch points for real and optimized networks.From this analysis, we explore how limiting the degree of unbalance U in an optimal artificial network yields asymmetries that better match biological networks.

Asymmetric vessel length distributions of real networks
We analyze MRI images of the human head and torso as well as micro tomography images from wild-type mouse lung.Both data sets break from strict symmetry.As shown in Fig. 7, the network-wide distribution for λ R is skewed toward symmetry (λ R ≈ 1), while the distribution for λ L is more uniform, representing a greater contribution from very asymmetric branching (λ L < 1).These results are representative of general features for length distributions in real biological networks.The fact that the optimized networks in Sec.V B do not exhibit a similar distribution for λ L signals that important biological factors are missing.Because of the skew toward symmetry in sibling segment radii, we limit the hierarchical unbalance of optimized networks in Sec.V C 2.

Degree of balance necessary to match biological networks
Imposing a constraint on hierarchical balance leads to configurations that reflect more realistic asymmetry in branching lengths.Hierarchical balance, which equalizes the number of service volumes that each sibling segment supplies, is related to the blood flow that is required to deliver resources and effectively limits the asymme-

Sibling Radius Ratio (λR)
Sibling Length Ratio (λL) Mouse Lung Human Torso  6) and ( 7), respectively] in mouse lung and human torsos.(a) Radii ratios are skewed toward symmetry (λ R ≈ 1), although they are not always perfectly symmetric.(b) Length ratios are not skewed toward symmetry (many ratios have λ L < 1), contrary to symmetric models.
try of sibling radii.In Fig. 8 we show results for several thresholds for the constraint on hierarchical balance.Decreasing the threshold U 0 yields more realistic distributions for λ L , as shown in Fig. 9.Because these networks are embedded in 2-D, decreasing U 0 can also result in more crossings between segments at different levels.By comparing Figs.9a and 9b, we see that the constraint on hierarchical balance leads to similar results independent of the weight of average path length H in configuration fitness.Instead of contributing significantly to fitness, H is effectively optimized through hierarchical balance.
While enforcing hierarchical balance leads to more realistic branching and length asymmetry distributions, it is not necessary to have a maximally balanced hierarchy.In Fig. 10, we show that lowering the threshold U 0 reduces network fitness.As the constraint on hierarchical balance U 0 decreases, the average fitness also decreases.Since the distribution of λ L is approximately uniform around U 0 ≈ 0.7 and below, the best value for the hierarchical constraint is U 0 ≈ 0.7 because this yields the fittest networks that have uniform distributions for length asymmetry.

VI. DISCUSSION
With our determination of branch point positions and exploration of distinct hierarchical configurations, we can remark on several consequences that follow from the general properties of optimized networks.Organizing the lengths between branch points to fill 2-or 3-D space with capillaries inevitably leads to asymmetries and unbalanced networks [30].Strictly symmetric and balanced networks are either inefficient in materials or not spacefilling.For example, in the H-tree all children branch orthogonally from the parent, resulting in inefficient paths.Other networks with more efficient paths lead to capillaries that are equidistant from the source, which could cover the surface of a sphere but not fill its volume.For the optimal, space-filling networks that we explore, we impose a constraint that pushes the network toward hierarchically balanced branching structures but does not require maximum balance.One can imagine other interesting metrics for hierarchical balance, but we concentrate on how a maximum degree of unbalance U 0 affects the structure of the network.This guarantees a minimum level of balance in the hierarchy but still allows freedom in the search for optimal networks, as well as nonuniformity in the hierarchical balance.
We construct a seed configuration that builds a network to ensure maximal hierarchical balance while maintaining efficient contiguity of subtrees.Configurations that tend to be hierarchically balanced, such as the BHC configuration (where the constraint is implicit in the construction algorithm) or optimized configurations that limit unbalance, do not show a strong skew toward symmetric branching in lengths.This hierarchical balance may result from gradual, incremental growth as an individual organism matures and ages.Nearby vessels grow to supply resources to new tissue, resulting in contiguous subtrees and favoring routes that reduce path lengths and avoid a single, meandering artery that branches directly to capillaries.Other computational models approach the growth and optimization of space-filling networks in different ways.Although there are many algorithms to generate structure that do not intentionally optimize network architecture or space-filling properties, near-optimal configurations may emerge spontaneously from certain simple rules.Examples of such pattern formation processes and associated algorithmic rules include models for both angiogenesis [18,31], as well as vasculogenesis (in terms of chemotaxic [32,33], mechanical substratum [34], and cellular Potts models [35,36]).However, these models do not adequately address our focus on branching length asymmetries for efficient, hierarchical, space-filling networks.Specifically, the pattern formation model for angiogenesis does not incorporate consistent space-filling service volumes, only space-filling arterial structure.The arterial structure fills some regions so that they are devoid of capillaries, while multiple tips converge to the same location elsewhere.The models for vasculogenesis do not optimize the development of a hierarchical branching network.However, dynamic vascular remodeling [29] can form structures both with and without closed loops while maintaining a uniform distribution of capillaries, although the optimal structures also suffer from large arteries branching directly to capillaries.We extend these models to understand the asymmetric lengths of adjacent segments in vascular networks and how these relate to space-filling service volumes.
Because of the many different factors and interactions that influence the structure of the cardiovascular system, our basic model can be expanded in many directions.Radius information can be incorporated into optimized networks by requiring flow to be uniform in all terminal service volumes.By including radius information, blood flow as well as more appropriate structural and energetic costs can lead to revised optimization principles, which require the calculation of the weighted Fermat point (e.g., see [20]) and has been explored previously in a limited, local context [22,37,38].Note that lowering the threshold U 0 tends to increase the minimum number of branching levels between the heart and capillaries.Less drastic hierarchical unbalance implies that the ratios of parent-child radii β = r c /r p should be near 1 (symmetric branching).This translates the global, topological property into a local branching quantity.
We do not expect that increasing the dimensionality of our networks to 3-D would change the qualitative results for branching asymmetry in length (λ L ) with hierarchical balance (specifically U 0 ).However, the numerical location for an optimal trade-off between fitness and balance may shift.Studies of large vessels (near the heart) show these vessels to be planar [39]), but the planarity cannot always hold across the entire network if tips must fill a 3-D space.Still, in the absence of obstacles, all optimization conditions enforce planarity in 3-D for branch points in their local context.Introducing regions where the network is prohibited (e.g. through bones, organs or from self-avoidance) constrains the Fermat point to the surface of a sphere or some other shape [40,41].
While the topological change of allowing loops introduces many complications to the properties of flow and hierarchical labels [42], such a modification can be beneficial in understanding reticulated vascular structures.Loops are especially important when considering network robustness (i.e., resilience to damage) within organs and leaves [16,17]) or pathological growth in tumors [43,44].These types of network properties can be included in future models.
Locally, the position of a branching junction minimizes the sum of vessel lengths in our model.Globally, we impose a threshold on the minimum hierarchical balance, which reduces the differential blood flow into sibling segments.Although real vascular networks consist almost entirely of bifurcations (although there is rapid, asymmetric branching from the aorta to capillaries through coronary arteries), the iterative approach described in Sec.II C can lead to low numbers of bifurcating junctions for some candidate networks.
Limiting the degree of unbalance in the hierarchy does not continue to shift the distribution of λ L away from symmetry (λ L ≈ 1) below U 0 ≈ 0.7, which suggests that there is an appropriate trade-off between the hierarchical balance threshold U 0 and configuration fitness F that does not require perfect symmetry for an efficient network structure.The increased cost of the network in Fig. 10 is similar for both curves, implying that the increase mostly comes from total network length L.
The large number of distinct bifurcating hierarchies necessitates that we carefully choose and execute the algorithms for searching the space of possible configurations.Consequently, we construct a favorable starting point and concentrate computational resources on regions that are most likely to contain optimal configurations.Using the numerical implementations in sections II and IV, we identify optimal networks and study the length properties of individual segments within the context of a network with space-filling terminal service volumes.
Our results have many implications for how vascular networks fill space efficiently.We exhaustively explore fitness landscapes for small networks and carefully guide the sampling of the space of hierarchies for large networks in order to determine near-optimal configurations.Our results show that strict hierarchical balance is not optimal for the architecture of cardiovascular networks.Furthermore, there is a trade-off between hierarchical balance (which is related to symmetric branching in radius at the local level) and the distribution for branching in lengths that shows the connection between the space-filling and efficiency requirements of the network.
By incorporating radius and flow information, as well as growth patterns that incorporate obstacles and loops, we can continue to build on present models to better understand vascular architecture and gain insights for its effects on resource delivery, metabolic scaling, aging, and repair after damage.
Appendix A: Similarity measure to compare hierarchical groupings between configurations While collapsing the landscape of measures to a single dimension informs us about the typical distribution of configurations, it retains no information about the relation of the hierarchies between different trees.To address this issue, we define a measure of similarity to compare how two hierarchies group the same set of tips.This measure is normalized such that similar hierarchies and groupings of service volumes have a similarity score near 1, while hierarchies that group service volumes in very different ways have a similarity score near 0. To meet these guidelines, we perform a simple count of the number of identical subtree groupings between two hierarchies and normalize by the maximum possible number that could be shared if the trees were identical.In accordance with these properties, define the similarity S(A, B) between two configurations A and B as where I s is the indicator function (1 if statement s is true and 0 otherwise) and subtree refers to the set of tips in that particular subtree.
Configurations that have a worse fitness measure are less similar to the optimal configuration, as shown in Fig. 11.However, note that similarity S(T i , T opt ) is not a monotonic function when rank i is defined by the configuration's measure.For example, consider hierarchy A, which may be very similar to a hierarchy B, which itself is very similar to C. Then it is possible that A and C are less similar to each other than each is to B, yet both are ranked higher than B with respect to a particular measure.This also means that optimal configurations are not always a single swap or regraft away from all near-optimal configurations, i.e. local minima are possible.The average similarity in Fig. 11b does not approach 1, indicating that the subtree grouping of service volumes can be very different between networks that are nearly optimal (F ≈ F opt ).Some of the stratification into distinct levels of similarity is apparent in Fig. 11a for the single trajectories (the "Single BHC Seed" and "Best Random Seed").The solid line shows the average fitness over 100 trajectories from the BHC for a system with 36 service volumes.Individual paths are for a network of 48 service volumes.One path begins with the BHC configuration and the other begins with a random seed (the best performing of 5 random seeds).] bounded in circular, square, and 1 × 4 rectangular areas are all similarly skewed toward symmetry.However, the distributions in BHC networks are all skewed away from symmetry.
[14] Yunlong Huo and Ghassan S.  [44] V. M. Savage, A. B. Herman, G. B. West, and K. Leu, "Using fractal geometry and universal growth curves as diagnostics for comparing tumor vasculature and metabolic rate with healthy tis-sue and for predicting responses to drug therapies," Discrete and Continuous Dyn.Syst.B 18, 1077-1108 (2013).

Figure 1 :
Figure 1: (a) Schematic of the simplest bifurcating tree network, showing the heart (hollow red triangle) and two service volumes (filled brown circles) with labels for the lengths of each associated segment.(b) Three possible locations for a bifurcation junction.The rightmost configuration shows the Fermat point of △123 that minimizes the sum of segment lengths.(c) The two distinct bifurcations (left) collapse to a single trifurcation (center) and set to the geometric median of the four endpoints (right).(d) Comparing C (a measure of some length property of the network) for each of the 15 configurations for four service volumes shows that the bottom right configuration is optimal with respect to C. The relative order and position of both the tips and branch points in Fig. 1d do not correspond to relative positions in space.

Fitness Landscape F ( 1 Figure 2 :
Figure 2: (Color) Ranked landscapes for unique consolidated configurations from exhaustive exploration of bifurcating trees of fixed area and fixed N tips (a) Exhaustive landscapes for fitness based only on total network length L [F (1, 0)].(b) Exhaustive landscapes for fitness based total network length L and average path length H [F (1, 9)].We choose the weights (C L , C H ) = (1, 9) so that the contribution from H is not dominated by the contribution from L.

Figure 3 :
Figure 3: (Color) Two classes of networks: (a) The optimal configuration that minimizes total network length [L in Eq. (1)] of the 15 possible trees (corresponding to Fig. 1d) consists only of bifurcations.(b) The optimal configuration that minimizes average path length between each service volume and the heart [H from Eq. (2)] of the 15 possible trees consists of a single hub.The regions of varying background color define the Voronoi cells corresponding to individual service volumes.

Figure 4 :
Figure 4: (Color) Average fitness landscapes for total network length L [F (1, 0)] over 100 trajectories for a single network of each size.
Figure 5: (Color) Optimal configurations for two fitness measures.We choose the weights (C L , C H ) = (1, 9) so that the contribution from H is not dominated by the contribution from L.

Figure 8 :
Figure 8: (Color) Optimal configurations for several constraints on hierarchical balance and two optimization weights for fitness F .The BHC seed is the same as in Fig. 5.We choose the weights (C L , C H ) = (1, 9) so that the contribution from H is not dominated by the contribution from L.

Figure 9 :Figure 10 :
Figure 9: (Color) Distributions for several thresholds of hierarchical unbalance U 0 .(a) Distributions of λ L optimized solely for total network length L [F (1, 0)].(b) Distributions for λ L optimized for total network length L and average path lengths H [F (1, 9)].All plots are averaged over 200 realizations of service volume distributions.

Figure 11 :
Figure 11: (Color) Partial search trajectories in the the space of hierarchies.(a) Ranked fitness landscapes for the partially explored hierarchical space during optimization for total network length L [F (1, 0)].(b) Corresponding similarity measures.All quantities are calculated with respect to the most optimal configuration found over all runs.The solid line shows the average fitness over 100 trajectories from the BHC for a system with 36 service volumes.Individual paths are for a network of 48 service volumes.One path begins with the BHC configuration and the other begins with a random seed (the best performing of 5 random seeds).

Figure 12 :
Figure 12: (Color) Network properties with increased size.(a) The number of unique configurations N conf ig for each fixed number of service volumes is narrowly distributed relative to the increased number of configurations from introducing an additional service volume.(b) The average number of service volumes increases directly proportional to the total area.The solid grey line represents the standard deviation about the mean.

Figure 14 :
Figure 14: (Color) The distributions for λ L in optimized networks [F (1, 0) and F (1, 9)] bounded in circular, square, and 1 × 4 rectangular areas are all similarly skewed toward symmetry.However, the distributions in BHC networks are all skewed away from symmetry.

Figure 15
Figure 15: (Color) Distributions for the ratio of child-to-parent length γ.(a) Distributions of γ in mouse lung and human head and torso.(b) Distributions of γ for several thresholds of hierarchical unbalance U 0 , optimized solely for total network length L [F (1, 0)].