General framework for branching angle optimization and asymmetry

Systematic patterns in branching angles or asymmetric ratios [28, 45–47] (see Results, Analysis of asymmetry patterns in vascular data) suggest that a constraint that likely arises through natural selection, the nature of the growth process, or both. Most hypotheses about the force of natural selection on vasculature have focused on principles that reduce the cost of materials and growth while also providing efficient flow mechanisms [27, 34–37]. These biological principles could apply at each branching junction (locally) during growth or across the whole network (globally) through some larger bauplan. In this study, we initially focus on the local optimization aspects of these principles. Then, we consider applying different spatial (regional) constraints on the branching structure, including simulating networks that have random branching within constraints that range from intermediate to global spatial levels of the network.

For the local optimization of bifurcating branching geometry, following a similar approach as in the previous studies by Zamir and Murray [33, 34, 37], we assume that the radius (r_i) and the unshared endpoints of the vessels, i.e., the vertices (V_i), are fixed, whereas the branching junction at the shared endpoint (J) varies (Fig 1C and 1D). In this framework, we derive the optimal placement of the branching junction (J) by constructing and minimizing a cost function based on each biological principle that is hypothesized to increase the fitness. In all our derivations, we follow Krogh’s model that regards blood vessels as cylinders [48, 49]. In addition, both human and mouse data provide evidence that all of the vessels involved in a branching junction lie within a single plane (S1 Fig). Thus, for our derivations of optimal branching geometry, we assume that the branching junction lies in the plane determined by the unshared endpoints of all vessels.

To try to elucidate the cause of the high degree of length asymmetry observed in data (see Results), we also take a new approach and derive the length asymmetry from the solutions of these optimization principles. That is, by finding optimal branching angle solutions, we determine the location of the branching junction that in turn uniquely determines vessel lengths and leads to predictions for optimal length asymmetry. The formula that associates the length asymmetry with the optimal branching angles is introduced in the Appendix A.

As mentioned above, we consider two general principles. First, because construction and maintenance can be expensive to the body, minimizing total material is a potential driving factor for the structure of the vascular system. Such a principle will be referred to as material-cost optimization (MC optimization). The total material cost across the bifurcation is the linear sum of the material cost for each vessel because material cost is an additive quantity over different vessels.

Next, a viable design for the vascular system requires efficient flow mechanisms—such as minimal power loss—to transport nutrients and oxygen to cells. This scheme will be called power-cost optimization (PC optimization). In order to establish a cost function that represents the total power loss across a bifurcation, we need to use appropriate fluid-dynamical concepts. Because the effective impedance is not a linearly additive quantity for combining vessels at a branching junction, PC optimization requires a more complicated cost function than MC optimization. In this respect, our derivation differs from and corrects previous branching angle optimizations by Zamir [33, 34]. As another new element to our approach, we further propose an optimization scheme that can incorporate the power cost due to the downstream impedances beyond just a single branching junction.

In the following sections, we first introduce the branching angle optimization solutions for the MC and PC optimizations for a local branching junction. Next, we explain a new scheme for the PC optimization that incorporates information about the flow properties of the downstream vessels. Lastly, we relax the optimization principles and present simulations of random branching to explore how constraints—local to intermediate to global—can alone affect the characteristics of the vascular network.

Material-cost (MC) optimizations

There are two types of material that are needed for the vascular system, blood vessels (endothelial cells) and blood (plasma and white and red blood cells). The amount of material necessary for vessel construction primarily depends on the surface-area of the vessel 2πrl. In contrast, the material devoted to the blood is proportional to the blood volume (πr²l). MC optimizations can be built upon these two different characteristics (surface-area (MC-SA) or volume (MC-V)). We consider both cases here.

To consider total surface-area and volume as distinct structural constraints that need to be minimized to conserve construction material, as in Murray and Zamir [33, 34, 37], we form a generic cost function H in which lengths are weighted according to the corresponding optimization. Explicitly, H = ∑_i h_il_i where h_i = 2πr_i for surface-area and h_i = πr_i² for volume. In general, the constant geometric factors like 2π can be ignored because they cancel from every term when the derivative of H is set equal to zero. As introduced above, we assume that the radius and the unshared endpoints for each vessel are fixed when performing the optimization calculation. For fixed radii (or equivalently for fixed cost terms h_i) and the vertices, the MC optimization problem is equivalent to the weighted Fermat Problem introduced by Greenberg and Robertello [50–52].

