The role of vascular complexity on optimal junction exponents

We examine the role of complexity on arterial tree structures, determining globally optimal vessel arrangements using the Simulated AnneaLing Vascular Optimization algorithm, a computational method which we have previously used to reproduce features of cardiac and cerebral vasculatures. In order to progress computational methods for growing arterial networks, deeper understanding of the stability of computational arterial growth algorithms to complexity, variations in physiological parameters (such as metabolic costs for maintaining and pumping blood), and underlying assumptions regarding the value of junction exponents is needed. We determine the globally optimal structure of two-dimensional arterial trees; analysing how physiological parameters affect tree morphology and optimal bifurcation exponent. We find that considering the full complexity of arterial trees is essential for determining the fundamental properties of vasculatures. We conclude that optimisation-based arterial growth algorithms are stable against uncertainties in physiological parameters, while optimal bifurcation exponents (a key parameter for many arterial growth algorithms) are affected by the complexity of vascular networks and the boundary conditions dictated by organs.

In living organisms, the bifurcation exponent associated with arteries and arterioles is often measured to deviate from three 6 , which is not fully understood, although several factors are known to lead to γ opt = 3 in single-vessel analyses. Reference 6 reviews a range of studies that measure γ in various organs and species, with γ ranging from ∼ 2.1 to ∼ 3.5 in systemic arterial trees (with a greater range in the pulmonary vasculature). There are two key lines of argument to explain why γ < 3 . The effects of pulsatile flow, elastic wall vessels, and turbulence all contribute to additional power dissipation in an arterial segment, leading to a reduction in the optimal junction exponent to γ opt = 2.33 6 . Alternatively, it has been argued that cross-sectional area is conserved at bifurcations, i.e. γ opt = 2 7 . Curiously, in some organs, γ is measured to be slightly greater than three 6,8 . To our knowledge, no explanation of this effect is available, since corrections to flow in single artery analyses to include turbulence, pulsatile flow, and elastic wall vessels, lead to γ opt < 3 . This suggests that single vessel analysis may be insufficient to explain γ opt > 3 and thus that complex networks representing large numbers of arteries may be needed to properly analyse vascular trees.
We propose that, in order to fully understand the optimal branching exponents in vascular trees, it is essential to take into account the complex structure of the entire arterial network in an organ, and the boundary conditions imposed by the organism on the arteries that enter organs (particularly on the flow and artery radius). A single vessel is part of a much larger arterial tree for an organ, that is in turn part of an organism, and the role of this additional complexity on optimal supply networks is poorly understood. Two factors define boundary conditions for arterial growth algorithms: (1) The metabolic demand of the organ determines the blood flow to the organ. (2) The radius of the primary artery supplying that organ is determined by a compromise between the whole organism and the organ.
There are a huge number of possible ways to connect arterioles to arteries, yet not all are optimal. Our goal is to examine the complex multiscale structures of vascular networks that emerge from optimisation considerations both computationally and analytically. The process of optimisation in complex trees can be difficult to reproduce computationally. When optimising, the space of possible connections between the capillaries and the arteries needs to be explored. The number of possible combinations of vessels associated with these connections is enormous, and all vessels contribute to metabolic and pumping costs. It is not possible to search through every combination for all but the smallest trees. Moreover, deterministic minimisation methods (e.g. steepest descent) run the risk of finding local rather than global minima. Thus, stochastic optimisation algorithms (such as simulated annealing) are needed 9 .
Stochastic optimisation algorithms are highly advanced optimisation approaches that use random variables to quickly search through the configuration space of a problem to find (near) globally optimal solutions 10 . In particular, simulated annealing is guaranteed to approach the global optimum as computational times are increased 11 . We previously introduced the Simulated AnneaLing Vascular Optimisation (SALVO) algorithm to find the globally optimal structure of arteries using simulated annealing to overcome these problems 8,9 . The power of the numerical SALVO algorithm is that the generated trees are globally optimised and therefore represent the lowest possible total power cost associated with the whole vasculature. While the globally optimised solution represents an idealised evolutionary endpoint, insight into the compromise associated with optimising the competing costs of the different metabolic requirements associated with complicated vasculatures can be gained by this analysis. A summary of SALVO can be found in the "SALVO" section, with a detailed introduction in Ref. 8 .
Our work goes beyond previous analyses 5,6,[12][13][14][15] , by optimising entire trees, rather than a single arterial bifurcation. The constraints on flow and radius of root vessels in real organs are also taken into account. Further to this www.nature.com/scientificreports/ analysis, we use the SALVO algorithm 8,9 to determine the globally optimal bifurcation exponent, which allows us to check any analytical expressions that we have calculated against numerically exact values. This paper is organised as follows: In the "Methods" section we introduce the methodology used in this paper. The "Analytical results" section presents analytical expressions for the properties of large vasculatures. The "Numerical results" section presents results from the SALVO algorithm applied to two-dimensional planes. Finally we present a "Discussion and conclusions" section.

