The role of vascular complexity on optimal junction exponents
TThe role of vascular complexity on optimal junctionexponents
Jonathan Keelan and James P. Hague School of Physical Science, The Open University, MK7 6AA, UK * [email protected]
ABSTRACT
We examine the role of complexity on arterial tree structures, determining globally optimal vessel arrangements using theSimulated AnneaLing Vascular Optimization (SALVO) algorithm, which we have previously used to reproduce features ofcardiac and cerebral vasculatures. Fundamental biophysical understanding of complex vascular structure has applications tomodelling of cardiovascular diseases, and for improved representations of vasculatures in large artificial tissues. In order toprogress in-silico methods for growing arterial networks, we need to understand the stability of computational arterial growthalgorithms to complexity, variations in physiological parameters such as tissue demand, and underlying assumptions regardingthe value of junction exponents. We determine the globally optimal structure of two-dimensional arterial trees; analysingsensitivity of tree morphology and optimal bifurcation exponent to physiological parameters. We find that, for physiologicallyrelevant simulation parameters, arterial structure is stable, whereas optimal junction exponents vary. We conclude that thefull complexity of arterial trees is essential for determining the fundamental properties of vasculatures. These results areimportant for establishing that optimisation-based arterial growth algorithms are stable against uncertainties in physiologicalparameters, while identifying that optimal bifurcation exponents (a key parameter for many arterial growth algorithms) aresensitive to complexity and the boundary conditions dictated by organs.
Vascular systems are highly complex, multiscale systems, with dominant physics that changes with length scale. Within atypical organ, vascular trees connect major arteries of ∼ ∼ − µ m. The physics of large vessels is often dominated by pulsatile flows and turbulence, while small vessels aremicrofluidic . Evolution must account for this complex and multiscale physics when optimizing arterial networks, sinceefficient vasculatures are crucial for supplying oxygen to tissues.Computational techniques and analytic expressions for describing these complex and multiscale networks have applicationsin physiology, tissue engineering, and medical diagnosis. Beyond the desire to understand the fundamental biological propertiesof vascular networks, deviations from optimal flow conditions could be a sign of underlying disease, and vascularization is akey issue limiting the growth of large engineered tissues.The challenge is that the complex and multiscale structures of vascular networks are difficult to reproduce in-silico .Within organs, there can be hundreds of thousands of arterioles dependent on every major artery. The number of possiblecombinations of vessels associated with these connections is enormous. Since it is not possible to deterministically searchall these combinations for all but the smallest trees, stochastic methods are needed. We previously introduced the SALVOalgorithm to find the globally optimal structure of arteries using simulated annealing to overcome these problems .Vascular networks are primarily constructed from bifurcations , which can be characterized by defining a bifurcationexponent (also known as radius exponent and junction exponent), γ . The radii of the two output vessels, r out , A and r out , B , arerelated to the radius of the input vessel, r in , via, r γ in = r γ out , A + r γ out , B . (1)Murray carried out a single-vessel analysis, which when combined with conservation of flow, shows that when flow invessels is approximately laminar (Poiseuille flow) the optimal junction exponent, γ opt = . In Murray’s analysis, there are twocompeting contributions to the metabolic demand of vessels: the power dissipated during Poiseuille flow, and the metaboliccost of maintaining a volume of blood. The former is minimized for wide vessels and the latter for narrow vessels, so the actualradius is a compromise.In living systems, the bifurcation exponent is often measured to deviate from three , which is not fully understood, althoughseveral factors are known to modify γ opt in single-vessel analyses. Inclusion of pulsatile flow, elastic wall vessels, and turbulence1 a r X i v : . [ q - b i o . T O ] J u l ontribute to reduction of the optimal junction exponent to γ opt = . . A key assumption leading to allometric scaling laws isthat cross-sectional area is conserved, i.e. γ opt = .Values of γ , measured in many vascular networks, are larger than expected from single-segment analyses. In some organs, γ is found to be slightly greater than three . To our knowledge, no explanation of this effect is currently available, sincecorrections to flow in single artery analyses to include turbulence, pulsatile flow, and elastic wall vessels, lead to γ opt <
3. Thissuggests a role for complexity in vascular network analysis.We propose that, in order to fully understand the optimal branching exponents in vascular trees, it is essential to take intoaccount the complexity of the entire arterial network in an organ, and the boundary conditions imposed by the organism. Asingle 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 additionalcomplexity is poorly understood. The metabolic demand of the organ determines the blood flow to the organ. The radius of theprimary artery supplying that organ is determined by a compromise between the whole organism and the organ. These twoproperties define boundary conditions for arterial growth algorithms.In this paper we carry out a theoretical and numerical analysis of the optimal bifurcation exponent for large and complexarterial trees with physiologically measured boundary conditions. The work goes beyond previous analyses , byoptimizing entire trees, rather than a single arterial bifurcation. The constraints on flow and radius of root vessels in real organsare also taken into account. Further to this analysis, we use the SALVO algorithm to determine the numerically exact globallyoptimal bifurcation exponent.
The arterial tree is divided into straight segments and bifurcations, with Poiseuille flow assumed within segments. The powercost for a single arterial tree segment experiencing Poiseuille flow is, W j = m b π r j l j + µ f j l j π r j (2)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 of blood,and µ the dynamic viscosity of blood. The power cost associated with bifurcations is neglected.The total cost, W , of an arterial tree is the sum of these individual segment costs, W = ∑ j ∈{ segments } W j . (3) Murray’s law is derived by optimizing total cost in a single segment (Eq. 2). By differentiating with respect to r j , ∂ W j ∂ r j = m b π f j l j − µ f j l j π r j . (4)When ∂ W j / ∂ r j =
0, the optimal r j can be found. This leads to a relation for flow in terms of r j f j = m / b π µ / r j (5)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 rootsegment respectively, and α is the length–radius exponent. This slightly modifies the preceding argument, so that, f j = m / b π ( µ ) / r j (cid:114) + α − α = f root r r j , (6)where f root is the flow through the root segment. a b c (a) (b) a ab c dd Figure 1.
Two types of update, that move node coordinates and swap the parent segments of nodes respectively, are requiredfor ergodicity. The figure shows a summary of these updates. In panel (a) node a is moved. In panel (b) parents of two nodes(nodes b and d) are swapped. The parent of node b is node a, and the parent of node d is node c. After the swap, the parent ofnode b is node c and the parent of node d is node a.
The power of the numerical SALVO algorithm is that the generated trees are globally optimized and therefore represent thelowest possible power cost, allowing the effects of evolutionary optimization to be investigated. While the globally optimizedsolution represents an idealized evolutionary endpoint, insight into the compromise associated with optimizing the competingcosts of the different metabolic requirements associated with maintaining a complicated vasculature will be gained by thisanalysis.The (SALVO) algorithm developed in earlier papers can be used to generate arterial trees in the 2D plane. In this section,an outline of this algorithm is given. The algorithm is similar to the approach for growing cardiac and cerebral vasculature ,with some differences relating to the use of fixed nodes to supply tissue.Equation 3 is the cost function at the core of the SALVO algorithm. In this paper we study an idealized two-dimensional(2D) piece of ‘tissue’. Leaf nodes are fixed in place, and there is no supply penalty. Similarly, there is no penalty for tissuepenetration of large vessels, since all vessels lie within the 2D tissue. The root node of the tree is fixed to the corner of a squareregion of side a . In contrast to earlier use of the algorithm, a Poisson disc process is used to place leaf nodes . The whole 2Dregion was accessible by nodes, with only metabolic- and flow-related penalties consistent with Eq. 2.On each iteration, modifications to the binary tree are attempted by either (1) moving a node or (2) changing the treestructure by swapping node connections. These updates are sufficient to ensure ergodicity. Updates are summarized in Fig. 1and Table 1. The root node is never updated. In this version of the algorithm, leaf nodes are never moved. We use simulatedannealing to optimize the cost function .Acceptance of the updates is determined according to the probability, P θ , θ + = min (cid:40) exp ( − ∆ W ( θ , θ + ) T ) , (cid:41) (7)where ∆ W ( θ , θ + ) = W ( θ + ) − W ( θ ) is the change in cost associated with modifying the tree from configuration ( θ ) toconfiguration ( θ + ) , respectively. T is the annealing temperature, which is slowly reduced using the common exponentialschedule, T θ + = ε T θ where θ is the iteration number, ε = exp ( ln T − ln T Θ ) / Θ , Θ the total number of iterations, and T ( T Θ )are the initial (final) temperatures.SALVO was implemented in C++ making full use of the 2011 standard library (g++ version 7.4.0 compiled with the -O3flag). Generated trees were analyzed using Python 3.7. A Threadripper 2990WX processor was used for the calculations, withcalculations for different γ , Ω and tree sizes sent to different threads by a Python script. The optimization of a 500 node treetakes approximately 30 minutes, and a 5000 node tree takes approximately 11 hours on a single thread for 10 updates. Arteries can be grouped together, so that each group comprises arteries with identical properties (e.g length, diameter, flow). Ina real arterial system, this would not be true, but it would still be possible to group arteries with similar lengths, radii, and flowstogether. ithout loss of generality, the total power can be rewritten as, W = ∑ j ∈{ r , l , q } N ( r j , l j , f j ) (cid:32) m b π r j l j + µ f j l j π r j (cid:33) , (8)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, f n = n f leaf (9)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: f n = f leaf ( r n / r leaf ) γ (10)thus, r n = r leaf ( f n / f leaf ) / γ = r leaf n / γ . (11)Experimental data suggest that the length of an arterial segment is proportional to a power of the radius, l n = l leaf ( r n / r leaf ) α = l leaf ( f n / f leaf ) α / γ , (12)where the value of the exponent α is typically close to 1.0 .By substituting Eqs. 12, the power required to maintain blood flow through a segment depends only on the flow f n , W n = W ( f n ) = m b π r l leaf ( f n / f leaf ) ( + α ) / γ + µ l leaf π r f n ( f n / f leaf ) ( α − ) / γ . (13)Thus, the dimensionless metabolic ratio , defined as Ω = m b π r / µ f , along with N , controls location in parameterspace. W n = C (cid:16) Ω n ( + α ) / γ + n +( α − ) / γ (cid:17) (14) = Cn +( α − ) / γ (cid:16) Ω n / γ − + n − / γ (cid:17) (15) C = µ f l leaf π r . Both C and Ω are defined in terms of the leaf node.A similar ratio for the root node, Ω root = m b π r / µ f can be defined for convenient contact with experiment. Thevalues r root and f root are often known from experiment, e.g. Doppler ultrasound, and N can be estimated. This ratio can berelated to Ω via Ω root = N / γ − Ω .The total power requirement is, W = ∑ n N n W n . (16) N n is the number of segments with flow n f leaf , and simplifies the function N ( r i , l i , f i ) . For any tree structure, N is always thenumber of leaf nodes, so N = N . There is always a single root node with total flow N f leaf , so N N =
1. The remaining N n aredependent on the structure of the tree. At each bifurcation, flow conservation requires that n in = n out , + n out , , so n is an integer.Total power is linear in length scale, so the structure of the power landscape (including the location of any minima withrespect to γ ) is independent of a . It is the global minimum with respect to γ that sets the structure of the tree, and when locatingthe minimum, ∂ W / ∂ γ =
0, so the factors of l leaf simply cancel, 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 .There are two special cases: fully symmetric and fully asymmetric. In the first case, identified as a fully symmetric tree, theflow is split evenly at each bifurcation. For the case which we shall identify as fully asymmetric, a single leaf node emerges ateach bifurcation and the rest of the flow passes down the other bifurcation. We will explore these special cases in the followingtwo sections. (a) γ op t Ω asymmetric, N =2.047x10 symmetric, N =1.023x10 symmetric, N =2.047x10 symmetric, N =1.049x10 symmetric, N =1.074x10 SALVO, N =2.160x10 (b) γ op t Ω α =0.89 α =1.00 α =1.05 α =1.15 Figure 2. (a) Deviations from Murray’s law ( γ opt =
3) depend strongly on changes in the metabolic ratio, but are insensitive tothe structure of the tree. The figure shows a comparison of γ opt vs Ω for fully symmetric, asymmetric and numerical trees. (b)Optimal bifurcation exponent is insensitive to changes in α . In a fully symmetric tree, all of the segments with flow n exist at the same bifurcation layer. Each layer, m , has 2 m segments,where m is the number of bifurcations upstream of that layer ( m = M layers.The total power cost can be determined by substituting the definitions N n = m if n = M − m otherwise N n =
0, into Eq. 16, W = C M ∑ m = m (cid:16) Ω ( M − m )( + α ) / γ + ( M − m )( +( α − ) / γ ) (cid:17) . (17)Thus, by summing the geometric series, the total power cost for a fully symmetric tree is, W = M C (cid:32) Ω − M ( − ( + α ) / γ ) − − − ( + α ) / γ + M ( +( α − ) / γ ) − − − ( +( α − ) / γ ) (cid:33) (18) 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, except there are N leaf nodes so N = N .Substitution into Eq. 16 gives, W = C (cid:32) N ∑ n = ( Ω n ( + α ) / γ + n +( α − ) / γ ) + ( Ω + )( N − ) (cid:33) (19)So the total power cost for an asymmetric tree is W = C (cid:16) Ω H ( − ( + α ) / γ ) N + H ( − ( +( α − ) / γ )) N (cid:17) + C ( Ω + )( N − ) , (20)where H ( r ) n is the generalized harmonic function, ∑ nk = / k r . The optimal value of γ is obtained by solving ∂ W / ∂ γ = ∂ W / ∂ γ = Ω , which can change due to physiologicalboundary conditions on flow and radius at the input vessels. These constraints may be due to limits in the size of the largestvessel and changing flow demands of tissue. Figure 2(a) shows the optimal value of γ . When Ω = α = γ opt =
3) is recovered. γ op t N Ω =0.1 Ω =0.5 Ω =0.9 Ω =1 Ω =1.1 Ω =1.5 Ω =2 Ω =5 Ω =10 Figure 3.
Deviations from Murray’s law are largest for small trees and strongly dependent on changes in the metabolic ratio.The figure shows γ opt vs N for a fully asymmetric tree. Table 1.
Simulation parameters and their ranges.Name symbol rangebifurcation exponent γ Ω N µ . × − Pa stissue size a Θ (10 for checks)SA initial ‘temperature’ T − SA final ‘temperature’ T Θ − Js − short move distance d move . d move . γ opt is qualitatively unchanged by the structure of the tree. Results for asymmetric and symmetric trees with N = . × follow essentially the same functional forms. The optimal bifurcation exponent for the asymmetric tree is closer to γ = γ = α (Fig. 2(b)). Thisstructural effect potentially has implications for the value of γ opt in organs, since α can vary with organ type, with estimatesranging from 0 . − .
15. In practice, changes in γ opt for this variation in α are far less than the error for measurements of γ and changes in α can essentially be neglected.Deviations from Murray’s law are larger for smaller trees and strongly dependent on changes in the metabolic ratio. Thelarger the tree, the closer to Murray’s law γ opt becomes. Figure 3 shows variation of γ opt with N for fully symmetric trees. Forvascular tree sizes of between 10 and 10 segments, which are typical in organs, γ opt ranges between 2 and 4. The generation of globally optimal trees using a numerical algorithm helps to test analytic expressions, and provides additionalmorphological measures that can be used to understand arterial networks. In this section, we use SALVO to investigate the roleof vascular complexity and physiological boundary conditions on the properties of globally optimal trees. Several properties ofthe numerically generated trees are investigated. We determine the sensitivity of globally optimal trees to Ω and γ . Throughexamination of W tot , we compute γ opt for complex trees. For each value of γ and Ω investigated, arterial trees with up to 5000nodes were generated. Table 1 summarizes the parameters used for the numerical calculations.Three sectors of the parameter space have qualitatively different tree structures (Fig. 4): (1) for γ (cid:46) , Ω > γ (cid:38) , Ω < γ (cid:46) , Ω < γ (cid:38) , Ω > (cid:46) γ (cid:46) Figure 4.
The structure of the globally optimal vasculature varies with γ and Ω . Trees have size N = Ω (cid:29) γ < Ω (cid:28) γ > α ≈
1, the power in a segment is W n = Cn ( Ω n / γ − + n − / γ ) . For γ >
3, the exponents (whichinvolve 3 / γ −
1) have opposite sign to those for γ <
3. So after the substitutions Ω = / Ω (cid:48) , γ = γ (cid:48) / ( γ (cid:48) − ) , C (cid:48) = Ω C , W n = C (cid:48) n ( Ω (cid:48) n / γ (cid:48) − + n − / γ (cid:48) ) , and the sum has an equivalent structure. The substitution is determined by identifying where1 − / γ = / γ (cid:48) −
1. Since the prefactor C (cid:48) scales the entire sum, then the minima of W and thus the results for γ , Ω and γ (cid:48) , Ω (cid:48) are identical. This symmetry is only approximate if α (cid:54) = γ (cid:46) Ω (and γ (cid:38) Ω ), the tree structure is highly asymmetric, with long trunks snaking throughleaf node sites (see top left panels in Fig. 4). This is due to the domination of the n − / γ term due to Poiseuille flow for low γ ,and the n − / γ (cid:48) metabolic cost term for large γ . Thus terms with large n (i.e. thick trunks) are favored.For γ (cid:46) Ω (and γ (cid:38) Ω ), long leaf segments connect root and leaf nodes (see top right panels in Fig.4). This is due to the domination of the n / γ − term due to metabolic maintenance of blood for low γ , and the n / γ (cid:48) − Poiseuilleterm for large γ . The term related to metabolic maintenance of blood (with n / γ − ) dominates. Thus, terms with small n (i.e.leaf nodes) are favored.For the biologically relevant regime, 2 < γ <
4, trees have a symmetric structure. No single term in W dominates. There issurprisingly little variation between the tree structures in this region.To quantify the effect of varying γ and Ω on the network structure, we have examined average segment length, path lengthand radius asymmetry. The radius of an arterial segment is given by r j , and the length by l j . Average length is defined as =2163 l/ a L / a r c > / ( r c < + r c > ) γ N =3968 γ Ω =0.1 Ω =0.5 Ω =1.1 Ω =2 Ω =5 Ω =101 2 3 4 5 Figure 5.
The tree morphology is highly sensitive to variation in the bifurcation exponent, but relatively insensitive tovariation in Ω within the region of interest between γ = γ =
4. There is essentially no sensitivity to tree size. l = ∑ l j / N . The average summed path length from root to leaf node is L = (cid:104) ∑ path l j (cid:105) . Radius asymmetry is measured using (cid:104) r c > / ( r c < + r c > ) (cid:105) (where r c > ≥ r c < ). The sensitivities of these quantities to variations in γ and Ω are presented in Figure 5.In the typical range of biological tissue (2 < γ < γ . Allmorphological properties are insensitive to variation in Ω . Average segment length is short and path length is long in this region,consistent with the branching structures seen for intermediate γ in Fig. 4. Bifurcation symmetry is in the range 0 . − .
62, sobifurcations are moderately symmetric. Although Ω leads to minor changes in tree morphology in this regime, we note it canaffect γ opt and thus the tree morphology via γ as a secondary effect.In the regions γ < γ > Ω is responsible for huge variations in the tree morphology, and γ can also producelarge variations in the various morphological and structural properties of the tree. Path length drops outside this region toapproximately a / √ γ < , Ω (cid:28)
1, theasymmetry increases dramatically. For all other regions of the parameter space, the asymmetry drops.Morphological measurements are essentially insensitive to changes in N , consistent with additional segments adding moredetail, but not qualitatively changing the tree structure. Panels on the left of Fig. 5 show results for N = N = γ opt can be determined without ambiguity from the minimum in W . Figure 6(a) shows howthe total power cost varies with γ . There is a clearly defined global minimum for all values of Ω shown. γ opt can be found byfitting 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 fromanalytic expressions. Numerical values of γ opt for various values of Ω vs N are shown in Fig. 6(b), and compare favorably toFig. 3. Several numerical values are compared with the analytic results in Fig. 2(a), also showing good agreement for bothsymmetric and asymmetric trees.A power-law relationship, l = Ar α , is found for the median segment length in terms of segment radius calculated using / o p t D 1 e x t r m o p t E Figure 6. (a) A well defined global minimum in total power cost means that optimal bifurcation exponent γ opt can bedetermined without ambiguity. The figure shows total power cost as a function of bifurcation exponent γ for several values of Ω . (b) The relationship of γ opt to N and Ω , numerically determined using SALVO, is qualitatively similar to the relationshipdetermined from analytic expressions. The figure shows γ opt vs N for several Ω .SALVO (Fig. 7(a)). The variation of the value of l / r root with r / r root is analyzed using Python 3. The expression l = Ar α isfitted to the median value. Figure 7(a) shows the fit. Regions shaded light blue show the range between the 25th and 75thpercentiles in the length histograms. For segments selected from trees with N > . < γ < . Ω = .
9, the fit hasexponent α = . ± . . As in experiment, there is a strong scatter on the length.The length–radius exponent, α , is consistent with experimental values for trees grown with realistic 2 . < γ < .
5, but canbecome effectively negative when long leaf segments dominate outside this region (Fig. 7(b)). To calculate the length–radiusrelation, segments are binned from trees with N > γ and Ω values. Where exponents are negative, the relationonly poorly follows a power law, and errors on α are large. The relation is well followed within the region 2 . < γ < .
5, andthis leads to smaller error bars. Overall, errors on α are relatively large, and could be reduced by making calculations foradditional N . In this paper we determined analytic expressions, and carried out numerical calculations, for the properties and structures ofglobally optimal vascular trees, with the aim of understanding how overall complexity and physiological boundary conditionscontribute to the optimal junction exponent and other structural properties of arterial trees. Analytic expressions were derivedfor the special cases of maximally symmetric and asymmetric arterial trees. The parameter space of the arterial trees was morefully explored by making numerical calculations with SALVO, enabling globally optimal vasculatures to be found for arbitrarytree morphology. Tree structures, morphological properties, sensitivity to dimensionless parameters and optimal bifurcation(junction) exponent are calculated.Our analytic expressions are consistent with numerical calculations, and predict that γ opt is insensitive to tree symmetry, sowe propose that the analytic expressions derived here are applicable to a wide range of vasculatures. Analytic expressions canbe used for much larger trees, and would, therefore, be useful for predicting the properties of vasculatures within a range oforgans where the number of vessel segments and overall complexity exceed the capabilities of current computers. We expectthat it will be possible to extend the analytic expressions to include pulsatile flow and turbulence, and will investigate thispossibility in future studies.We predict that tree complexity is a significant contributor to the bifurcation exponents in living organisms. The deviationswe find from 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 Ω (cid:54) =
1. This may occur since allorgans, with their dramatically varying demands, are connected to the same major vasculature. We expect that large variationsof γ opt with increasing complexity will also occur if a more detailed analysis including pulsatile flow and turbulence is carriedout.We predict that arterial tree complexity can lead to optimal bifurcation exponent, γ opt >
3, a situation which can be found inexperiment, and is of interest since inclusion of turbulence and pulsatile flow in single artery analyses leads to γ opt <
3. Large r / r root l / r r oo t D P H G L D Q Ax = 0.887 ± 0.088 -1-0.5 0 0.5 1 1.5 2 2 2.5 3 3.5 4 (b) α γΩ =0.1 Ω =0.5 Ω =2.0 Ω =10.0 Figure 7. (a) A power law relationship is found for the median segment length in terms of segment radius calculated usingSALVO. The figure shows median values of l / r root vs r / r root , a power law fit (dashed line), and the 25th and 75th percentiles(light blue shading). To calculate the length–radius relation, segments are binned from trees with N > . < γ < . Ω = .
9. (b) The length–radius exponent, α , is close to one for trees grown with 2 . < γ < .
5. To calculate the length–radiusrelation, segments are binned from trees with N > γ and Ω values.values of γ are measured in e.g. the brain vasculature ( γ = . , retina ( γ = . , γ = . ± . ) and other mammalianvasculatures where γ can range as high as 4 . Such large γ are not predicted by single segment analyses including effects relatedto pulsatile flow, elastic vessel walls and turbulence ( γ = . . Complexity and boundary conditions provide an additionalcontribution that can account for larger values of γ opt .We predict that tree structures within the physiological regime are only sensitive to γ ; outside the physiological regimestructures are also highly sensitive to Ω ; and for all regimes tree structures are insensitive to N . Changes in N do not qualitativelychange the morphology of the tree, but add more detail. Outside the regime 2 < γ <
4, structure can change dramatically with Ω . 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 setthe radii in the generated trees . Similarly, allometric scaling arguments require knowledge of γ , and variations of γ opt could modify such approaches. γ opt is quite hard to measure experimentally, leading to values with large uncertainties, and weconsider the calculation of such values to be a useful application of our technique.For values of γ consistent with living systems, we find power law exponents in our computational trees that are consistentwith the value α ∼ . < α < . . We find a similar rangeof values in our numerical calculations, and with improved description of the flow, the accuracy of the predictions could beimproved. Values of α are also useful as input to other calculations.Future work to include additional physics, such as pulsatile flow, turbulence and vessel elasticity, would lead to acomputational model with enhanced predictive power. These improvements to the treatment of flow through vessels couldbe incorporated into both the analytic expressions derived in this paper, and into the cost function of SALVO without havingto change the core algorithm. Once analytical expressions are modified to include this additional physics, we suggest thatparameters such as m b could be determined from empirical results.The significant structural changes visible at γ ∼ γ ∼ Contributions
JPH carried out the analytical calculations and led the study. JK carried out the numerical calculations. Both authors contributedto writing of the manuscript and analysis of the data. ata availability
The datasets generated and analyzed during the current study are available in the ORDO repository:doi.org/10.21954/ou.rd.12220490 (note, these will be added at proof stage, data available on request).
Acknowledgments
The authors have no competing interests. JK would like to acknowledge EPSRC grant no. EP/P505046/1.
References Nakamura, Y. & Awa, S. Radius exponent in elastic and rigid arterial models optimized by the least energy principle.
Physiol. Rep. , e00236, DOI: 10.1002/phy2.236 (2014). Keelan, J., Chung, E. M. L. & Hague, J. P. Simulated annealing approach to vascular structure with application to thecoronary arteries.
R. Soc. Open. Sci. , 150431 (2016). Keelan, J., Chung, E. M. L. & Hague, J. P. Development of a globally optimised model of the cerebral arteries.
Phys. Med.Biol. , 125021 (2019). Changizi, M. A. & Cherniak, C. Modeling the large-scale geometry of human coronary arteries.
Can. J. Physiol. Pharma. , 603–611 (2000). Murray, C. D. The physiological principle of minimum work: I. the vascular system and the cost of blood volume.
Proc.Natl. Acad. Sci. , 207–14 (1926). West, G. B., Brown, J. H. & Enquist, B. J. A general model for the origin of allometric scaling laws in biology.
Science , 122–126 (1997). Zamir, M., Medeiros, J. A. & Cunningham, T. K. Arterial bifurcations in the human retina.
J. Gen. Physiol. , 537–548(1979). Wischgoll, T., Choy, J. S. & Kassab, G. S. Extraction of morphometry and branching angles of porcine coronary arterialtree from CT images.
Am. J. Physiol. Hear. Circ. Physiol. , H1949–55, DOI: 10.1152/ajpheart.00093.2009 (2009). Zamir, M. Fractal dimensions and multifractility in vascular branching.
J. Theor. Biol. , 183–90, DOI: 10.1006/jtbi.2001.2367 (2001).
Horsfield, K. & Woldenberg, M. J. Diameters and cross-sectional areas of branches in the human pulmonary arterial tree.
Anat. Rec. , 245–251 (1989).
Bridson, R. Fast poisson disk sampling in arbitrary dimensions.
Proc. ACM SIGGRAPH (2007).
Aarts, E., Korst, J. & Michiels, W. Simulated annealing. In Burke, E. K. & Kendall, G. (eds.)
Search methodologies ,187–210 (Berlin: Springer, 2005).
Kamiya, A. & Takahashi, T. Quantitative assessments of morphological and functional properties of biological trees basedon their fractal nature.
J. Appl. Physiol. , 2315–2323, DOI: 10.1152/japplphysiol.00856.2006 (2007).
Habib, M., Al-Diri, B., James, L., Hunter, A. & Steel, D. Constancy of retinal vascular bifurcation geometry across thenormal fundus and between venous to arterial bifurcations.
Invest. Opth. Vis. Sci. (2006). Al-Diri, B., Hunter, A., Steel, D. & Habib, M. Manual measurement of retinal bifurcation features.
Conf. Proc. IEEE Eng.Med. Biol. Soc.
Schreiner, W. & Buxbaum, P. F. Computer-optimization of vascular trees.
IEEE Trans. Biomed. Eng. , 482–491, DOI:10.1109/10.243413 (1993). Schreiner, W. et al.
Optimized arterial trees supplying hollow organs.
Med. Eng. Phys. , 416–429 (2006). Karch, R., Neumann, F., Neumann, M. & Schreiner, W. Staged growth of optimized arterial model trees.
Ann. Biomed.Eng. , 495–511 (2000)., 495–511 (2000).