To solve these minimization problems, we present a different method than that used to obtain Zamir and Murray’s solutions. Our method relies on distance metrics without invoking a coordinate system (Appendix B) and is therefore simpler and more general. The method is straightforward when one realizes that the location of the branching junction is uniquely defined by the parent vessel length (l₀) and the angle (φ₁) between the parent vessel and the edge determined by the unshared endpoints of itself (V₀) and one of its daughter vessels (V₁) (see Fig 1D). Based on this, we obtain the full solution of optimal branching angles by finding the stationary and singular points of H with respect to l₀ and φ₁ throughout the entire space.

In this framework, we recognize that the first order derivatives of H are discontinuous and thus undefined at the unshared endpoints (V₀, V₁, V₂). Hence, singularities (values of infinity) are attained at these endpoints. Moreover, we find that the stationary solution for minimizing H exists only in the interior of the triangle defined by the unshared endpoints. In the latter case, the optimal branching angle solution is (1)

Note that our branching angles are defined differently than in Zamir [33, 34], so these equations are equivalent but not identical. Recognizing our definitions are relative to the parent rather than centerline extended from the parent, these expressions can be translated into Zamir’s solutions by subtracting our angles θ₁ and θ₂ from π. Moreover, in our 2D (planar) optimization scheme, we use the conventional counter-clockwise direction to define the angles so that the trigonometric functions have consistent signs (Fig 1D). Correctly defining the directionality of the angles is needed for our mathematical derivation, but for comparing to empirical data, we compute only the magnitude of the branching angles. Although it is more challenging to define directionality of angles for the real, empirical data in three dimensions, it is still straightforward to measure the magnitude of the angle between the straight lines defined by the endpoints of vessels (see Materials and Methods), allowing us to directly compare with the solutions from the 2D optimization.

Notice that the branching angle solution for the stationary case does not exist when the right side of the above expressions are not in the interval [-1,1]. For instance, when h₀ = 2.1 and h₁ = h₂ = 1, then cos θ₀ = 1.21 and θ₀ is not defined. Values outside of the allowed interval occur frequently when substituting values from real vasculature, so these equations should not be blindly applied without considering the allowed regions and intervals. Moreover, there are cases for which an optimal branching solution can be computed from the above equations, but the resulting branching junction does not lie inside the triangle defined by the fixed vessel endpoints. In this case, the computed optimum does not correspond to the true optimum, which actually occurs at one of the fixed vessel endpoints (i.e., the vertices of the triangle). Indeed, for both of these cases where the stationary solution either cannot be computed or the computed answer lies outside the pre-defined triangle, the optimal solution is actually attained at one of the singularity points, i.e., the vertices V_i.

Based on these observations, we find the conditions on the cost parameters and the geometry of the fixed endpoints that correspond to degenerate branching solutions. For example, when the cost of one vessel exceeds the total cost by the other vessels (h_i ≥ h_j + h_k), then the optimal branching occurs at the unshared endpoint of the vessel i (V_i) to eradicate this particularly costly vessel. On the other hand, when the solution θ_i defined by (Eq 1) is less than the angle V_jV_iV_k (i.e., cos θ_i > cos V_jV_iV_k), which we refer to as the triangle condition, then the branching junction collapses onto the vertex V_i. The details of the proof are given in the Appendix C.

These degeneracies are not specified in previous work [33, 34, 37] because vessels are assumed to have volume flow rates that are proportional to r³. With this presumption and by the conservation of volume flow rate across each branching junction, those studies assume a strict relationship between radii of the parent and daughter vessels (Generalized Murray’s law: r₀^d = r₁^d + r₂^d, where d ∈ [2,3] [36]) that avoids the degeneracy (S2 Text). Moreover, previous work does not explicitly consider cases in which the branching angle solution cannot satisfy the triangle conditions defined above. Importantly, we find that these ignored cases and conditions correspond to the vast majority of values calculated from empirical branching vessels (S2 Fig).

Power-cost optimization for a single branching junction (PC-0)

Another biological principle and property to optimize is the power loss for pumping blood from the heart to the capillaries. This principle is tantamount to minimizing the total power for circulating blood or the total power lost that represents additional power beyond what is used to move the blood. In this way, as much of the additional power as possible is then redistributed and devoted to other needs such as foraging and reproduction [2, 5]. Much of the power that is required to push blood through the parent and daughter vessels is lost due to dissipation, especially in the smaller vessels that dominate the numbers and energetics of the whole vascular network. We calculate the power dissipated by a single vessel in terms of the volume flow rate of the blood () and the impedance (Z) via the formula [27, 35, 53]. To correctly compute the total power loss of all vessels connected by a single branching junction—a parent and daughter vessels, we must employ rules of fluid dynamics to combine impedances into an equivalent impedance. This equivalent impedance must in turn be used in the cost function for the power minimization (see below), and as we show, the correct version is different than the simple linear addition rule used by Zamir for both structural (correctly) and flow (incorrectly) constraints [33, 34].