Methods
In this section, we discuss some of the basic formalism related to the calculation of optimal γ , and its relation to the SALVO algorithm.
Power cost. In order to make it possible to calculate the properties of large and complex arterial trees, the arterial tree is divided into straight segments and bifurcations, such that segments join bifurcations. Thus each bifurcation becomes a node within a tree network. In the following, we shall interchangably use the terms node and bifurcation. The largest vessel from which blood arrives shall be referred to as the root node. The smallest vessels at the end of the tree shall be identified as leaf nodes. A schematic of the terminology associated with the arterial network can be found in Fig. 1b.
Poiseuille flow is assumed within each segment. It is also assumed that maintaining blood has a metabolic cost proportional to its volume. With these assumptions, the power cost for pumping blood through a single arterial tree segment is 5 , where j denotes a segment, r j the segment radius, l j its length, f j its volumetric flow, m b the metabolic power cost to maintain a volume of blood, and µ the dynamic viscosity of blood. The first term in this expression represents the metabolic cost to maintain a volume of blood, and has units of energy expended per unit time (power). The second term represents the power required to pump blood through the segment. The power cost associated with bifurcations is neglected.
The total cost, W , of an arterial tree is the sum of these individual segment costs, Equation 3 will be referred to as the cost function.
Murray's law. Murray's law ( f ∝ r 3 ) is derived by optimising total power expenditure for pumping blood through a single segment, as given by Eq. (2) 5 . By differentiating Eq. (2) with respect to r j , When ∂W j /∂r j = 0 , the optimal relation between f j and r j can be found. This leads to Murray's law, In the following analysis, we will assume that l = l root r α /r α root , where l root and r root are the length and radius of the root segment respectively, and α is the length-radius exponent 6 . This slightly modifies the preceding argument, so that, where f root is the flow through the root segment.
SALVO. The (SALVO) algorithm developed in earlier papers for three dimensional arterial trees 8,9 can also be used to generate arterial trees in two-dimensional planes. In this section, an outline of this algorithm is given. The algorithm is similar to the approach for growing cardiac and cerebral vasculature 8,9 , with some differences relating to the use of fixed nodes to supply tissue. A detailed description of the SALVO algorithm can be found in Ref. 8 .
In a key difference to previous work, we study an idealised two-dimensional (2D) piece of 'tissue' , with fixed leaf-node positions. The root node of the tree is fixed to the corner of a square region of side a(= 10 mm) . A Poisson disc process is used to place leaf nodes 16 . The whole 2D region is accessible by nodes, with only metabolic-and flow-related penalties in the cost function, Eq. (3). Since the leaf nodes are evenly distributed by the Poisson disc process, and do not move during the update process, a penalty for under-and over-supply (see Ref. 9 ) is not required. Also, no penalties are required to stop penetration of large vessels into the tissue (as in Ref. 8  www.nature.com/scientificreports/ cavities within the tissue (such as the ventricles of the heart in Ref. 9 ) since all vessels lie within a continuous 2D tissue. Overall this decreases the complexity of the algorithm. On each iteration, modifications to the binary tree are attempted by either (1) selecting a node at random and then moving it or (2) selecting two nodes at random and then changing the tree structure by swapping the parents of those nodes. These updates are sufficient to ensure ergodicity. Updates are summarised in Fig. 2 and relevant parameters are summarised in Table 1. The root node is never updated. An example of moving a node is shown in Fig. 2a. In the update, the initial configuration at the top of the figure changes to the final configuration at the bottom of the figure, by moving the location of the node labelled a. The distance moved is short (0.05 mm) in 30% of updates, and long (0.5 mm) in 20% of updates. In the version of the algorithm used in this paper, leaf nodes are never moved. An example of swapping the parents of a node is shown in Fig. 2b. In this case, it is the parents of nodes labelled b and d that are swapped. In the initial configuration, node c is the parent of node d,    www.nature.com/scientificreports/ and node a is the parent of node b. After the update, node a is the parent of node d and node c is the parent of node b. The swapping of parent nodes is attempted on 50% of iterations. We use simulated annealing to optimise the cost function 17 . Within this framework, acceptance of the updates is determined according to the Metropolis condition, is the change in cost associated with modifying the tree from the configuration in iteration (θ) to the configuration proposed for iteration (θ + 1) , and W as defined in Equation 3 is the cost function at the core of the SALVO algorithm. In practice, a uniform random variate r ∈ [0, 1) is calculated and if P > r the proposed configuration is accepted. If the proposed configuration is rejected, then the configuration of iteration θ is carried forward to iteration θ + 1 . T θ is the annealing temperature, which is slowly reduced using the common exponential schedule, T θ +1 = ǫT θ where θ is the iteration number, ǫ = exp(ln T 0 − ln T ν )/ν , ν the total number of iterations, and T 0 ( T ν ) are the initial (final) temperatures. Updates are iteratively applied until the final temperature is reached. The algorithm is summarised in the flow diagram of Fig. 2c.
As we have discussed elsewhere 8 , limitation of the method to a few thousand nodes is related to the combinatorial (factorial) growth of the number of possible tree configurations. 5000 nodes trees are already quite detailed, and determination of γ opt for such trees already offers a substantial computational challenge: The required 10 8 updates take approximately 11 hours on a single thread of a Threadripper 2990WX processor, and we have to distribute the optimisation of large numbers of similar trees with different γ across the full 64 threads of the processor to find a single value of γ opt (and are then calculating for many different physiological parameters). To make these calculations with much larger numbers of nodes, very significant advances in computational power would be needed. The primary issue is that the number of updates needed to optimise the tree grows with the number of possible combinations. Therefore, the optimisation of larger trees is not just a matter of greater computational power, since the required computational power grows so rapidly. A possibility is application of alternative optimisation algorithms, e.g. genetic algorithms, ant-trail optimisation etc, although these are much harder to implement for this problem. Unfortunately, even if certain optimisation algorithms are faster, all will eventually be limited by the combinatorial growth of possible tree configurations. In the context of this paper, we need to identify the global minimum, so approximate strategies such as multiscale algorithms 18 or constrained constructive optimisation (CCO) 19 are not appropriate.

Analytical results
In this section, we discuss analytic approximations to the total power of large and complex vascular trees. This provides initial insight into the deviations in γ opt caused by tree complexity. We start by discussing a formalism for simplifying the total power calculation of large and complex arterial trees.

Formalism and simplifications.
Arteries can be grouped together, so that each group comprises arteries with identical properties (e.g length, diameter, flow). In a real arterial system, this would not be exact, but it would still be possible to group arteries with similar lengths, radii, and flows together. By using this grouping, the total power can be rewritten as, where N(r j , l j , f j ) is the number of arterial segments with identical radii, lengths and flows.
Under the restriction that the flow in all leaf nodes is identical and equal to f leaf , the flow in each segment is, where n is an integer and represents the total number of leaf nodes downstream of the segment. Comparing Eq. (1) with flow conservation, a radius-flow relation is identified: Thus, by substituting Eq. (9) into Eq. (10), the radius can be rewritten in terms of n and γ: Experimental data suggest that the length of an arterial segment is proportional to a power of the radius, where the value of the exponent α is typically close to 1.0 6,20 . By substituting Eqs. (11) and (12) into Eq. (2), the power required to maintain blood flow through a segment is found to depend only on the flow f n , www.nature.com/scientificreports/ Thus, the dimensionless metabolic ratio, defined as � = m b π 2 r 6 leaf /8µf 2 leaf , the bifurcation exponent, and the number of nodes, N, define the parameter space. The power associated with a segment is then, where C = 8µf 2 leaf l leaf /πr 4 leaf . Both C and are defined in terms of the leaf node properties. A similar ratio for the root node, � root = m b π 2 r 6 root /8µf 2 root can be defined for convenient contact with experiment. The values r root and f root are often known from experiment, e.g. Doppler ultrasound, and N can be estimated. This ratio can be related to via The total power required to supply the whole vascular tree is, N n is the number of segments with flow nf leaf , and simplifies the function N(r i , l i , f i ) . For any tree structure, N is always the number of leaf nodes, so N 1 = N . There is always a single root node with total flow Nf leaf , so N N = 1 . No node has flow greater than Nf leaf , so N n>N = 0 . The remaining N n are dependent on the structure of the tree. At each bifurcation, flow conservation requires that n in = n out,1 + n out,2 , so n is an integer. Total power is linear in length scale, so the location of any minima in the power with respect to γ is independent of a. It is the global minimum with respect to γ that sets the structure of the tree, and when locating the minimum, ∂W/∂γ = 0 , so the factor of l leaf in C simply cancels, thus making the solution independent of a. Changes in r leaf can be absorbed into the ratio m b /µ and thus are similar to changing the metabolic requirements of the organ 9 .
There are two special tree structures: the fully symmetric and fully asymmetric trees. In the first case, identified as a fully symmetric tree, the flow is split evenly at each bifurcation. For the case which we shall identify as fully asymmetric, a single leaf node emerges at each bifurcation and the rest of the flow passes down the other bifurcation. We will explore these special cases in the following two sections.
Fully symmetric vascular tree. In a fully symmetric tree, all of the segments with flow n exist at the same bifurcation layer. Each layer, denoted by the integer m, has 2 m segments, where m is the number of bifurcations upstream of that layer ( m = 0 at the root segment). Within a layer, all segments have the same flow, and thus the same radius and length. The tree has a total of M layers, so 0 ≤ m ≤ M . Therefore N n = 2 m if n = 2 M−m , N n = 0 otherwise, and the total power cost in Eq. (16) becomes, By summing this geometric series, the total power cost for a fully symmetric tree is found to be, Fully asymmetric tree. The total power cost of the fully asymmetric tree may be calculated by noting that each discrete flow is represented once for all n, so N n = 1 , for n < N . The exception being that there are N leaf nodes so N 1 = N.
Substitution into Eq. (16) gives, So the total power cost for an asymmetric tree is where H (r) n is the generalized harmonic function, n k=1 1/k r .
Optimal bifurcation exponent. The optimal value of γ is obtained by numerically solving ∂W/∂γ = 0 for Eqs. (18) and (20) . Results are shown in Fig. 3 for symmetric and asymmetric trees, for various �, α and N.
The optimal bifurcation exponent is strongly dependent on the metabolic ratio, , which can change due to physiological boundary conditions on flow and radius at the input vessels. These constraints may be due to limits in the size of the largest vessel in the tree imposed by the physiology of the whole organism and the flow  www.nature.com/scientificreports/ demands of the organ. Figure 3a shows the optimal value of γ . When = 1 and α = 1 the result of Murray's law ( γ opt = 3 ) is recovered. γ opt is qualitatively unchanged by the structure of the tree. Results for asymmetric and symmetric trees with N = 2.047 × 10 3 follow essentially the same functional forms. The optimal bifurcation exponent for the asymmetric tree is closer to γ = 3 than the symmetric tree. Also shown in Fig. 3a are numerical values from SALVO, which will be discussed later.
The optimal bifurcation exponent is modified away from γ = 3 by changes in the length exponent, α (Fig. 3b). This structural effect potentially has implications for the value of γ opt in organs, since α can vary with organ type, with estimates ranging from 0.89-1.15. In practice, changes in γ opt for this variation in α are far smaller than the error for measurements of γ and changes in α can essentially be neglected.
Deviations from Murray's law are largest for small trees and strongly dependent on changes in the metabolic ratio. The larger the tree, the closer to Murray's law γ opt becomes. Fig. 4 shows variation of γ opt with N for fully symmetric trees. For vascular tree sizes of between 10 3 and 10 6 segments, which are typical in organs, γ opt ranges between 2 and 4.

Numerical results
The generation of globally optimal trees using a numerical algorithm helps to test analytic expressions, and provides additional morphological measures that can be used to understand arterial networks. In this section, we use SALVO to investigate the role of vascular complexity and physiological boundary conditions on the properties of globally optimal trees. Several properties of the numerically generated trees are investigated. We determine the dependence of globally optimal tree structures on and γ . Through examination of W , we compute γ opt for complex trees. For each value of γ and investigated, arterial trees with up to 5000 nodes were generated. Table 1 summarises the parameters used for the numerical calculations. We note that we carried out checks on convergence with a subset of trees by using an anneal schedule with a larger number of steps ( ν = 10 9 ), finding no major changes to the tree structure.
Tree morphology. There are three regions of the parameter space with qualitatively different tree structures, examples of which can be seen in Fig. 5. In the figure, the vessel widths are normalised to the root radius to improve visibility. Trees are generated for N = 100 and various and γ:    In the star regions ( γ 2, � ≫ 1 and γ 4, � ≪ 1 ), long leaf segments connect root and leaf nodes (see top right panels in Fig. 5). This is due to the domination of the n 3/γ −1 term (that represents metabolic maintenance of blood) for low γ , and the n 3/γ ′ −1 Poiseuille term for large γ . Thus, terms with small n (i.e. leaf nodes) are favoured. Trees in both star regions are very similar, which is not a coincidence, and can be explained by examining the structure of Eq. (15). When α ≈ 1 , the power in a segment is W n = Cn(�n 3/γ −1 + n 1−3/γ ) . For γ > 3 , the exponents (which involve 3/γ − 1 ) have opposite sign to those for γ < 3 . So after the substitutions � = 1/� ′ , γ = 3γ ′ /(2γ ′ − 3), C ′ = �C , W n = C ′ n(� ′ n 3/γ ′ −1 + n 1−3/γ ′ ) , and the sum has an equivalent structure. The substitution is determined by identifying where 1 − 3/γ = 3/γ ′ − 1 . Since the prefactor C ′ scales the entire sum, then the minima of W and thus the results for γ , � and γ ′ , � ′ are identical. This symmetry is only approximate if α = 1.
Trees in the asymmetric regions ( γ 2, � ≪ 1 and γ 4, � ≫ 1 ), have a highly asymmetric structure, with long trunks snaking through leaf node sites (see top left panels in Fig. 5). This is due to the domination of the n 1−3/γ term due to Poiseuille flow for low γ , and the n 1−3/γ ′ metabolic cost term for large γ . Thus terms with large n (i.e. thick trunks) are favoured. A similar argument to that given for the star regions explains why the trees in both asymmetric regions have very similar structures. To quantify the effect of varying γ and on the network structure, we have examined average segment length, path length, radius asymmetry and Hausdorff dimension (Fig. 6). Average length is defined as l = l j /N . The average summed path length from root to leaf node is L = � path l j � . Radius asymmetry is measured using �r c> /(r c< + r c> )� (where at each bifurcation the larger and smaller of the radii of child vessels are labelled r c> and r c< , such that r c> ≥ r c< ). The Hausdorff dimension is calculated using a box counting method.
In the physiological region, the dominant factor controlling morphological properties is γ . Within that range, morphological properties do not depend strongly on . Average segment length is short and path length is long in this region, consistent with the branching structures seen for intermediate γ in Fig. 5. Bifurcation symmetry is in the range 0.58-0.62, so bifurcations are moderately symmetric. Although leads to minor changes in tree morphology in this regime, we note it can affect γ opt and thus the tree morphology via γ as a secondary effect.
In the star and asymmetric regions, is responsible for large variations in the tree morphology, and γ can also produce large variations in the various morphological and structural properties of the tree. Path length drops outside this region to approximately a/ √ 2 consistent with a large number of straight paths from the root node to leaf nodes. The radius asymmetry is increased in the asymmetric region relative to the physiological region, and is decreased in the star region. For all other regions of the parameter space, the asymmetry drops. www.nature.com/scientificreports/ Morphological measurements are not strongly dependent on changes in N, consistent with additional segments adding more detail to the tree, but not qualitatively changing the tree structure. Panels on the left of Fig. 6 show results for N = 2163 and panels to the right for N = 3968.
The Hausdorff dimension, d Haus , is also calculated (lowest panels of Fig. 6). This sits in the range between d = 1 and d = 2 . The dimension of the trees with the physiological range 2 < γ < 4 are lower than the spidery trees found outside this range.
A power law, l = Ar α , is found to relate the median segment length calculated using SALVO to the segment radius. Figure 7a shows data gathered from trees with 5000 > N > 2000 , 2.75 < γ < 3.25 , = 0.9 . Since there are many segments in the tree, and many trees in the analysis (to improve statistics), there is a distribution of segment lengths corresponding to each radius, just as found in experiments 6 . To indicate the spread of this distribution, we shade the range of lengths that sit between the 25th and 75th percentiles of the length distribution in light blue. To determine the parameter, α , we fit l = Ar α to the median value, finding the exponent to be α = 0.887 ± 0.088 , consistent with experimental values 6 .
The length-radius exponent, α , is consistent with experimental values for trees grown with 2.5 < γ < 3.5 , but can become effectively negative when long leaf segments start to dominate outside this region (Fig. 7b). To calculate the length-radius relation, segments are binned from trees with N > 2000 and specific γ and values before fitting a power law. Where exponents are negative, the relation only poorly follows a power law, and errors on α are large. The power law relation is well followed within the region 2.5 < γ < 3.5 , and this leads to smaller error bars. Overall, errors on α determined from fitting the power law are relatively large.
Optimal bifurcation exponent. The optimal bifurcation exponent γ opt can be determined without ambiguity from the minimum in W . Figure 8a shows how the total power cost varies with γ . There is a clearly defined global minimum for all values of shown. γ opt can be found by fitting a quadratic form to the bottom of the minimum.
The variation of γ opt with and N, numerically determined using SALVO, is qualitatively similar to the results from analytic expressions. Numerical values of γ opt for various values of vs N are shown in Fig. 8b, and compare favourably to Fig. 4. Several numerical values are compared with the analytic results in Fig. 3a, also showing good agreement for both symmetric and asymmetric trees.
We note that the location of the minimum in the cost function is very stable to changes in the random number seed, which determines the random number sequence used both in the Poisson disc process that initialises the leaf nodes, and as part of the Monte Carlo algorithm at the heart of simulated annealing. In practice a different seed was used for each point in Fig. 8a with no discernible fluctuation in the curves.

Sensitivity to uncertainty in physiological parameters. Uncertainty in physiological parame-
ters. There are two free physiological parameters that act as input to SALVO, and γ . In this section we discuss the uncertainty in these values.
Large uncertainties on the value of � = m b π 2 r 6 leaf /8µf 2 leaf , dominate errors in γ . Experimentally measured values for γ in systemic arterial trees (collated in Ref. 6 ) range from γ ∼ 2.1 to ∼ 3.4 . Fractional errors on these values are typically of order 15% for vessels of size > 1 mm and (significantly) less than 7% for vessels of size < 1 mm. The uncertainty in is itself dominated by uncertainty in the value of m b . Estimates of the parameter m b can vary by up to a factor ∼ 2 22 (corresponding to ∼ 40% fractional error). We expect radius and flow measurements to be significantly more reliable than the m b estimate. Thus we estimate to also vary by a factor 2 due to experimental uncertainty (or ��/� ∼ √ 2 − 1 ∼ 0.41). www.nature.com/scientificreports/ Uncertainty in tree morphology. Morphological properties depend on only two parameters, and γ . The uncertainties in morphological parameters can be related to variations in γ and in the usual way as, �O/O ∼ |dO/dγ | 2 (�γ /γ ) 2 + |dO/d�| 2 (��/�) 2 , where O represents one of the four calculated morphological properties. However, since the fractional uncertainties in measured γ values are much smaller than those in , we calculate the sensitivity based on Figure 9 shows the sensitivity of morphological properties to uncertainty in , for ∼ 1 . To estimate |dO/d�| we used the difference between morphological measures using values of = 0.9 and = 1.1 to estimate �O/�� , and display the average of these two measures in the figure with error bars representing the uncertainty in O due to . None of the morphological properties are very sensitive to the uncertainty in . There is little uncertainty within the physiological range 2 < γ < 4 . The largest uncertainties are found outside that range, but they do not lead to the possibility of qualitative changes to conclusions. www.nature.com/scientificreports/ Uncertainty in optimal bifurcation exponent. The key result of this paper is the optimal bifurcation exponent. This is only dependent on , so we estimate that the uncertainty of γ opt is, We estimate uncertainties of γ opt due to uncertainties in to be 10% . Figure 10 shows the sensitivity of γ opt to the uncertainty in . For tree sizes of N = 5000 , uncertainty in γ opt is approximately 6%, rising to ∼ 10% for the smallest trees of ∼ 100 segments. Thus the results are robust against the difficulties of estimating metabolic constants.

Discussion and conclusions
In this paper we determined analytic expressions, and carried out numerical calculations, for the properties and structures of globally optimal vascular trees, with the aim of understanding how overall complexity and physiological boundary conditions contribute to the optimal junction exponent and other properties of arterial trees. Analytic expressions were derived for the special cases of maximally symmetric and asymmetric arterial trees. The parameter space of the arterial trees was explored further by making numerical calculations with SALVO, enabling globally optimal vasculatures to be found for arbitrary tree morphology. The dependencies of tree structures, morphological properties, and optimal bifurcation (junction) exponent on physiological parameters are calculated.
The analytic expressions derived here are consistent with numerical calculations, and predict that γ opt does not vary strongly with tree symmetry, so we propose that the analytic expressions derived here are applicable to a wide range of vasculatures. Analytic expressions can be used for much larger trees than numerical optimisations, and would, therefore, be useful for predicting the properties of vasculatures within a range of organs where the number of vessel segments and overall complexity exceed the capabilities of current computers. We expect that it will be possible to extend the analytic expressions to include pulsatile flow and turbulence, and will investigate this possibility in future studies.
We predict that tree complexity is a significant contributor to the bifurcation exponents measured in living organisms. The deviations originating from tree complexity are of similar size to those predicted by including turbulence and pulsatile flow in previous analyses. These deviations are particularly significant if physiological boundary conditions lead to = 1 . This may occur since all organs, with their dramatically varying demands, are connected to the same major vasculature: so flow and radius associated with the arteries supplying organs may reflect compromise within an organism. We expect that large variations of γ opt with increasing complexity will also occur if a more detailed analysis including pulsatile flow and turbulence is carried out.
We predict that arterial tree complexity can lead to optimal bifurcation exponents, γ opt > 3 , a situation which can be found in experiment, and is of interest since inclusion of turbulence and pulsatile flow in single artery analyses leads to γ opt < 3 . Large values of γ are measured in e.g. the brain vasculature ( γ = 3.2) 8 , retina ( γ = 3.1 23 , γ = 3.9 ± 0.12 24 ) and other mammalian vasculatures where γ can range as high as 4 6 . Such large γ are not predicted by single segment analyses including effects related to pulsatile flow, elastic vessel walls and turbulence ( γ = 2.3) 6 . Tree complexity and organism imposed boundary conditions provide an additional contribution that can account for larger values of γ opt .  Figure 10. Sensitivity of γ opt to uncertainty in is found to be less than 10% for various and N. www.nature.com/scientificreports/ We predict that tree structures within the physiological regime, 2 < γ < 4 , are weakly dependent on all parameters except γ ; outside the physiological regime structures also depend strongly on ; and for all regimes tree structures are independent of N. Changes in N do not qualitatively change the morphology of the tree, but add more detail. Outside the regime 2 < γ < 4 , structure can change dramatically with .
For 2.5 < γ < 3.5 , we find length exponents in our computational trees that are consistent with the value α ∼ 1 obtained experimentally. Experimental values range from 0.85 < α < 1.21 6 . We find a similar range of values in our numerical calculations, and with improved description of the flow, the accuracy of the predictions could be improved. Values of α could be useful as input to other calculations.
Accurate values of γ opt are particularly relevant to computational techniques used for growing very large arterial trees in-silico, such as constrained constructive optimization (CCO). Such algorithms rely upon a fixed bifurcation exponent to set the radii in the generated trees 19,25,26 . Similarly, allometric scaling arguments require knowledge of γ 7 , and variations of γ opt could modify such approaches. γ opt is quite hard to measure experimentally, and we consider the calculation of such values to be a useful application of our technique.
Future work to include additional physics, such as pulsatile flow, turbulence and vessel elasticity, would lead to a computational model with enhanced predictive power. These improvements to the treatment of flow through vessels could be incorporated into both the analytic expressions derived in this paper, and into the cost function of SALVO without having to change the core algorithm. Once analytical expressions are modified to include this additional physics, we suggest that parameters such as m b could be determined from empirical results. The significant structural changes visible at γ ∼ 2 and γ ∼ 4 would also be interesting areas for further study, since the rapid changes in the tree morphology are reminiscent of a phase transition. These changes are on the edge of the physiologically relevant regime. Confirmation of a phase transition would require the identification of an order parameter and the signatures of critical behaviour.
Finally, we hypothesise that evolutionary compromises may favour closer adherence to the predictions of single segment analyses in organs with large flow demands to the detriment of less flow-hungry organs. Additional studies could be carried out to test this hypothesis. Overall, the computational and analytical approaches introduced here lead to a range of predictions regarding the structures of vascular trees, that provide interesting links to experimental and theoretical approaches.
Fundamental biophysical understanding of complex vascular structure has applications to modelling of cardiovascular systems and diseases. Computational techniques and analytic expressions for describing complex and multiscale networks have potential applications in computer modelling of physiology 8,9,18 , medical imaging 18 and diagnosis of cardiovascular disease 27 . Beyond the desire to understand the fundamental biological properties of vascular networks, deviations from optimal flow conditions could be a sign of underlying disease 27 . Another application of arterial growth algorithms is the correction of gaps in computed tomography and magnetic resonance angiography on small length scales 18 .

Data availability
The datasets generated and analysed during the current study are available in the ORDO repository: https ://doi. org/10.21954 /ou.rd.12220 490.