Asymmetries arising from the space-filling nature of vascular networks
aa r X i v : . [ phy s i c s . m e d - ph ] A ug Asymmetries arising from the space-filling nature of vascular networks
David Hunt and Van M. Savage
Department of Biomathematics, University of California at Los Angeles, Los Angeles, CA 90095, USA (Dated: September 25, 2018)Cardiovascular networks span the body by branching across many generations of vessels. Theresulting structure delivers blood over long distances to supply all cells with oxygen via the relativelyshort-range process of diffusion at the capillary level. The structural features of the network thataccomplish this density and ubiquity of capillaries are often called space-filling . There are multiplestrategies to fill a space, but some strategies do not lead to biologically adaptive structures byrequiring too much construction material or space, delivering resources too slowly, or using toomuch power to move blood through the system. We empirically measure the structure of realnetworks (18 humans and 1 mouse) and compare these observations with predictions of modelnetworks that are space-filling and constrained by a few guiding biological principles. We devisea numerical method that enables the investigation of space-filling strategies and determination ofwhich biological principles influence network structure. Optimization for only a single principlecreates unrealistic networks that represent an extreme limit of the possible structures that could beobserved in nature. We first study these extreme limits for two competing principles, minimal totalmaterial and minimal path lengths. We combine these two principles and enforce various thresholdsfor balance in the network hierarchy, which provides a novel approach that highlights the trade-offsfaced by biological networks and yields predictions that better match our empirical data.
PACS numbers: 87.10.Vg 89.75.Da 89.75.Fb 89.75.Hc 89.75.Kd
I. INTRODUCTION
The vital functions of the cardiovascular system are thedistribution of oxygen and nutrient resources throughoutthe body, as well as the collection and filtration of wasteby circulating blood. Transfer of resources and waste oc-curs primarily at the capillary level via diffusion throughnearby tissue. Consequently, this smallest level of thenetwork must reach all living cells in order to maintainthem, filling the entire space of the body. In modelsdeveloped by Krogh for effective diffusion of oxygen [1],cells cannot survive beyond a maximum distance from acapillary. This defines a service volume of cells that is as-sociated with each capillary, which has a typical size thathas been observed to vary across species based on cellularmetabolic rate [2, 3]. The constraint on maximum dis-tance from capillaries necessitates that the final levels ofthe cardiovascular network are also space-filling through-out the body. In this paper we investigate the relationbetween this space-filling property and basic optimiza-tion principles such as the minimization of costs fromconstruction material and pumping power. Specifically,we highlight how this relation influences the asymmetriesin sizes and flow rates of sibling segments as measured inempirical data.A central focus of our investigation of cardiovascularsystems is the space-filling properties of networks, butthese properties are also of great interest in many othercontexts. General space-covering hexagonal patterns ap-pear in nature in the cell structure of beehives as wellas in economic theories for market areas [4]. Trees (thewoody, perennial plants) have been studied for both howforests fill an area [5], as well as how the vascular struc-ture within an individual plant determines the hydraulics of resource delivery [6, 7]. Apollonian networks [8] de-scribe the space-filling packing of spheres of various sizes,which are similar in the cardiovascular system to consid-ering the volumes of different subtrees of the network.Efficiently filling space in two dimensions is importantfor information visualization [9]. In addition to these ap-plications, Kuffner and LaValle study space-filling treenetworks (i.e., networks that branch and have no closedloops) to determine a route from one location to another[10], but without the biological constraints that we im-pose here. For cardiovascular networks, this motion plan-ning is analogous to how the vascular structure routesblood. Efficient routing of blood through a hierarchy iscentral to models that investigate allometric scaling ofmetabolic rate with body mass [2, 3, 11–14], which buildon metabolic scaling by Kleiber [15]. Determining spe-cific space-filling strategies will inform these models tobetter describe the cardiovascular system.Developmental processes (i.e., growth) as well as evo-lutionary pressures, such as efficiency in material andenergy use, shape the structure of cardiovascular net-works. Filling a volume or surface efficiently with theterminal nodes of a branching network is nontrivial, es-pecially when the distribution system must reliably de-liver blood at each stage of development. The systemmust also be robust to changes in vessel lengths and thenumber of hierarchical levels that result from growth andobstructions from damage [16, 17]. We propound thatthese structural challenges lead to the asymmetric fea-tures in both the local relations of segments, as well asin the global paths from the heart to each service volumethat we observe in empirical data.At the local level, the ratio of lengths between par-ent and child segments may vary across the network, de-viating from strict self-similarity. Additionally, siblingsegments may vary in length, which we call asymmetricbranching . Within our numerical method, these featuresarise as a result of optimizing branch point positions rel-ative to adjacent branch points across the network. Vari-ation in the length and number of levels in paths meansthat the network is also not symmetric or balanced inthese global properties. By allowing these asymmetriesand explicitly ensuring space-filling structure, we expandother models that are strictly balanced in network hier-archy and perfectly symmetrical in local quantities.Asymmetries observed in real systems motivate ourinvestigation of the space-filling properties and asym-metries in cardiovascular networks. These observationsshow that such networks have a tendency for the lengthsof sibling segments to be distributed less symmetricallythan is the case for radii. The empirical data in Sec. Vcome from two types of sources. Images collected throughmicrotomography of mouse lung comprise one data set.The mice were part of a study with different manipula-tions of matrix GLA protein (which causes the vascula-ture to be under- or over-branched [18]), but we focuson the wild-type specimen for our analysis. The otherdata set, collected through MRI, excludes the lungs andinstead focuses on the central vasculature in the humanhead and torso [19]. We utilize the custom software Angi-cart [19] to analyze these two distinct vascular data sets.Because of the noninvasive nature of these data acqui-sition and analysis techniques, future studies have theopportunity to track the development of cardiovascularsystems as individual organisms grow and age, includingrepair after the system incurs damage (i.e., wound heal-ing). Such studies can elucidate the effects that patternsof growth and changes from damage have on the final,mature state of the network.In this paper, we study the optimization principles thatcorrespond to evolutionary pressures for efficiency in ma-terial cost and power loss. Our focus is the influence ofspace-filling patterns on length asymmetry distributionswithout the explicit inclusion of radius information. Thelist of candidate networks includes all possible hierar-chical (topological) connections between the heart andall capillaries. For each hierarchy and unique permu-tation of pairings between terminal vessels and servicevolumes (see Fig. 1d), we must determine the positionsof the branch points. The combination of the hierarchy,service-volume pairings and branch point positions de-fines the configuration of a candidate network. For thesereasons, we must search through many candidate config-urations to determine the most efficient structure. Wequantify the fitness of each candidate network using in-dividual segment lengths between branch points as wellas full path lengths between each capillary and the heart.To perform a reliable comparison between candidateconfigurations, it is crucial to determine branch pointpositions in a consistent way. We determine these posi-tions iteratively for the entire network in order to identifythe global optimum. While the local process of choosingbranch points that minimize total vessel lengths (or sim- ilar features) is relatively straightforward to iterate overthe network, any single branch point and its relation toits neighbors relies indirectly on updates that are appliedelsewhere. This dependence emerges from the fact thateach end of a vessel is connected to a branch point, whichupon repositioning affects the lengths of all vessels that itjoins. The uniqueness of the
Fermat point — the branchpoint position that minimizes the sum of the lengths ofvessels at a single junction — is already well established(for example, see [20]). This allows us to carefully con-struct an algorithm (described in Sec. II C) that reliablyrelaxes all branch point positions into the global opti-mum.After determining the positions of branch points fora given hierarchy, we compare distinct configurations tofind the optimum network. The search through configu-rations is also a central problem in phylogenetics, wherethe goal is to construct phylogenies to identify similargroups of species and trace the development of genesthrough speciation. Even in the case of genes that controlbiological traits, a loosely analogous space-filling phe-nomenon emerges in the form of species filling the nichesin the environment. With our specific goal of completespatial covering of network tips, we develop strategies inSec. IV for exploring the space of hierarchies that aresimilar to those used on phylogenetic trees.The organization of the subsequent sections is as fol-lows. In Sec. II, we describe the basic assumptions forour space-filling network model, including the details ofthe local optimization of branch point locations and theglobal paths through the network. In Sec. III, we in-troduce the specific quantitative network properties thatwe use to compare the fitness of candidate networks. Weintroduce the properties of the space of tree hierarchiesand our implemented exploration strategies in Sec. IV.In Sec. V, we detail the results from the several layers ofoptimization that we implement, and discuss the insightsthat they offer in Sec. VI.
II. CONSTRUCTION OF ARTIFICIALVASCULAR NETWORKS
To better understand the connection between the localasymmetries of individual vessel lengths and the globalconstraint on space-filling capillaries, we optimize candi-date networks that are embedded in two spatial dimen-sions (2- D ) in silico with respect to specific optimizationprinciples. We explore these optimized artificial networksand quantify their branching length asymmetries to com-pare with our empirical data. Our model’s simplificationof the cardiovascular network focuses on the lengths ofsegments as defined by the straight line between adja-cent branch points. Reticulated structures occur withinleaves to mitigate damage [16, 17] and within animals asanastomoses (or pathologically as fistulas). However, wefocus on the vast majority of the cardiovascular systemthat distributes resources through a hierarchical, tree-likestructure, in which no segments or subtrees rejoin down-stream to form closed loops before reaching the capillar-ies. This is sufficient for the focus of our investigation ofthe asymmetric, space-filling structure that distributesthe resource-rich blood from the heart throughout thebody.The space-filling property of the cardiovascular net-work constrains the hierarchical structure of the networkand the positions of branch points. Here we describeour process for the construction of individual networksand the space of possible networks through the follow-ing steps: defining a distribution of space-filling servicevolumes in the space of the body, identifying all uniquehierarchies and pairings between tips of hierarchies anddistinct service volumes, and determining the positionsof branch points for each hierarchy and pairing. A. Space-filling service volumes
Because real systems do not organize or grow on a reg-ular (symmetric, isotropic) grid, we position service vol-umes randomly within the space they fill. Constructionof service volumes begins by choosing a random pointwithin the body volume that represents the location of acapillary. We then randomly choose other points (capil-lary locations) so that none lie within a predefined con-stant distance from another capillary location. After de-termining a set of capillary locations that span the 2-Darea, the entire body is partitioned into Voronoi cells fedby the closest capillary. In this way, each capillary be-comes associated with a specific service volume, and thesum of the service volumes fills the entire space (see Fig.3 or 5).
B. Space of hierarchically distinct trees andpairings with service volumes
Because multiple branching levels connect the servicevolumes to the heart, there are many possible hierarchi-cal orderings of branching junctions across these differ-ent levels. For example, there are two unique hierarchieswhen there are four service volumes: the top three con-figurations in Fig. 1d are the same perfectly balancedhierarchy, while the remaining trees have the same unbal-anced hierarchy. The distinguishing feature is the pair-ings of tips in the hierarchy to specific service volumes(1-4).There are many distinct pairings of terminal tips withthe set of fixed service volumes for each bifurcating tree.For four service volumes, there are a manageable 15unique bifurcating trees (shown in Fig. 1d). Before de-termining branch point positions for small networks, weexhaustively search through all possible hierarchies andpairings with service volumes. For large networks, thenumber of distinct hierarchies and pairings ((2 n − n distinct service volumes) is prohibitively large, so wesample the space as described in Sec. V B.We do not disqualify configurations if one vessel pathcrosses with another (these would likely separate in threedimensions), and there is no exchange of resources or in-teraction in blood flow at such locations. Crossings arenot observed for networks that minimize only total net-work length without a constraint on hierarchical balance,but they often occur for optimal configurations with astrong constraint on hierarchical balance. C. Optimization of branch point positions for afixed hierarchy and pairing
We now detail our algorithm for the optimization ofthe positions of branch points that connect a distribu-tion of service volumes to the heart through a hierarchi-cal network. Within our algorithm, the position of eachbranch point depends solely on the location of the adja-cent branch points in the network. Distant vessels affecteach other indirectly, but not through any direct long-range process. Using the limited, local information givenby the neighborhood of a branch point, each junction isassigned a uniquely-defined position that minimizes thesum of the Euclidean distances to each neighboring junc-tion. This is equivalent to the Fermat point of the tri-angle formed by the two downstream ends of the childvessels and the one upstream end of the parent vessel(see Fig. 1b). The Fermat point of a triangle is a specialcase of the more general geometric median , the uniquepoint that minimizes the sum of distances to an arbi-trary number of other fixed points. We follow the algo-rithm presented in [21] to avoid errors in determining thegeometric median. Assigning branch point positions asgeometric medians effectively minimizes the constructioncosts for the local network structure.We construct our networks from simple bifurcations,but using the Fermat point to assign branch point po-sitions can lead to coincident (degenerate) bifurcations,as shown in Fig. 1c. Degenerate branch points are con-solidated at the geometric median of the upstream end-point of the parent and three or more downstream end-points of the associated children segments. In this way,two degenerate bifurcations become a trifurcation and,more generally, n degenerate bifurcations become a sin-gle ( n + 1)-furcation. Networks that are hierarchicallydistinct in their bifurcating structure can become iden-tical networks by collapsing bifurcations. Through ex-haustive explorations (described in Sec. V A), we findthat this marginally reduces the number of possible con-figurations (see results in Fig. 12a). Because we have no a priori filter to identify which bifurcating trees are re-dundant, we must consider the entirety of the space of la-beled, rooted, bifurcating trees throughout the algorithmto identify sufficiently optimal configurations. With po-sitions defined in a consistent way, we can now comparethe properties of distinct hierarchies to determine which Global Optimization via Search throughthe Space of Hierarchies (a)
Local Optimization of BranchPoint Position (b)(c) (d)
Figure 1: (a) Schematic of the simplest bifurcating tree network, showing the heart (hollow red triangle) and twoservice volumes (filled brown circles) with labels for the lengths of each associated segment. (b) Three possiblelocations for a bifurcation junction. The rightmost configuration shows the Fermat point of △
123 that minimizesthe sum of segment lengths. (c) The two distinct bifurcations (left) collapse to a single trifurcation (center) and setto the geometric median of the four endpoints (right). (d) Comparing C (a measure of some length property of thenetwork) for each of the 15 configurations for four service volumes shows that the bottom right configuration isoptimal with respect to C . The relative order and position of both the tips and branch points in Fig. 1d do notcorrespond to relative positions in space.is the best for a particular space-filling strategy. III. SELECTION CRITERIA FOR BIOLOGICALNETWORKS
All characteristics of an organism that affect fitnessand are heritable are under selection. A key question iswhich features of the vascular network are under selec-tion. Here we define specific fitness measures that aretied to the structure of the network configuration thatallows us to rank candidate networks and determine theoptimal configuration.
A. Global length properties of space-fillingconfigurations
Here we introduce a standardized fitness measure thatallows us to compare candidate networks for their suit-ability to transport blood and resources. Each indepen- dent measure for a network’s fitness relates to a physicalquantity that likely guides the evolution of the cardio-vascular system toward a more efficient network. Specif-ically, the system’s cost in construction material and themaintained volume of the blood relates to the total net-work length — the sum of the lengths of all segments.Competing with this minimization of materials, the dis-sipation of the heart’s pumping power relates to the pathlengths between each capillary and the heart. The powerdissipated by smooth (Poiseuille) flow through a segmentis directly proportional to the length of the segment [22].In the absence of radius information, reducing the cost ofpumping blood is equivalent to reducing the total pathlengths that blood travels.We define these two fitness measures — one dealingwith total network length and the other with individualpath lengths — as L = X all segments iin network ℓ i (1) H = X all paths pin network X all segments iin path p ℓ i (2)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 combinesthese two measures by their respective weights C L and C H : C ( C L , C H ) ≡ C L L + C H H. (3)This cost function connects minimization of material andpower dissipation to study optimal networks that arespace-filling.Because an increase in cost corresponds to a decrease infitness, we place this approach in an evolutionary frame-work by inverting and scaling our cost function C to bea relative fitness function F : F ( C L , C H ) ≡ C min C ( C L , C H ) (4)where C min is the most optimal network under consider-ation. By this definition, the optimal configuration hasa fitness F = 1 and less optimal configurations have afitness F < B. Equal Distribution of Resources throughHierarchical Balance
Because the network tends to exhibit nearly symmetricbranching in radius and must distribute resources equallyto each capillary in the body, the network hierarchy can-not be overly unbalanced, with one segment having manymore tips to supply than its sibling. In accordance withthis argument, empirical data do not show major arter-ies branching directly into capillaries. We address thisconstraint by selecting for networks with more balancedhierarchies.A hierarchy is better balanced if there are roughlyequal numbers of tips supplied downstream by each sib-ling segment. Conversely, a hierarchy becomes morepoorly balanced as the disparity grows between the num-ber of tips. In this sense, we define the degree that ahierarchy is unbalanced U as U = 1 − min all siblingpairs (i , j) ( N ( i ) tips N ( j ) tips ) (5)where N ( i ) tips is the number of distinct downstream servicevolumes supplied through segment i and segment j isalways the sibling with the most downstream tips. Inour algorithm, we select against hierarchies for which thedegree of unbalance U is greater than some threshold U .We eliminate configurations above this threshold beforeoptimizing branch points and calculating fitness. IV. GLOBAL OPTIMIZATION IN THE SPACEOF HIERARCHIES
To determine the optimal hierarchy and its connectiv-ity, we search the space of rooted, labeled, bifurcatingtrees. The positions of the branch points are fixed by theprocess in Sec. II C. The globally optimal network of allconfigurations maximizes the fitness F while satisfyingthe space-filling constraint on service volumes. As an ex-ample, the optimal configuration in Fig. 1d correspondsto 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 simi-larities to phylogenetic trees, for which software is avail-able to search through the space of hierarchies [24, 25].Since the available software is not tailored to our spe-cific goals of optimizing space-filling networks, we imple-ment our own algorithms. Because of the large numberof distinct bifurcating rooted trees (that grows factori-ally with size), efficient search strategies generally focuson regions with greater fitness. We develop strategiesto search through possible configurations and find space-filling networks that best satisfy the general biologicalconstraints from Sec. III. A. Navigating in the space of hierarchies
Our numerical technique guides the search by select-ing changes that increase configuration fitness. Makingsmall changes in the branching structure, such as a singleswap of two segments in the hierarchy or a regraft of onesegment to a spatially near part of the tree, yields newconfigurations. Because the change to the hierarchy issmall, using the positions of branch points from the pre-vious configuration saves time in optimizing the globalpositions in the new configuration.For local swaps in the hierarchy, we exchange a seg-ment with one child (including the associated down-stream subtrees) of the segment’s sibling (i.e., its nib-ling). There are 2( n −
2) possible nibling swaps for n ( ≥
2) discrete service volumes. However, niblingswaps do not address changes for segments that are dis-tant in the hierarchy but have small spatial separation.To account for these changes, we regraft single segmentsto spatially near branches of the hierarchy. We limitthe search of spatially proximal branch points to thosewithin twice the minimum service volume separation ofeach other. This restriction maintains a linear increasein the number of explored regrafts with the number ofservice volumes, in contrast to the factorial increase thatwould result from including all possible regrafts.
B. Seed for trajectories: Balanced HierarchyConstruction
We accelerate the identification of near-optimal net-works by choosing an initial configuration that avoidsmany suboptimal structures (e.g. configurations withmany repeated crossings or non-contiguous subtrees). Toimprove overall computation time, some approaches ex-plore permutations of pairing tips with service volumesunder a constant hierarchy [26]. Fortunately, the di-mensionality of the space for each branch point posi-tion never exceeds three in our study, which allows us toconstruct a favorable configuration directly through spa-tial partitioning. Such a favorable configuration avoidsless-fit configurations and satisfies the intuitive guidelinesthat branch points connect nearby subtrees (efficiency byproximity) and that sibling subtrees have similar num-bers of service volumes in accordance with symmetricbranching in radius.To ensure the maximal hierarchical balance for theseed, we begin with a single set that contains all ter-minal service volumes. This set is then partitioned intotwo subsets of equal size (or within 1 service volume ifthe number is odd, which guarantees U ≤ . C. Efficient search trajectories
We further accelerate our search by limiting the num-ber of nearby configurations considered at each step. Weaccomplish this through a carefully guided greedy searchthrough the space of hierarchies (effectively simulated an-nealing [27, 28] at zero temperature), which often findsa near-optimal configuration. A greedy strategy offersexpedited elimination of configurations that are far fromoptimal; the algorithm abandons configurations that sat-urate at a fitness lower than the current most optimalconfiguration. Our implementation allows five iterationsof the process in Sec. II C, then excludes configurationsthat fail to reduce the cost C [Eq.(3)] by at least 5% ofthe remaining difference from the current optimal config-uration. The algorithm with this exclusion scheme suc-cessfully identifies near-optimal configurations.Because the sampling process is not exhaustive, thesearch through the space of possible hierarchies is notguaranteed to yield a globally optimal configuration.However, performing a reasonably thorough search aswe outline here and conducting several runs from theBHC seed (in our simulations, at least 10 runs) increasethe likelihood of identifying a configuration that is near- optimal and shares many of the branching length asym-metries that an optimal configuration exhibits. With de-pendable algorithms for determining branch point posi-tions and for exploring the space of possible hierarchies,we can now investigate the length properties of space-filling networks under several basic space-filling strate-gies. V. RESULTS AND ANALYSIS
We now present the results of optimized networks andof the analysis on real vascular networks, including theproperties of the most optimal networks. To build in-tuition about the space of hierarchies, we first explorethe space exhaustively for small networks and establishthe distinct patterns that the two optimizations L and H produce. In comparing optimal configurations with ob-servations of real systems, we find better agreement byenforcing a constraint on the degree of unbalance U inthe hierarchies of candidate configurations. A. Exhaustive search for small networks
To become more familiar with the landscape of pos-sible configurations, we exhaustively explore the spaceof hierarchies and pairings for networks that are smallenough to quickly yield comprehensive results for a sin-gle realization of fixed service volumes. We collapse andreorganize the higher-dimensional space of branch pointswaps into a single dimension by ranking each configura-tion based on the fitness F . This reorganization involvesa normalization of rank so that the fittest configurationoccurs at 0 and the least fit occurs at 1. Rescaling therank is necessary even for networks of the same size andshape because different realizations may have differentnumbers of unique configurations after consolidating de-generate bifurcations (mentioned in Sec. II C), despitehaving the same number of service volumes. Larger net-works tend to have greater range with respect to bothcosts in L and H .The minimum distance between each service volumeis constant for each of the networks that constitute theensemble of realizations for the curves in Fig. 2a. Eachcurve represents the average fitness (relative to the fittestconfiguration for the particular set of capillary positions)over an ensemble of networks with a fixed number of tipsand constant total area. In generating the ensemble, weexclude those that arise with a different number of tipsthan what is desired until we accumulate 1000 configu-rations of the target size. Across curves, the total areaincreases to produce networks with more service volumesmore frequently.One might expect a large set of similarly fit, near-optimal networks, which would be represented by aplateau near the optimum. However, the sharp descentaway from the optimal configuration in Fig. 2a indicates Fitness Landscape F (1 ,
0) Fitness Landscape F (1 , M ea s u r e / ( B e s t M ea s u r e ) N tips = 4 N tips = 5 N tips = 6 N tips = 7 (a) M ea s u r e / ( B e s t M ea s u r e ) N tips = 4 N tips = 5 N tips = 6 N tips = 7 (b) Figure 2: (Color) Ranked landscapes for unique consolidated configurations from exhaustive exploration ofbifurcating trees of fixed area and fixed N tips (a) Exhaustive landscapes for fitness based only on total networklength L [ F (1 , L and average path length H [ F (1 , C L , C H ) = (1 ,
9) so that the contribution from H is not dominated by thecontribution from L .that there are few configurations that are near-optimal.From an evolutionary perspective, this implies that thevascular networks of organisms are under strong selec-tion. Furthermore, the slope near the optimum becomessteeper as more service volumes are introduced, so thatthe best configurations become more distinct from otherpossibilities as the number of service volumes grows.Considering the very large number of service volumesin real organisms, this again indicates that real vascularnetworks are under strong selection pressures for space-filling and efficiency.Optimal networks that have no constraint on hierarchi-cal balance fall into two general classifications dependingon the relative weights of total network length L and av-erage path length H in the fitness measure F . As shownin the simple examples of Fig. 3, network fitness mea-sures that are weighted to minimize L yield bifurcatingtrees, while measures that are weighted to minimize H yield hubs. Bifurcating trees better correspond to realnetworks, suggesting that total network length L playsa larger role than average path length H . Since a sin-gle hub is not observed (and not expected from materialcosts) in real systems, we do not consider configurationsthat ignore total network length L . However, optimizingonly for L leads to meandering, bifurcating paths, whichbecome shorter and more direct by including both costs( L and H ). Furthermore, additional global information isnecessary 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 anal-ysis is best suited for optimality that always includes asignificant contribution from total network length L anda weaker contribution from average path length H . B. Trajectories for sampling larger networks
With better intuition about the space of hierarchiesfrom small networks, we now explore the space for largernetworks with more service volumes. The branchingproperties of larger networks give more applicable resultsto connect particular space-filling strategies with the ob-servations of real cardiovascular systems. We first sum-marize the properties of optimized networks without anyconstraint on hierarchical balance ( U = 1).Because the search through the space of hierarchies isnot exhaustive for large networks, we cannot show rankedlandscapes averaged over ensembles with different servicevolume positions as we did for small networks. Instead,we show landscapes from a single realization of servicevolume positions that come from an ensemble of trajec-tories that start with the BHC configuration and end ata local optimum (Fig. 11a in App. B). The greedy al-gorithm samples fewer less-fit configurations, yielding ashallower slope near the optimum than the exhaustivelyexplored landscapes in Fig. 2a. Since the starting point Global Optimization forTotal Network Length L [ F (1 , H [ F (0 , (a) (b) Figure 3: (Color) Two classes of networks: (a) Theoptimal configuration that minimizes total networklength [ L in Eq. (1)] of the 15 possible trees(corresponding to Fig. 1d) consists only of bifurcations.(b) The optimal configuration that minimizes averagepath length between each service volume and the heart[ H from Eq. (2)] of the 15 possible trees consists of asingle hub. The regions of varying background colordefine the Voronoi cells corresponding to individualservice volumes.of the search (the BHC configuration) is already favor-able, we expect that the worst-ranked configuration ofthe partial search is already very near optimal.Searches through the space of hierarchies and the prop-erties of optimal configurations do not vary with differ-ent convex body shapes. Fig. 5 shows example opti-mized networks for a maximally symmetric body shape(see App. B for other shapes). The general trends oflong, meandering paths for solely minimizing total net-work length L and of more direct paths when including H are consistent across both isotropic, circular areas andelongated, rectangular areas.To characterize branching features of these large con-figurations, we quantify the asymmetric branching at-tributes with the two ratios λ L = ℓ c ℓ c (6) λ R = r c r c (7)choosing ℓ c ≤ ℓ c for the lengths of child 1 and 2 and r c ≤ r c for the radii (shown in Fig. 1a). Note thatperfect symmetry corresponds to λ L = λ R = 1 andsmaller values of λ L and λ R correspond to more asym-metric branching. Distributions for the branching asym-metry ratio in length λ L for various sizes and shapesare shown in Fig. 6. Branch points for which one ofthe child segments does not exist because of degener-acy with a service volume center do not contribute to F it n e ss F ( , ) N tips = 60 N tips = 52 N tips = 39 N tips = 34 N tips = 22 Figure 4: (Color) Average fitness landscapes for totalnetwork length L [ F (1 , λ L . There is little change in thefeatures of the distribution of λ L across different sizesof networks with U = 1. This trend persists for bothisotropic and anisotropic enclosing shapes (as Fig. 6 inApp. B shows). A summary for the cross-generationallength ratio is given in App. C.Many branch points in these networks coincide witha service volume, predicting large trunks that feed cap-illaries directly. Similar results appear in the study offlow through a dynamic, adaptive network [29]. How-ever, such a trend does not agree with the empirical data.Although there is asymmetry in adjacent segments atbranch points and a lack of strict balance in the hierar-chy along different paths, we observe that large arteriesdo not branch directly to capillaries and arrive at thesame expectation from the dynamics of blood flow. Themajor qualitative distinction between the BHC and theoptimized configurations with U = 1 is that the BHCis a network with a balanced hierarchy. Upon inspectionof empirical data in Sec. V C, we find that the branch-ing length asymmetries for the BHC configuration (givenin App. B) motivate an additional constraint on hier-archical balance during the search through the space ofhierarchies. BHC Seed Total Network Length L [ F (1 , L and MeanPath Length H [ F (1 , Figure 5: (Color) Optimal configurations for two fitness measures. We choose the weights ( C L , C H ) = (1 ,
9) sothat the contribution from H is not dominated by the contribution from L . λ L P ( λ L ) 〈 N tips 〉 = 47 〈 N tips 〉 = 89 〈 N tips 〉 = 175 〈 N tips 〉 = 336 Figure 6: (Color) Distributions of λ L for several sizes ofcircular areas, averaged over 100 realizations of servicevolume distributions and optimized solely for totalnetwork length L [ F (1 , C. Comparison of optimized networks withempirical data
The results in Sec. V B show that optimization fortotal network length or average path length with no con-straint on hierarchical balance leads to distributions of asymmetry in sibling vessel length that skew toward sym-metry ( λ L ≈ λ L that characterizes the local length asymmetries at branchpoints for real and optimized networks. From this anal-ysis, we explore how limiting the degree of unbalance U in an optimal artificial network yields asymmetries thatbetter match biological networks.
1. Asymmetric vessel length distributions of real networks
We analyze MRI images of the human head and torsoas well as micro tomography images from wild-typemouse lung. Both data sets break from strict symmetry.As shown in Fig. 7, the network-wide distribution for λ R is skewed toward symmetry ( λ R ≈ λ L is more uniform, representing a greatercontribution from very asymmetric branching ( λ L < λ L signals that important bio-logical factors are missing. Because of the skew towardsymmetry in sibling segment radii, we limit the hierar-chical unbalance of optimized networks in Sec. V C 2.
2. Degree of balance necessary to match biological networks
Imposing a constraint on hierarchical balance leads toconfigurations that reflect more realistic asymmetry inbranching lengths. Hierarchical balance, which equalizesthe number of service volumes that each sibling segmentsupplies, is related to the blood flow that is requiredto deliver resources and effectively limits the asymme-0
Sibling Radius Ratio ( λ R ) Sibling Length Ratio ( λ L ) λ R P ( λ R ) Mouse LungHuman Torso increasing asymmetry perfect symmetry (a) λ L P ( λ L ) Mouse LungHuman Torso increasing asymmetry perfect symmetry (b)
Figure 7: (Color) Observed radius and length branching asymmetry ratios [Eqs. (6) and (7), respectively] in mouselung and human torsos. (a) Radii ratios are skewed toward symmetry ( λ R ≈ λ L < U yields more realistic distri-butions for λ L , as shown in Fig. 9. Because these net-works are embedded in 2- D , decreasing U can also resultin more crossings between segments at different levels.By comparing Figs. 9a and 9b, we see that the constrainton hierarchical balance leads to similar results indepen-dent of the weight of average path length H in config-uration fitness. Instead of contributing significantly tofitness, H is effectively optimized through hierarchicalbalance.While enforcing hierarchical balance leads to more re-alistic branching and length asymmetry distributions, itis not necessary to have a maximally balanced hierarchy.In Fig.10, we show that lowering the threshold U reducesnetwork fitness. As the constraint on hierarchical balance U decreases, the average fitness also decreases. Sincethe distribution of λ L is approximately uniform around U ≈ . U ≈ . VI. DISCUSSION
With our determination of branch point positions andexploration of distinct hierarchical configurations, we canremark on several consequences that follow from the gen-eral properties of optimized networks. Organizing thelengths between branch points to fill 2- or 3- D spacewith capillaries inevitably leads to asymmetries and un-balanced networks [30]. Strictly symmetric and balancednetworks are either inefficient in materials or not space-filling. For example, in the H-tree all children branch or-thogonally from the parent, resulting in inefficient paths.Other networks with more efficient paths lead to capil-laries that are equidistant from the source, which couldcover the surface of a sphere but not fill its volume. Forthe optimal, space-filling networks that we explore, weimpose a constraint that pushes the network toward hi-erarchically balanced branching structures but does notrequire maximum balance. One can imagine other inter-esting metrics for hierarchical balance, but we concen-trate on how a maximum degree of unbalance U affectsthe structure of the network. This guarantees a min-imum level of balance in the hierarchy but still allowsfreedom in the search for optimal networks, as well asnonuniformity in the hierarchical balance.We construct a seed configuration that builds a net-work to ensure maximal hierarchical balance while main-taining efficient contiguity of subtrees. Configurations1 Total Network Length L [ F (1 , L and Mean Path Length H [ F (1 , U = 1 . U = 0 . U = 0 . Figure 8: (Color) Optimal configurations for severalconstraints on hierarchical balance and twooptimization weights for fitness F . The BHC seed is thesame as in Fig. 5. We choose the weights( C L , C H ) = (1 ,
9) so that the contribution from H isnot dominated by the contribution from L .that tend to be hierarchically balanced, such as the BHCconfiguration (where the constraint is implicit in theconstruction algorithm) or optimized configurations thatlimit unbalance, do not show a strong skew toward sym-metric branching in lengths. This hierarchical balancemay result from gradual, incremental growth as an indi-vidual organism matures and ages. Nearby vessels growto supply resources to new tissue, resulting in contiguoussubtrees and favoring routes that reduce path lengths andavoid a single, meandering artery that branches directlyto capillaries.Other computational models approach the growth andoptimization of space-filling networks in different ways.Although there are many algorithms to generate struc-ture that do not intentionally optimize network archi-tecture or space-filling properties, near-optimal configu-rations may emerge spontaneously from certain simple rules. Examples of such pattern formation processes andassociated algorithmic rules include models for both an-giogenesis [18, 31], as well as vasculogenesis (in termsof chemotaxic [32, 33], mechanical substratum [34], andcellular Potts models [35, 36]). However, these modelsdo not adequately address our focus on branching lengthasymmetries for efficient, hierarchical, space-filling net-works. Specifically, the pattern formation model for an-giogenesis does not incorporate consistent space-fillingservice volumes, only space-filling arterial structure. Thearterial structure fills some regions so that they are de-void of capillaries, while multiple tips converge to thesame location elsewhere. The models for vasculogenesisdo not optimize the development of a hierarchical branch-ing network. However, dynamic vascular remodeling [29]can form structures both with and without closed loopswhile maintaining a uniform distribution of capillaries,although the optimal structures also suffer from large ar-teries branching directly to capillaries. We extend thesemodels to understand the asymmetric lengths of adja-cent segments in vascular networks and how these relateto space-filling service volumes.Because of the many different factors and interactionsthat influence the structure of the cardiovascular system,our basic model can be expanded in many directions.Radius information can be incorporated into optimizednetworks by requiring flow to be uniform in all terminalservice volumes. By including radius information, bloodflow as well as more appropriate structural and energeticcosts can lead to revised optimization principles, whichrequire the calculation of the weighted Fermat point (e.g.,see [20]) and has been explored previously in a limited, lo-cal context [22, 37, 38]. Note that lowering the threshold U tends to increase the minimum number of branchinglevels between the heart and capillaries. Less drastic hier-archical unbalance implies that the ratios of parent-childradii β = r c /r p should be near 1 (symmetric branching).This translates the global, topological property into alocal branching quantity.We do not expect that increasing the dimensionality ofour networks to 3- D would change the qualitative resultsfor branching asymmetry in length ( λ L ) with hierarchicalbalance (specifically U ). However, the numerical loca-tion for an optimal trade-off between fitness and balancemay shift. Studies of large vessels (near the heart) showthese vessels to be planar [39]), but the planarity cannotalways hold across the entire network if tips must fill a3- D space. Still, in the absence of obstacles, all opti-mization conditions enforce planarity in 3- D for branchpoints in their local context. Introducing regions wherethe network is prohibited (e.g. through bones, organs orfrom self-avoidance) constrains the Fermat point to thesurface of a sphere or some other shape [40, 41].While the topological change of allowing loops intro-duces many complications to the properties of flow andhierarchical labels [42], such a modification can be ben-eficial in understanding reticulated vascular structures.Loops are especially important when considering network2 λ L for Total NetworkLength L [ F (1 , λ L for Total NetworkLength L and Average PathLength H [ F (1 , λ L P ( λ L ) U = 0.00 U = 0.05 U = 0.10 U = 0.20 U = 0.30 U = 0.40 U = 0.50 (a) λ L P ( λ L ) U = 0.00 U = 0.05 U = 0.10 U = 0.20 U = 0.30 U = 0.40 U = 0.50 (b) Figure 9: (Color) Distributions for several thresholds of hierarchical unbalance U . (a) Distributions of λ L optimizedsolely for total network length L [ F (1 , λ L optimized for total network length L and averagepath lengths H [ F (1 , U 〈 F 〉 F (1, 9) F (1, 0) Figure 10: (Color) Average fitnesses F for ensembles of100 optimized configurations. Configurations are morefit if a greater hierarchical unbalance U is allowed.robustness (i.e., resilience to damage) within organs andleaves [16, 17]) or pathological growth in tumors [43, 44].These types of network properties can be included in fu-ture models. Locally, the position of a branching junction minimizesthe sum of vessel lengths in our model. Globally, we im-pose a threshold on the minimum hierarchical balance,which reduces the differential blood flow into sibling seg-ments. Although real vascular networks consist almostentirely of bifurcations (although there is rapid, asym-metric branching from the aorta to capillaries throughcoronary arteries), the iterative approach described inSec. II C can lead to low numbers of bifurcating junc-tions for some candidate networks.Limiting the degree of unbalance in the hierarchy doesnot continue to shift the distribution of λ L away fromsymmetry ( λ L ≈
1) below U ≈ .
7, which suggeststhat there is an appropriate trade-off between the hierar-chical balance threshold U and configuration fitness F that does not require perfect symmetry for an efficientnetwork structure. The increased cost of the networkin Fig. 10 is similar for both curves, implying that theincrease mostly comes from total network length L .The large number of distinct bifurcating hierarchies ne-cessitates that we carefully choose and execute the algo-rithms for searching the space of possible configurations.Consequently, we construct a favorable starting point andconcentrate computational resources on regions that aremost likely to contain optimal configurations. Using thenumerical implementations in sections II and IV, we iden-tify optimal networks and study the length properties ofindividual segments within the context of a network withspace-filling terminal service volumes.Our results have many implications for how vascular3networks fill space efficiently. We exhaustively explorefitness landscapes for small networks and carefully guidethe sampling of the space of hierarchies for large net-works in order to determine near-optimal configurations.Our results show that strict hierarchical balance is notoptimal for the architecture of cardiovascular networks.Furthermore, there is a trade-off between hierarchicalbalance (which is related to symmetric branching in ra-dius at the local level) and the distribution for branch-ing in lengths that shows the connection between thespace-filling and efficiency requirements of the network.By incorporating radius and flow information, as well asgrowth patterns that incorporate obstacles and loops, wecan continue to build on present models to better under-stand vascular architecture and gain insights for its ef-fects on resource delivery, metabolic scaling, aging, andrepair after damage. Appendix A: Similarity measure to comparehierarchical groupings between configurations
While collapsing the landscape of measures to a singledimension informs us about the typical distribution ofconfigurations, it retains no information about the rela-tion of the hierarchies between different trees. To addressthis issue, we define a measure of similarity to comparehow two hierarchies group the same set of tips. Thismeasure is normalized such that similar hierarchies andgroupings of service volumes have a similarity score near1, while hierarchies that group service volumes in verydifferent ways have a similarity score near 0. To meetthese guidelines, we perform a simple count of the num-ber of identical subtree groupings between two hierar-chies and normalize by the maximum possible numberthat could be shared if the trees were identical. In accor-dance with these properties, define the similarity S ( A, B )between two configurations A and B as S ( A, B ) ≡ σ ( A, B )max { σ ( A, A ) , σ ( B, B ) } σ ( X, Y ) = X subtree min network X X subtree nin network Y ( I m ⊆ n + I n ⊆ m )where I s is the indicator function (1 if statement s is trueand 0 otherwise) and subtree refers to the set of tips inthat particular subtree.Configurations that have a worse fitness measure areless similar to the optimal configuration, as shown in Fig.11. However, note that similarity S ( T i , T opt ) is not amonotonic function when rank i is defined by the con-figuration’s measure. For example, consider hierarchy A,which may be very similar to a hierarchy B, which it-self is very similar to C. Then it is possible that A and C are less similar to each other than each is to B, yetboth are ranked higher than B with respect to a par-ticular measure. This also means that optimal config-urations are not always a single swap or regraft awayfrom all near-optimal configurations, i.e. local minimaare possible. The average similarity in Fig. 11b does notapproach 1, indicating that the subtree grouping of ser-vice volumes can be very different between networks thatare nearly optimal ( F ≈ F opt ). Some of the stratificationinto distinct levels of similarity is apparent in Fig. 11afor the single trajectories (the “Single BHC Seed” and“Best Random Seed”). Appendix B: Network size and shape
In Fig. 12a we show the distribution of the number ofdistinct hierarchies after consolidating degenerate bifur-cations for these ensembles of fixed numbers of servicevolumes. The variance in the number of unique config-urations increases with network size, but a dominatingcontribution that increases the number of configurationscomes from adding a service volume. Not surprisingly,the average total number of service volumes h N tips i scaleslinearly with the total area (see Fig. 12b). In Fig. 14we show distributions for branching length asymmetry λ L . All distributions for optimized networks that haveno constraint on the hierarchy (i.e. U = 1) exhibit astrong skew toward symmetry. Appendix C: Parent-child length ratio γ In Fig. 15a we show the network-wide distribution forthe ratio of child-to-parent lengths γ = ℓ c /ℓ p for child c with parent p . Although there is the tendency that γ <
1, some child segments have a relatively shorter parent.Although slight, an increased threshold for U shifts morechild segments to be shorter than their associated parent(see Fig. 15b). Independent of the threshold, the nonzerovariance of this distribution shows that γ is not constantthroughout the network. ACKNOWLEDGMENTS
We would like to thank Kristina I. Bostrm, MD,PhD and Yucheng Yao, MD, PhD for sharing theirdata of mouse lung vasculature. We would also like tothank Daniel Ennis, PhD for sharing his data of humanhead and torso vasculature and Mitchell Johnson for hiswork in developing the software to analyze vascular im-ages. We are grateful to Eric Deeds, PhD and TomKolokotrones, MD, MPH, PhD for engaging in stimu-lating discussions about the work. We would also liketo thank Elif Tekin for her patient help in refining thepresentation in this article.4 F it n e ss F ( , ) Single BHC SeedBest Random SeedAverage BHC Seed (a) F , Partial Search)0.50.60.70.80.91 S i m il a r it y S ( T , T op t ) Single BHC SeedAverage BHC SeedBest Random Seed (b)
Figure 11: (Color) Partial search trajectories in the the space of hierarchies. (a) Ranked fitness landscapes for thepartially explored hierarchical space during optimization for total network length L [ F (1 , N config -4 -2 P ( N c on fi g ) N tips = 4 N tips = 5 N tips = 6 N tips = 7 (a) 〈 N ti p s 〉 (b) Figure 12: (Color) Network properties with increased size. (a) The number of unique configurations N config for eachfixed number of service volumes is narrowly distributed relative to the increased number of configurations fromintroducing an additional service volume. (b) The average number of service volumes increases directly proportionalto the total area.The solid grey line represents the standard deviation about the mean. [1] August Krogh, “The supply of oxygen to the tis-sues and the regulation of the capillary circulation,”J. Physiol. , 457–474 (1919). [2] G. B. West, J. H. Brown, and B. J. Enquist, “A generalmodel for the origin of allometric scaling laws in biology,”Science , 122–126 (1997). BHC Seed Total Network Length L [ F (1 , L and Mean Path Length H [ F (1 , Figure 13: (Color) Optimal configurations for several shapes and two fitness measures. Both isotropic (circular) andanisotropic (rectangular) areas exhibit similar space-filling strategies for near-optimal configurations. We choose theweights ( C L , C H ) = (1 ,
9) so that the contribution from H is not dominated by the contribution from L . [3] V. M. Savage, E. J. Deeds, and W. Fontana,“Sizing up allometric scaling theory,”PLoS Comput. Biol. , e1000171 (2008).[4] T¨onu Puu, “On the genesis of hexagonal shapes,”Networks and Spatial Economics , 5–20 (2005).[5] Geoffrey B West, Brian J Enquist, andJames H Brown, “A general quantitativetheory of forest structure and dynamics,”Proc. Nat. Acad. Sci. U.S.A. , 7040–7045 (2009).[6] VM Savage, LP Bentley, BJ Enquist, JS Sperry,DD Smith, PB Reich, and EI von Allmen, “Hy-draulic trade-offs and space filling enable better pre-dictions of vascular structure and function in plants,”Proc. Nat. Acad. Sci. U.S.A. , 22722–22727 (2010).[7] J. S. Sperry, D. D. Smith, V. M. Savage, B. J. En-quist, K. A. McCulloh, P. B. Reich, L. P. Bent-ley, and E. I. von Allmen, “A species-level modelfor metabolic scaling in trees i. exploring bound-aries to scaling space within and across species,”Functional Ecology , 1054–1065 (2012).[8] Jos´e S Andrade Jr, Hans J Herrmann, Roberto FSAndrade, and Luciano R Da Silva, “Apollonian net- works: Simultaneously scale-free, small world, eu-clidean, space filling, and with matching graphs,”Phys. Rev. Lett. , 018702 (2005).[9] John Stasko, Richard Catrambone, Mark Guzdial, andKevin McDonald, “An evaluation of space-filling informa-tion visualizations for depicting hierarchical structures,”Int. J. Human-Computer Stud. , 663–694 (2000).[10] James J Kuffner and Steven M LaValle,“Space-filling trees: A new perspective on in-cremental search for motion planning,” in Intelligent Robots and Systems (IROS), 2011 IEEE/RSJ International Conference on (IEEE, 2011) pp. 2199–2206.[11] Jayanth R Banavar, Amos Maritan, and Andrea Ri-naldo, “Size and form in efficient transportation net-works,” Nature , 130–132 (1999).[12] Jayanth R Banavar, John Damuth, AmosMaritan, and Andrea Rinaldo, “Supply–demand balance and metabolic scaling,”Proc. Nat. Acad. Sci. U.S.A. , 10506–10509 (2002).[13] Peter Sheridan Dodds, “Optimal form ofbranching supply and collection networks,”Phys. Rev. Lett. , 048702 (2010). | | | | | | | | | || | | | | | | | | || | | | | | | | | |0 0.2 0.4 0.6 0.8 1 λ L P ( λ L ) BHC F (1, 0) F (1, 9)CircleSquare1 × Figure 14: (Color) The distributions for λ L in optimizednetworks [ F (1 ,
0) and F (1 , × toward symmetry. However, the distributions inBHC networks are all skewed away from symmetry. [14] Yunlong Huo and Ghassan S. Kassab, “In-traspecific scaling laws of vascular trees,”J. R. Soc. Interface , 190–200 (2012).[15] M. Kleiber, “Body size and metabolism,”Hilgardia , 315 (1932).[16] Francis Corson, “Fluctuations and redun-dancy in optimal transport networks,”Phys. Rev. Lett. , 048703 (2010).[17] E. Katifori, G. J. Sz¨oll˝osi, and M. O. Magnasco, “Dam-age and fluctuations induce loops in optimal transportnetworks,” Phys. Rev. Lett. , 048704 (2010).[18] Y. Yao, S. Nowak, A. Yochelis, A. Garfinkel,and K. I. Bost¨om, “Matrix GLA protein, an in-hibitory morphogen in pulmonary vascular develop-ment,” J. Biol. Chem. , 30131–42 (2007).[19] M. G. Johnson, D. Ennis, and V. M. Savage, “Quantify-ing statistical self similarity in humanvessel networks us-ing data extracted automatically from MRI,” PLoS Com-put. Biol. (to be published) (2015).[20] Y. Shen and J. Tolosa, “The weighted fermat triangleproblem,” Int. J. Math. Math. Sci. (2008).[21] Germ´an A Torres, “A weiszfeld-like algorithm for a weberlocation problem constrained to a closed and convex set,”arXiv preprint arXiv:1204.1087 (2012).[22] M. Zamir, The Physics of Coronary Bloof Flow (Springer Science+Business Media, Inc., 2005).[23] D. Sankoff and P. Rousseau, “Locating the ver-tices of a steiner tree in an arbitrary metric space,”Math. Programming , 240–246 (1975).[24] Fredrik Ronquist, Maxim Teslenko, Paul van der Mark,Daniel L Ayres, Aaron Darling, Sebastian H¨ohna, BretLarget, Liang Liu, Marc A Suchard, and John P Huelsen-beck, “Mrbayes 3.2: efficient bayesian phylogenetic in-ference and model choice across a large model space,”Systematic Biol. , 539–542 (2012).[25] St´ephane Guindon, Jean-Fran¸cois Dufayard, Vin-cent Lefort, Maria Anisimova, Wim Hordijk, and Olivier Gascuel, “New algorithms and meth-ods to estimate maximum-likelihood phyloge-nies: assessing the performance of phyml 3.0,”Systematic Biol. , 307–321 (2010).[26] Griffin Weber, Lucila Ohno-Machado, andStuart Shieber, “Representation in stochas-tic search for phylogenetic tree reconstruction,”J. Biomed. Info. , 43–50 (2006).[27] A. Dress and M. Kr¨uger, “Parsimonious phyloge-netic trees in metric spaces and simulated annealing,”Adv. Appl. Math. , 8–37 (1987).[28] L. A. Salter and D. K. Pearl, “Stochastic search strat-egy for estimation of maximum likelihood phylogenetictrees,” Systematic Biol. , 7–17 (2001).[29] D. Hu and D. Cai, “Adaptation and opti-mization of biological transport networks,”Phys. Rev. Lett. , 138701 (2013).[30] Jayanth R Banavar, Melanie E Moses, James HBrown, John Damuth, Andrea Rinaldo, Richard MSibly, and Amos Maritan, “A general ba-sis for quarter-power scaling in animals,”Proc. Nat. Acad. Sci. U.S.A. , 15816–15820 (2010).[31] H. Meinhardt, “Morphogenesis of lines and nets,”Differentiation , 117–123 (1976).[32] G. Serini, D. Ambrosi, E. Giraudo, A. Gamba,L. Preziosi, and F. Bussolino, “Modeling theearly stages of vascular network assembly,”EMBO J. , 1771–1779 (2003).[33] A. Gamba, D. Ambrosi, A. Coniglio, A. de Can-dia, S. Di Talia, E. Giraudo, G. Serini, L. Preziosi,and F. Bussolino, “Percolation, morphogenesis,and burgers dynamics in blood vessels formation,”Phys. Rev. Lett. , 118101 (2003).[34] D. Manoussaki, S.R. Lubkin, R.B. Vemon,and J.D. Murray, “A mechanical model forthe formation of vascular networks in vitro,”Acta Biotheoretica , 271–282 (1996).[35] R. M. H. Merks and J. A. Glazier, “Dynamic mechanismsof blood vessel growth,” Nonlinearity , C1 (2006).[36] R. M. H. Merks and P. Koolwijk, “Modelingmorphogenesis in silico and in vitro: Towardsquantitative, predictive, cell-based modeling,”Math. Modelling of Natural Phenomena , 149–171 (2009).[37] M. Zamir, “Nonsymmetrical bifurcations in arterialbranching,” J. Gen. Physiol. , 837–845 (1978).[38] M. Zamir, The Physics of Pulsatile Flow (Springer NewYork, 2005).[39] T. Wischgoll, J. S. Choy, and G. S. Kassab, “Ex-traction of morphometry and branching anglesof porcine coronary arterial tree from ct images,”Am. J. Physiol. Heart Circ. Physiol. , H1949–1955 (2009).[40] A. N. Zachos, “An analytical solution of the weightedfermat-torricelli problem on the unit sphere,”arXiv preprint arXiv:1408.6495 (2014).[41] A. N. Zachos, “Exact location of the weightedfermat–torricelli point on flat surfaces of revolution,”Results Math. , 167–179 (2014).[42] Y. Mileyko, H. Edelsbrunner, C. A. Price, and J. S.Weitz, “Hierarchical ordering of reticular networks,”PLoS ONE , e36715 (2012).[43] A. B. Herman, V. M. Savage, and G. B. West, “A quan-titative theory of solid tumor growth, metabolic rate andvascularization,” PLoS One , e22973 (2011). γ P ( γ ) HumanMouse (a) γ P ( γ ) U = 0.00 U = 0.05 U = 0.10 U = 0.20 U = 0.30 U = 0.40 U = 0.50 (b) Figure 15: (Color) Distributions for the ratio of child-to-parent length γ . (a) Distributions of γ in mouse lung andhuman head and torso. (b) Distributions of γ for several thresholds of hierarchical unbalance U , optimized solelyfor total network length L [ F (1 , [44] V. M. Savage, A. B. Herman, G. B. West, andK. Leu, “Using fractal geometry and universalgrowth curves as diagnostics for comparing tumorvasculature and metabolic rate with healthy tis- sue and for predicting responses to drug therapies,”Discrete and Continuous Dyn. Syst. B18