For any collection of vessels, total volume flow rate is defined by flow through a single cross section that cuts through all vessels at the same branching level. Equivalent impedance Z_eq for the collection of vessels is defined by , based on fluid dynamics (analogue to Ohm’s Law), where Δp is the pressure difference across the extreme endpoints through which blood flows in and out of the entire collections of vessels. The total power loss for the collected vessels is then calculated from , where all of this is based on standard rules of fluid dynamics.

For a single branching junction, the optimization of power loss is derived for the case that there is a source with a constant rate of flow, , entering the parent vessel. Based on our power loss equation (), the PC optimization is then equivalent to minimizing the equivalent impedance Z_eq of the branching junction. By following the direction of the flow, we calculate the equivalent impedance by noting that vessels at the same level (i.e., sibling vessels) are in a parallel configuration, whereas the vessels across levels (i.e., parent and daughter vessels) are in a series configuration. Pressure drop across sibling vessels may be asymmetric but to simplify the calculation of an equivalent impedance, we further posit that the pressure drop across each daughter vessel is approximately the same, i.e., Δp₁ = Δp₂ ≔ Δp_d. Representing total volume flow rate of the daughter vessels by and equivalent impedance of the daughter vessels by Z_d,eq, we have Δp_d/Z_d,eq = Δp₁/Z₁ + Δp₂/Z₂. Canceling the pressure terms gives . By invoking conservation of fluid, we can combine this expression with the parent vessel in series to obtain for the equivalent impedance of all the vessels that connect at a single branching junction.

Assuming smooth, Poiseuille flow, the impedance of a single vessel is given by , where μ is the viscosity of the blood [24, 27, 53, 54]. Thus, for each vessel we implement a cost function of the form hl, where h is the cost per length, i.e., . Putting all of this together to find a solution for the PC-0 optimization (power loss at a single branching junction), our goal is to find the position of the branching junction that minimizes

By numerically calculating Z_eq as a function of the junction point and by using high resolution in space, we generate heat maps that illustrate the behavior of Z_eq. These heat maps reveal that the branching junction always collapses onto one of the vertices—unshared vessel endpoints—for any values of vessel radii (S3 Fig). Based on our numerical evidence, we show analytically that power loss and equivalent impedance, Z_eq, attain minima only at a vertex of our original triangle that is defined by the unshared endpoints of the vessel. The specific vertex at which the junction collapses is determined both by the vessel radii and the relative locations of the unshared endpoints of the vessels. The proof follows from two steps. First, we calculate the equivalent impedance at each vertex to see how specific cases of cost parameters (h₀,h₁,h₂) determine which vertex is the best location for the branching junction J. Second, we show that when the junction is located within the triangle, Z_eq is always greater than the minimum value of the impedances when the junction is at a specific vertex. Together, this proves that the minimum of Z_eq is attained at one of the vertices (Appendix D). Therefore, the PC-0 optimization leads to a degenerate branching geometry by completely eliminating the vessel that is most costly.

Thus, minimizing power loss at a single branching junction leads to either no branching at all throughout the entire cardiovascular network or to a single hub at the heart with a long, individual vessel to each terminal tip—the most downstream vessel of the network. This architecture would result in extremely large numbers of vessels directly connected to the heart because terminal tips correspond to capillaries (numbering in the millions or billions) for complete networks. This is unrealistic as sequential branching is the most noticeable and perhaps most important feature of the vascular system or any resource-distribution network. Therefore, we consider how to modify the PC minimization to attain more biologically realistic results that include vascular branching and thus lead to predictions for the branching angles and vessel lengths. In the next section, we describe how to adapt the power-minimization calculation to include larger sections of the network that expand beyond a single branching junction to encompass downstream vessels and branching junctions.

Power-cost optimization beyond a single branching junction (PC-1)

The solutions above thus lead to either a single path or to the heart being a single hub that separately connects to each capillary through a direct and independent vessel. Either configuration for vascular architecture would prevent the vascular system from efficiently distributing blood to all the downstream vessels, capillaries, and cells. To overcome these problems, we include larger sections of the vascular network by incorporating impedances of downstream vessels. To do this generically, we recognize that downstream vessels are in series with each of the daughter vessels, so we can represent the downstream impedance by adding constant terms c₁ and c₂ to the impedance of the daughter vessels. Thus, the equivalent impedance of the bifurcation becomes

For the case of vessels above the capillaries, these constant terms represent the impedance from all of the downstream vessels. For the capillary case, these constant terms are still not zero, however, because they represent the minimum impedance of a capillary, which is not allowed to be zero. We now investigate different geometries for a single branching junction for which the downstream impedances c₁ and c₂ are constant.

We simplify this problem using two general principles. First, when the impedance of the parent vessel and the branching daughter vessels are matched, no pulsatile reflections occur and the power loss at the bifurcation is minimized [27, 53]. Second, we assume the simple case that the siblings have identical impedances and each sibling has the same number of downstream vessels. From these two assumptions, we find that the equivalent downstream impedance is much larger than single vessel impedances except for vessels close to the capillaries (Appendix E).

Because the number and small surface-areas of capillaries will again likely dominate the power loss for the network, we solve the optimization problem when the downstream impedances are large, i.e., c_i ≫ Z_i. Consequently, we expand as a series to first order and obtain the following approximation:

We note that the constant term can be ignored because its derivative is zero and hence does not influence the optimization. As a result, we want to find the location of the branching junction that minimizes

The coefficients and are always less than 1, so the cost per length for the daughter vessels are diminished by these rescaling constants, thus reducing the likelihood the solution will collapse onto a daughter vessel endpoint. Defining the non-dimensionalized ratio k ≔ c₁/c₂, the optimization function becomes

This further implies that the cost of the daughter vessel with the larger downstream impedance is diminished less than the cost of the other daughter vessel. This pushes the optimal branching junction towards the daughter vessel with smaller downstream impedance.

The new optimization function for power loss above (i.e., PC-1) has the same form as the cost function for the material-cost (MC) optimizations, except that the costs per length for daughter vessels are rescaled by terms that depend on k. We therefore simplify the notation and define . With these definitions for , we can immediately use our results for Eq (1) to obtain the stationary solution for which the branching junction occurs inside the triangle of the unshared endpoints.

As for Eq 1 and the MC optimizations, the stationary solution does not exist or does not provide the minimum and the degenerate solution occurs at vertex V_i for the following cases: 1) or 2) . These two conditions correspond to six inequalities in terms of k( = c₁/c₂) if we take all combination of i, j, and k. Solving these, we obtain the values of k that result in solutions within the triangle versus those that collapse on a vertex. Mapping these values into the c₁c₂-plane yields predictable lines that separate the plane into regions classified as non-degenerate (unction within the triangle defined by the vertices) or degenerate (collapse onto a vertex) solutions. Indeed, by using approximations to solve the above limiting case (), we predict the full solution space. The validity of our approximate solutions for this PC-1 optimization problem is further explored in Results (Analytical solutions for power-cost optimization beyond single branching (PC-1)) by comparing them with numerical solutions.

Expanding from local to global constraints for the random placement of branching junctions

The above optimization schemes are based on the local consideration of branching junctions: at each branching the unshared endpoints (V_i) are fixed and the branching junction (J) is attained within the triangle of these endpoints (Fig 1D). To explore the effects of the size of the constraint region and also to observe the effects of random branching, possibly resulting from developmental processes, at different spatial (regional) constraints, we now consider relaxing the locality assumption to various degrees from fully local to fully global.

We simulate the branching network by randomly placing branching junctions within spatial regions that range from local boundaries (i.e., the positions of adjacent branching junctions) to global boundaries (i.e., within some maximum distance of the center of the network but otherwise unrestricted). For all the simulations, the hierarchical ordering of vessels, the location of the terminal tips, and most upstream branching nodes (i.e., the source) are the same as in the empirical data. We created 100 realizations for each simulation type to get the average behavior of network characteristics. All these different simulation types are illustrated with a simplified example network that has 3 branching levels (Fig 2).

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 2. Comparison of real data for vascular networks versus random simulations of branching junctions.

The real and simulated networks (via local to global spatial constraints) are separated by different rows. A schematic small network is given to describe how different simulations are performed. The vessels and the fixed endpoints of the real branching network are represented in red. Vessels that result from random branching simulations are in black. The healthy mouse lung network and the simulated mouse lung networks are shown within a minimum spherical boundary that contains all branching data from the real network. Here, the red nodes for each figure correspond to the real data, whereas the black nodes correspond to the simulated data. Note that the terminal tips and the most upstream node (i.e., the source) are determined from real data and fixed throughout all simulations. The resulting asymmetry ratio distributions for length and branching angles are provided for the real network and for each of the simulations. The statistical comparisons of random branching simulations with empirical data are given in Table 1.

https://doi.org/10.1371/journal.pcbi.1005223.g002

For the fully local case, we randomly place a branching junction inside each triangle of unshared endpoints of three connected vessels, corresponding to the same spatial restriction as in our optimizations above (Fig 2). This random branching model would likely correspond to branchings determined by locations of highest shear stress [32] or local maxima in gradients of vascular endothelial growth factor (VEGF) [30] that signal branching at local and intermediate scales. For the fully global case, we only require that the branching junctions are randomly positioned within a minimum sphere that contains all nodes from the real data [55] (Fig 2) and that the network terminates at the appropriate (most extreme upstream and downstream) endpoints in the network.

For intermediate degrees of spatial constraint between the fully local and fully global cases, we consider two possibilities. These intermediate states are constructed such that they involve sequential updates of the branching junction. The first intermediate randomly branching network simulation (intermediate 1) starts by randomly positioning a branching junction within the local triangle corresponding to the most upstream vessel (i.e., source) of the network and the endpoints of its two daughter vessels. Based on the new location of this branching junction, the endpoints of the daughter vessels are then updated and used to define new triangles in the next step in which the daughters become the parents. The simulation continues this updating process by working down through the network until it reaches the terminal tips (Fig 2, S4 Fig). Notably, this simulation leads to branching junctions that are approximately confined to the plane. The other intermediate randomly branching network (intermediate 2) starts with the terminal tips (most downstream vessels of the network) and builds backwards to the first branching node (i.e., the source). We assign the position of each branching junction by creating a spherical boundary around the two fixed downstream endpoints (i.e., V₁, V₂, Fig 1D) such that the center of the sphere is at the midpoint of V₁V₂ and the sphere has a radius of the length |V₁V₂|. We then randomly position the branching junction at a point that can occur anywhere within the volume of this three-dimensional sphere. Thus, for this simulation a branching junction may not always fall within the plane defined by the vessel endpoints, as it does for the first simulation for intermediately-constrained random branching. For this case, the upstream endpoint (i.e., V₀) does not affect the location of the branching junction, reducing the degree of locality compared with the constraint for the simulation for intermediate 1.

Analysis of asymmetry patterns in vascular data

To characterize the branching asymmetry of the vascular structure, we calculate asymmetry ratios between siblings. In particular, the asymmetry ratios for radius and length are λ_r = r₁/r₂ and λ_l = l₁/l₂, respectively, where the value of the sibling with the larger radius or length is always chosen to be in the denominator (r₁ ≤ r₂, l₁ ≤ l₂) [21, 27, 56]. Thus, a ratio of 1 indicates perfect symmetry, whereas smaller values indicate more asymmetrical branching. We further provide a similar measure to quantify local asymmetry in sibling branching angles as λ_θ = θ₁/θ₂ again with θ₁ ≤ θ₂ (Fig 1C).

It is easy to see that asymmetry in vessel lengths is related to asymmetry in branching angles. When the downstream ends of the daughter vessels are equidistant from both the upstream and downstream ends of the parent vessel, siblings have identical length and branching angles, resulting in symmetry λ_l = λ_θ = 1. However, if daughter vessels are not equidistant from the upstream end, even symmetric sibling lengths can result in asymmetric branching angles or vice versa. Therefore, the value of λ_θ represents a combination of the asymmetry in lengths and the asymmetry of the alignment of daughter vessel endpoints.

After executing Angicart on high-quality tomographic images, we quantified the extent to which the analyzed networks are asymmetric by plotting the frequency distributions of λ_r, λ_l, and λ_θ (Fig 3A and 3B). Even though asymmetric radius branching exists, data for both mouse and human show distributions of λ_r that are skewed towards 1, meaning that the radius ratio is skewed towards perfect symmetry. In contrast, the length asymmetry ratio (λ_l) is distributed almost uniformly, suggesting a high degree of asymmetry in sibling lengths. The frequency distribution of branching angle asymmetry (λ_θ) is biased towards the right, corresponding to perfect symmetry in both networks and similar to results for asymmetry in radii (statistical calculations are given in Table 1). That is, a parent vessel tends to branch into sibling vessels that are separated by equal branching angles. Intriguingly, the similarity of the pattern observed in radius asymmetry might suggest a correlation between the radius and the branching angle, providing motivation to investigate how branching angle depends on radius.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 3. Histograms or frequency distributions of the asymmetry ratios for radius (λ_r), length (λ_l), and branching angles (λ_θ) of vascular networks.

(a) mouse lung (1 individual) and (b) human head and torso (18 individuals). Note that radius and branching angle asymmetry ratios are both skewed towards perfect symmetry, whereas the length asymmetry ratio shows no skew and reveals much more asymmetry. (c) Histograms of branching angles for combined data of human and mouse networks appear to be unimodal both for θ₀ and for θ₁ & θ₂ with peaks at 1.51 and 2.21 radians, respectively.

https://doi.org/10.1371/journal.pcbi.1005223.g003

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Table 1. Statistical comparison of material-cost (MC) optimizations and random spatial constraints with empirical data.

https://doi.org/10.1371/journal.pcbi.1005223.t001

In addition to the asymmetry ratios explored above, we plot the histogram for the raw data on branching angles. Human and mouse networks show similar patterns, so we combine data for these two networks in our histogram plot (Fig 3C). Analyzing the branching angle between the sibling vessels (θ₀) and the branching angle between the parent and daughter vessels (θ₁ and θ₂) separately, we find unimodal distributions that peak at 1.51 and 2.21 radians, respectively. This shows that planar branching [57, 58] with orthogonal daughter vessels is frequent across the networks. Individual branching angle plots for mouse and human networks are given in S5A and S5B Fig.

Altogether, these network-level patterns for vessel radius, length, and branching angles hold across 18 different human subjects, different species (human and mouse), and different tissues (head and torso versus lung). Moreover, the radius and length asymmetries are consistent with findings in plants [45–47]. These results suggest that very general and ubiquitous selection pressures and developmental processes may shape the architecture of the vascular system across taxa, from humans to mice to plants, as well as across tissues, from lungs to head and torso.

Optimal branching patterns for material-cost (MC) optimization

In this section, we introduce the MC optimization results that include surface-area (MC-SA) and volume (MC-V) constraints. Here, we only focus on non-degenerate branching solutions—solutions that do not collapse to a vertex—and compare those with real data. The fraction of the non-degenerate and degenerate cases is provided in the S6 Fig.

Network-level comparison

Taking the radii and vessel endpoints from our real human and mouse vascular networks, we use the material-cost (MC) optimization solution provided above (Eq 1) to compute the optimal branching angle and length asymmetry at each branching junction. To compare predicted values to the real network structures presented in the previous section, we plot the distributions of λ_l and λ_θ as well as the raw branching angle distributions for θ₀ and θ₁ & θ₂. We find that all these properties are visually consistent across both networks and for both volume and surface-area constraints.

In particular, distributions of optimal λ_l are close to uniform, whereas optimal λ_θ distributions are skewed towards perfect symmetry (Fig 4A and 4B). Both of these match the general patterns of asymmetry in branching observed in Fig 3. However, the degree of skewness in optimal λ_θ is sharper than the real λ_θ distributions, especially for the MC-SA. The statistical analysis of all these plots (including mean, median, skewness) are presented in Table 1.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 4. Histograms or frequency distributions of optimal asymmetry ratios for length (λ_l) and branching angle (λ_θ) derived from material-cost (MC) optimizations.

Surface-area (MC-SA) results are shown as solid lines and volume (MC-V) results are shown as dashed lines for (a) mouse lung and (b) human head and torso.

https://doi.org/10.1371/journal.pcbi.1005223.g004

Next, we plot optimal θ₀ and θ₁ & θ₂ histograms for the combined datasets (Fig 5). We observe that optimal calculations yield unimodal distributions as in histograms of the real data. Optimal θ₀ shows a mode around 1.79 radians, hence the optimal calculations are shifted towards larger values compared to the actual θ₀ distribution that have a mode at 1.51 radians. In contrast, the peak for the uniform distribution θ₁ & θ₂ at 2.24 radians almost matches the actual peak at 2.21 radians. The separate figures for human and mouse networks for each constraint are given in the S5 Fig. The full statistical analysis of branching angle histograms for the individual networks and the combined datasets are provided in the S3 Table.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 5. Histogram of optimal branching angles for combined data of human and mouse networks for material-cost (MC) optimizations.

All histograms appear to have unimodal characteristics both for θ₀ and for θ₁ & θ₂ with respective peaks at (a) 1.79 and 2.25 for the surface-area constraint and (b) 1.79 and 2.24 for the volume constraint.

https://doi.org/10.1371/journal.pcbi.1005223.g005

Junction-level comparison

The previous section compares the network-level patterns of the optimal calculations and the empirical data. Here, we provide comparisons at local branching junctions for the MC optimizations. We plot actual versus optimal branching angles and calculate the Pearson correlation coefficient, which would be 1 if our predictions were always perfect. From the plots of branching angle for MC-V (Fig 6) and MC-SA (S7 Fig) optimizations, we find that the predictions and the empirical data are weakly correlated (p-values <0.01). Moreover, we see that the volume constraint yields better agreement with data than the surface-area constraint at the junction-level. However, the correlation coefficients for both constraints indicate that the predicted optimal branching angles are a weak predictor of the actual branching angles at the local junction level. Consequently, it seems our theory needs further refinement or replacement.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 6. Junction-level comparison of optimal versus actual branching angles for the volume constraint of material-cost optimizations (MC-V).

(a) mouse lung and (b) human head and torso. The Pearson correlation coefficients and p-values are calculated for each plot.

https://doi.org/10.1371/journal.pcbi.1005223.g006

Analytical solutions for power-cost optimization beyond single branching (PC-1)

Now, we introduce the results for the solution of the PC-1 optimization. As described in the section presenting the PC-1 optimization scheme, we derive the approximate solution by considering a limiting case of the downstream impedances c₁ and c₂. With this method our solution predicts regions in the c₁c₂-plane separated by lines over most of the range of values. The regions in the c₁c₂-plane correspond to combinations of values that lead to branching geometries that are categorized as follows: collapse to daughter endpoint, collapse to parent endpoint, or no-collapse (i.e., a non-degenerate branching point).

Using the vessel endpoint and radius information from the real data, we predict the linear equations that form the boundaries between these regions, and we label them according to the categories of solution described above (Fig 7). As a further check, we compute numerical solutions at each discretized point in the c₁c₂-plane and mark different categories of solution by different colors. We examine two examples, corresponding to symmetric (Fig 7A and 7B) and asymmetric (Fig 7C and 7D) parameter values. Both show that our approximate analytical solution matches the numerical solution extremely well. On closer inspection, it is clear that the different regions and categories are not separated via purely linear functions across the entire plane, but instead the boundaries are curved for small values of c₁ and c₂ (Fig 7B and 7D). This result reveals a mismatch between the analytical and numerical solution at the smallest scales, i.e., vessels close to the capillaries, which is exactly where our solutions should fail based on the limits of our approximation scheme.

Download:

• PPT
  PowerPoint slide
• PNG
  larger image
• TIFF
  original image

Fig 7. Comparison of approximate solutions with numerical solutions for the PC-1 (power-cost (PC) optimization beyond single branching).

Approximate solutions define linear boundaries on the c₁c₂-plane between different categories of the solution space: collapse to daughter endpoint, collapse to parent endpoint, and no-collapse. The different categories calculated from numerical simulation are marked by different colors as indicated in the figure. (a) An example of symmetric branching in vessel radius with parameter values: |V₀V₁| = |V₀V₂| = |V₁V₂| = 1, r₀ = 1.20, r₁ = 1, and r₂ = 1, where c₁ and c₂ take values in the range [0,20]. (b) Zoomed version of (a) into the plane [0,2] × [0,2] with the same resolution as in (a). (c) An example of asymmetric branching in vessel radius with parameter values: |V₀V₁| = 0.8, |V₀V₂| = |V₁V₂| = 1, r₀ = 1.1, r₁ = 0.85, and r₂ = 1, where c₁ and c₂ take values in the range [0,20]. (d) Zoomed version of (c) into the plane [0,2] × [0,2] with the same resolution as in (c).

https://doi.org/10.1371/journal.pcbi.1005223.g007

Network-level results of randomly branching networks with local to global constraints

Lastly, we compare our empirical results for network-level characteristics of branching angles and length asymmetries to results from simulated randomly branching networks with local to global constraints. Because mouse and human networks yield similar results, we only present results for the mouse in the main text. Random simulation results for the human subjects are given in S8 Fig.

As in the empirical data, random simulations with local to intermediate constraints yield uniform distributions for asymmetric length ratio (λ_l), whereas the branching-angle asymmetry ratio (λ_θ) is skewed towards symmetry (Fig 2). In contrast, the fully global constraint generates branching networks with distributions skewed towards symmetry for both λ_l and λ_θ, inconsistent with empirical data (Fig 2). Further statistical comparisons with empirical data and optimal branching results are presented in the next section.

Comparison of optimal branching, random branching, and empirical data

To provide a rigorous comparison across different MC optimizations (surface-area and volume) and spatial constraints (local to global) with the empirical data, we compute the statistical properties—mean, standard deviation (SD), skewness and standard error for each—of the resulting asymmetry distributions in Table 1. Additionally, we use the Kullback-Leibler divergence measure (KL) to quantify the distance between the empirical asymmetry distributions and the optimization or random-simulation asymmetry distributions. We determine p-values by performing bootstrap samples up to half the size of the real data [59]. By definition, a p-value is equal to 1.00 for a comparison of the real data with itself because that implies an exact match in the distributions.

There are several conclusions based on these results. First, the random branching simulations that are globally constrained do not produce results that are statistically similar with mouse lung or human head and torso in terms of length asymmetry. In addition, the intermediate 2 constraint on random branching—where spheres are used to determine the branching locations, Fig 2—poorly matches with real data in terms of both length and branching angle asymmetries for the human head and torso network.

Except for the global constraint, results for the mouse lung network reveal that the random branching simulations perform as well as the MC optimizations in terms of the general characteristics (i.e., the first few moments—mean, SD, skewness) of the distributions. In contrast, for the human head and torso, the MC optimizations yield overall better agreement with the real data than the random simulations, especially for length asymmetry. This finding is consistent with the KL significance test. Based on the KL p-values, we observe that all of the random branching simulations do not do a good job of matching the real data for human networks. However, it suggests that vascular branching derived from MC-SA does reasonably well at recreating length asymmetry patterns that are similar to real human data. In contrast, the mouse lung compares more favorably with the random-branching simulations at local and intermediate scales than it does with the MC optimizations (Table 1).

Processing of angiographic images, vessel extraction with Angicart software, and resolution of data

In our study we analyze the cardiovascular structure of mouse lung and human head and torso by processing the three-dimensional stacks of images via the software Angicart [29]. Mouse lung images are collected through the microtomography (micro-CT) [44], whereas human head and torso images from 18 different subjects are obtained through MRI [29, 43]. The detailed image acquisition for each dataset are given in [44] and [29]. Here, we present how these images are processed to extract the vascular geometry. Before processing the tomographic images in Angicart, the images are first downsampled to reduce the noise and decrease the processing time. Angicart relies on the manual input of the intensity threshold parameter between 0 and 1 to identify the set of voxels belonging to the vascular network. For each set of data, the intensity threshold has been chosen by the visual inspection of blood vessels, conducting the sensitivity analysis as well as identifying the percentage of the possibly misclassified vessels as presented in Newberry et al. [29]. All these steps are done to smooth or reduce the noise in the data. The details on the preprocessing steps, set of thresholds, as well as the version of the software in regards to these criterions are given in the Supplementary Material. Performing Angicart on the angiographic images extracts the branching network topology with the vessel characteristics such as positional coordinates, radius, length, number of children. The output of the Angicart for each dataset is available in S1 Dataset.

We also note that different imaging modalities can lead to different imaging qualities and that the level of spatial scales that they can identify varies. Because micro-CT provides higher resolution images than MRI, the vessel sizes identified in the mouse lung are substantially smaller than the vessel sizes identified in the human head and torso data (>10 microns versus >1mm). These differences have the advantage of allowing us to investigate the branching geometries for large versus small vessels as well as pulsatile versus viscous dissipation flow regimes [4, 27, 35]. Moreover, the difference in the quality of different medical imaging techniques inevitably affects the amount of noise in the data. Possibly misidentified vessels from the Angicart output are characterized by defining non-deformed vessels—more than 10% of the volume of voxels of the vessel lies inside a distance of radius+1 voxel from the centerline of the corresponding vessel. We find that the fraction of carefully classified, non-deformed vessels is greater in the mouse lung data compared to the human head and torso data (S2 Table), meaning it is of higher quality on average, as expected from the image resolution properties of micro-CT and MRI. This realization makes it especially noteworthy that the mouse lung data exhibits the strongest signal of random branching, which we can confidently say is not due to effects of noise but to actual branching patterns (see Results, Comparison of optimal branching, random branching, and empirical data).

Computing branching angles from extracted vessel skeleton

To calculate our three branching angles (Fig 1C) using our empirical vessel measurements, we first must specify the three intersecting lines that define the angles. One endpoint is common to all three lines, representing the branching junction (J) that connects the parent to the two daughter vessels. Therefore, the remaining choice is in how to define the other endpoint of all three vessels. In principle, this could be done using points along the vessel that are very close to the branching junction (J), the midpoint between the branching junction and the other vessel endpoint, or the exact endpoint of the other vessel. For a perfectly straight vessel, this choice will have no effect on the computed branching angle, but for curved or tortuous vessels, which are common in real vascular systems, these choices will lead to different values for the branching angles. For this study, we choose the exact endpoints of the vessels to define the straight lines that define the angles (S9 Fig). We do this for three reasons. First, we relate branching angles to vessel lengths, which are defined relative to the endpoints, so these relations will be most faithful if we use lines corresponding to the full length of the vessel. Second, we argue there is more of a constraint on the endpoints of vessels than the exact path they take to reach those points, which may include more developmental stochasticity. Third, there is no arbitrariness to the choice of vessel endpoints, whereas choosing 1 pixel versus 5 pixels away from the branching junction is more subjective. Using our choice of lines and endpoints, we calculate the magnitude of the angles between the straight lines defined by the positional coordinates—the endpoint coordinates of the vessels V₀, V₁, V₂, and J—at the bifurcation as shown in Fig 1C and 1D.