Simulated annealing approach to vascular structure with application to the coronary arteries.

Do the complex processes of angiogenesis during organism development ultimately lead to a near optimal coronary vasculature in the organs of adult mammals? We examine this hypothesis using a powerful and universal method, built on physical and physiological principles, for the determination of globally energetically optimal arterial trees. The method is based on simulated annealing, and can be used to examine arteries in hollow organs with arbitrary tissue geometries. We demonstrate that the approach can generate in silico vasculatures which closely match porcine anatomical data for the coronary arteries on all length scales, and that the optimized arterial trees improve systematically as computational time increases. The method presented here is general, and could in principle be used to examine the arteries of other organs. Potential applications include improvement of medical imaging analysis and the design of vascular trees for artificial organs.

Introduction

Arterial trees are vital for the efficient transport of oxygen and nutrients to tissue. Their anatomy has been studied for many centuries through the dissection of cadavers, inspection of corrosion casts, medical imaging techniques and computational models. It has been determined that individual arterial bifurcations follow optimality principles that lower metabolic demand locally [1-6], as demonstrated by the scaling laws followed by arterial trees [7-11]. More recently, there has been a high level of interest in models that mimic arterial growth (angiogenesis) using physical and physiological principles to simulate vascular anatomy. These models are created based on local optimization principles, where the anatomy of each branch in the arterial tree is governed by a compromise between maximizing fluid dynamical efficiency and minimizing the quantity of blood required. However, models of coronary vasculature, based on local optimization are not able to explain if the organization of major arteries is the result of fluid dynamical optimization across the ‘whole organ’ [12-16].

The relationship between the radii of vessels in individual bifurcations is well categorized by the equation: , where r is the radius of the parent artery, r those of the daughter arteries and γ is the bifurcation exponent [4,17,18], which has the value 3.0 in Murray's formulation. This relation arises from the combination of the continuity equation and the diameter–flow rate relation [9]. Several current methods for in silico growth of vascular trees into a simulated tissue substrate aim to optimize the local properties of individual bifurcations [15,19]. A procedure, known as constrained constructive optimization (CCO), starts by inserting a single artery into the tissue. A new vessel with a position chosen at random is then connected to the original artery and the link point is moved such that the energy of the arteries is minimized. New arteries are then iteratively added and optimized until a predetermined number of terminal sites have been added. The overall result is that CCO and similar methods create trees whose structure is predetermined by the order in which new arteries are added: if the order is changed, the final tree structure also changes. Morphologically, CCO reproduces a reasonable distribution of vessel sizes due to the application of Murray's law, but creates arterial branches that are more symmetric than those found in nature (especially for the largest arteries) [20] and significant extensions are required to generate vessels in hollow organs [21]. The variations in the structure and positions of larger arteries in CCO generated trees are problematic, as organs such as the heart exhibit only small differences in large artery structure over a population (aside from rare abnormalities). An alternative method known as global constructive optimization (GCO) attempts to overcome the problems of CCO by including a multiscale pruning update that is global in the sense that it acts simultaneously on a significant subset of the tree, but otherwise only includes updates allowing local modifications to the topology of the tree. As such, GCO is limited to sampling a subset of the allowed topologies of the arterial trees [22], so while it is expected to offer improvements over CCO, it carries no guarantee of reaching the global minimum. Due to the use of local downhill searches, GCO also has similar issues with hollow organs. To obtain a universal optimization technique to compute in silico arterial trees for arbitrary tissue structures, a different approach is needed.

Another method for the generation of large scale arterial trees uses extensive morphological databases [14,23]. These trees contain far more vessels than is feasible to generate with techniques such as CCO, as the topology of the tree is taken from experimental data. However, as detailed morphological databases do not exist for the vast majority of organs, the use of these techniques is impossible in the general case. Morphologically generated models provide trees suitable for large scale fluid dynamical studies and organ phantoms [24]. They achieve this by reproducing experimental data in a computationally accessible form. As such they have no predictive powers that can contribute to the understanding of the origins of arterial tree structure.

A separate class of models exist for use in modelling the growth of tumourous vasculature and the process of vascular remodelling, which involve the direct simulation of sprouting angiogenesis [25-28]. These models seek to reproduce the rapid and dynamic process of tumour vascularization, or the growth of vasculature in normal tissue as it grows. Using in silico simulation of sprouting angiogenesis to obtain the vascular structure in a fully grown organ would be extremely difficult, as full details of the distribution of tissue and oxygen demands would be needed for all stages of embryonic and childhood development. As organs such as the heart have very small levels of variation in vascular structure over the population, the structure itself is likely to be caused by a process different from that of tumour vascularization. We suggest that this process is an optimum seeking one, and that as such an optimization procedure is required to accurately model it. We examine if, regardless of the complex processes that guide angiogenesis during growth, the final structure of the coronary vasculature in adult mammals is near optimal.

Development of a method which reaches a morphologically accurate solution based solely upon optimization criteria would be useful in vascular research, allowing for the modelling of realistic vascular trees in organs lacking extensive morphological databases. The inability of CCO and extensions to find the global energy minimum, and the subsequent lack of consistent structure (particularly of the larger arteries), is problematic if organ specific vasculature is required. An approach capable of producing an arterial tree, which minimizes pumping power and blood volume, while providing adequate blood flow to critical regions would be invaluable in this regard. This paper goes beyond previous work by introducing a far more flexible and universal method for generation of ‘whole organ’ arterial trees, in any arbitrarily shaped tissue substrate, that obeys both local and global optimization criteria. To identify globally optimized arterial trees, we use a powerful computational technique known as simulated annealing (SA) [29]. Although SA is computationally expensive, correctly applied SA techniques have a key advantage of being mathematically and computationally proved to converge to a global energy minimum. To achieve this, our SA-based approach has the potential to sample all possible arterial tree configurations, ranging from perfectly symmetric, intricately bifurcating structures, to asymmetric trees characterized by a single trunk vessel. This is achieved by allowing: (i) repositioning of bifurcations and (ii) swapping the parent vessels of bifurcations between different parts of the tree. By introducing these forms of plasticity to our models, the entire parameter space of the tree can be explored, allowing the method to identify the best possible arterial configuration for supplying a particular organ. Full details of this novel method can be found at the end of the article. As an example application we determine the near optimal configuration of arteries for supplying the heart and compare our computer generated coronary vasculature with morphological data from real coronary arteries. Specifically, we determine that the observed anatomy of the coronary arteries is similar to that expected from near global minimization of total energy expenditure, and validate the approach against porcine data, finding a very high level of agreement with morphological data.

Methods

The main purpose of any arterial tree is to maintain adequate blood perfusion with minimal total metabolic expense. The suitability of an arterial tree for this purpose is governed by two considerations: (i) as blood is viscous, the power required to pump blood through the vasculature should be minimized and (ii) as energy is required to generate and maintain blood, the volume of blood required should be minimized. Murray's law achieves this for individual bifurcations, but the optimal organization of large numbers of connected bifurcations is far from obvious. The interplay between these competing concerns for thousands of arterial segments leads to a complex optimization problem. Note that in the following, bifurcations will be referred to as nodes, arteries will be referred to as ‘segments between nodes’, and terminal arterioles are referred to as ‘end nodes’.

Metabolic cost to maintain blood volume

The first component of the approach involves calculating the power needed to maintain the entire tree, which will be used as a value in the cost function. The power consumption of the tree can be split into two separate parts: the first is the metabolic cost of maintaining the blood volume and tissue associated with the tree, and the second is the power required to pump blood through the tree. The length and radius of each segment (vessel) i of the tree must be known to calculate the volume. By assuming a fixed bifurcation exponent, the radii are determined by the topology and only vessel lengths rely on the geometrical arrangement. To calculate the cost, volume must be multiplied by a constant, m, corresponding to a physiologically reasonable metabolic demand of the same quantity of blood and vascular tissue [30]. Thus, the metabolic cost due to the volume of the tree will be given by where m is taken to be 641.3 J s−1 mm−3 and Vtree is the volume of the entire tree.

Power cost to pump blood through vessels

To calculate the power needed to pump blood through the entire tree, we must know the pressure and volumetric flows inside each segment (vessel) of the tree, which can be found by first assuming that Poiseuille's law, Δp=QR, is followed inside the segments, where Δp is the pressure drop over the vessel, and Q is the flow. The assumption that flow is laminar inside the vessels is justified provided that the typical length of a vessel is much larger than the radius, and that pulsatile flow effects are negligible. Vessels within the simulated trees have a typical length radius ratio of 10, and while in the largest arteries of the tree pulsatile effects may still be present, these rapidly decay so that the vast majority lie within a non-pulsatile regime. We assume both Murray's law and that terminal node flows are constant to simplify calculation of the relevant fluid dynamical quantities: the only quantity which relies on the structure of the tree is the pressure. In a sense, the segments can be considered as connected set of resistors, with the resistance given by where r is the radius of the vessel, L its length and μ=3.6×10−3 Pa s the viscosity of blood. The pressures (and hence flows) for every node in the tree can then be found recursively. W, the power consumed by each segment i, is then calculated using Summing over all segments in the tree, the total power required to maintain the proper flow through the tree is

Ensuring tissue supply

The primary purpose of the vascular tree is to supply blood; thus it is important that terminal nodes are correctly dispersed inside the tissue. Initially, terminal nodes are randomly distributed inside the tissue, with each node having associated with it a sphere of influence for blood supply. The radius of this sphere is calculated using physiological values for the blood demand of the tissue. The density of myocardium is ρ=1.06×103 kg m−3 [31], and the flow demand is  [32] leading to a flow demand per cubic metre of heart tissue of qrequired=2×10−2 s−1. The total flow into the heart is Q0=4.16×10−6 m3 s−1 [33], which can be converted to total flow per node as Q=Q0/N, where N is the total number of arterioles (end nodes). The radius of the supply sphere is then calculated via 4πR3supply/3=Q/qrequired. The sphere can be thought of as a microcirculatory ‘black box’ [15], where the exact fluid dynamical details of the blood flow have been ignored. Spheres of blood supply associated with end nodes are stored in a voxel map (a voxel is a three-dimensional generalization of a pixel) of the tissue, where each terminal node adds exactly one to each voxel inside its sphere of supply (figure 1b). While blood demand is not constant within the myocardium at any single instance of time, the majority of fluctuations are high frequency oscillations which are assumed to be averaged out in the present model [34]. The terminal nodes are then allowed to move inside the tissue, where after each move a new voxel supply map is calculated, and the overlap (each voxel supplied by more than one sphere, or the dark red voxels in figure 1b) is used as a value in the cost function of the SA algorithm. In addition, all voxels not being supplied are given a cost, so that the overall penalty associated with having both unsupplied and oversupplied voxels is chosen to be where b is the value of the supply at the voxel and the sum is performed over all the voxels comprising the tissue. In practice, this cost is set to be much larger than all other costs, as any unsupplied tissue would die. Therefore, the terminal nodes spread evenly through the tissue early in the optimization. Other functions may be used, provided that the minimum in the function for each voxel occurs at b=1. The value Cs then defines the fitness of the tree to supply blood, and the penalty for over supplying voxels forms a sort of self-avoidance algorithm, where terminal nodes are encouraged to pack the tissue as densely as possible without overlapping. A benefit of this method is that it allows easy integration of medical imaging into the model, as well as providing an easy method for differentiating tissue with different blood supply demands.

Figure 1.

(a) The distance map for a spherical surface. Voxels outside the surface have value 0, with those inside the surface contributing a value relating to their distance from the surface. (b) Each arteriole supplies a spherical region shown by the lightly shaded squares. Where there is significant overlap between two spheres, there is a penalty. Unsupplied voxels also incur a penalty in the cost function.

Exclusion of large vessels from tissue

In order to create a realistic vascular tree, it must be possible to exclude some segments from penetrating the tissue. For instance, in the case of the heart, it would be unlikely to find a very large artery within the myocardium, and vessels may not penetrate the ventricles; rather, the larger arteries and arterioles lie on the surface of the heart, with only the smaller arterioles and capillaries being found inside the tissue. To mimic this structure, the approach makes use of a cut-off radius Rcutoff, whereby segments with radius larger than Rcutoff may not penetrate the tissue. In the calculations performed in this article, Rcutoff=0.01 mm. To determine which segments with radius greater than Rcutoff have penetrated the tissue we first take a distance transform of the tissue surface for each tissue voxel (figure 1a). This provides a second voxel map of the tissue, distinct from the blood supply map, giving a measure of the distance of a point from the surface when it is inside the tissue (outside of the surface, the value is zero). For each segment satisfying the radius criteria, a list of voxels that its centre-line penetrates is generated [35], along with a value for the length element of the segment present inside that voxel. A cost is then calculated based upon the value of the distance transform at each of the voxels according to where i, j and k are the Cartesian voxel coordinates taken from the centre-line of the segment. D is the value of the distance transform at that voxel coordinate. is the length of the segment spent inside the voxel. The sum is performed over all the voxels contained in the list calculated from the centre-line. This cost can then be used in the SA algorithm as a penalty that favours moving large segments out of the tissue.

Pressure constraints

In physiologically realistic trees, capillary networks should receive a constant pressure Pterm to function correctly. A new cost can be devised to ensure this. A suitable candidate is where the sum is performed over all terminal nodes, and P is the actual terminal node pressure. In practice, for trees which can be optimized on feasible time scales (i.e. of a few thousand nodes), the pressure drop from root to end node is less than 1% of the total pressure drop of a real arterial tree, with most of the pressure drop occurring over smaller arterioles than those considered here, so it is unnecessary to perform this calculation. When it becomes possible to grow larger trees, the pressure at the capillaries will need to be taken into consideration. This will add a significant computational cost.

Total cost function

We have now determined a form for all the relevant costs associated with an arbitrary tree configuration supplying arbitrary tissue shapes. We can therefore define a total cost which gives a numeric measure of the fitness of a given tree, where A indicates a weighting value which scales each relevant cost. There is no way to analytically determine what weights to use, and the selection of appropriate weights must be found experimentally; however, a few basic principles such as having a very high weight for the blood supply cost and a low weight for the end node pressure cost can guide the process. In principle, As should be infinite, as tissue without supply dies. In this work, we use A=1×104, A=0, As=1×1030 and A=100. In this way, As and A force the exclusion of vessels and uniform supply of tissue to act like constraints. While the exclusion cost should technically be infinite, as no arteries are found in the ventricles of living humans, it is advantageous to give it a large but finite value. This allows the optimization procedure to identify gradients, giving it extra information and speeding up convergence. This size of the constant for A may seem small in this regard; however, its value C is already raised to the power 6 in equation (2.6). As any scaling of the cost function does not effect the location of minima, we can absorb one of the weightings by scaling everything else. This allows a new cost function to be defined.

Simulated annealing

To select the fittest, most optimized trees, we use a powerful technique for optimization problems known as SA [36,37]. The primary difference between SA and a conventional downhill search is that SA also spends some time exploring solutions with higher cost function, and in this way can climb out of shallow valleys in the fitness function to explore other deeper regions.

The total cost C will play the role of energy in the SA algorithm, so that the probability of accepting a change to a tree of cost C, resulting in a tree of cost C is given as where is the probability of going from state i to state f, is the change in the cost function associated with going from state i to f, and T is the SA temperature parameter (not be confused with ambient temperature). The small probability to accept a higher cost tree during update allows the tree to climb out of local valleys in the cost function. The algorithm proceeds by making changes to the tree structure, calculating the change in cost function, and then either accepting or rejecting the change by comparing with a random number between 0 and 1. T starts large and is reduced slowly. If T has been reduced sufficiently slowly, then the global minimum of the cost function is guaranteed to be reached. In practice, the problem space is too large to achieve this in reasonable time, and slightly different trees with very similar cost are found if the algorithm is run with several random number seeds. The most important consideration is the lowest achieved cost. As such, if the structure of trees generated varies between different runs, we always display data from the run with the lowest cost function. As computational power increases, longer runs will be achievable leading to progressively better optimizations. The highest T used here is 1×1010, dropping during the algorithm to 10−5. Typically, a tree containing 1000 nodes will need 109 updates, with a doubling of nodes taking roughly quadruple the number of updates (up to approx. 6000 nodes with one month of CPU time). The large value of As means that the supply of tissue is determined by downhill search, while all other costs are minimized by SA.

Exploring the tree structure: translations and node swaps

The SA algorithm must have access to set of updates which allow it to alter the configuration of the tree. It is necessary to find changes that can be made to the topological and geometrical structure of the tree such that all possible solutions between perfectly symmetric structures and a single trunk vessel can be explored (i.e. the algorithm is ergodic). This is achieved by allowing: (i) repositioning of bifurcations, which is achieved by translating a node in space (figure 2a) and (ii) swapping the parent vessels of bifurcations between different parts of the tree (figure 2b). For all nodes but the root node, this move is valid, and performed consecutively it allows all possible tree topologies to be explored. If one of the two nodes is a direct parent of the other (i.e. while traversing up the tree from one of the chosen nodes, the other node is encountered) the move is rejected to avoid forming a closed loop. With these two updates, the entire parameter space of the tree can be explored, allowing the algorithm the opportunity to reach a globally optimal solution.

Figure 2.

Tree modification updates.

Strahler order

The Strahler (or stream) ordering method was first introduced to classify river systems, but can be applied to any bifurcating system. In standard Strahler ordering, nodes at the end of a tree (in this case the arterioles) are assigned a number 1. At a bifurcation, if two vessels (segments) of the same order meet, then the order of the parent vessel is 1 higher. However, if two vessels of different orders meet, the artery supplying these vessels has the largest order of the two. For example, if two arteries of order 1 meet, then the vessel supplying these arteries has order 2. If an artery of order 3 meets an artery of order 2, then the vessel supplying these arteries has order 3 (an example is shown in figure 3). Therefore, within this scheme, vessels with the lowest order are arterioles. The major vessels have the largest order. The Strahler order used here is then diameter adjusted following the approach in [38].

Figure 3.

Schematic of the Strahler ordering process.

Within the Strahler ordering scheme it is possible to identify continuous sections of vessels with the same order number. These are referred to as elements, so a single arterial element may pass through multiple bifurcations. Throughout this article, it is the properties of elements which will be calculated for direct comparison with [9]. We note that due to the early termination of the simulated trees, calculated order numbers are modified so that the root nodes have an order number equivalent to that of the largest arteries of real coronary arterial trees. For example, in the work of Kassab, the largest diameter defined Strahler order number is 11, corresponding to the input artery. For a computer generated tree of only 6000 nodes spanning order numbers 1–6, 5 must be added to each order number so that the orders of the root nodes (largest vessels) match and a direct comparison can be made. This is consistent with assuming that the smallest vessels in the computer generated tree correspond to vessels of order 6. Which is due to the absence of smaller vessels downstream of the smallest arteries in the in silico model.

Results

In this paper, globally optimized vessels are grown using an SA-based approach to supply a myocardial substrate, and validated through comparison with morphological data from the porcine arterial tree. We choose to examine the heart vasculature, as the structure of the large coronary arteries has been found to be similar between individuals [39] and the full arterial tree has been well characterized in porcine models [9]. For modelling the coronary arteries, we used the following parameters. (i) A tissue substrate representing an ellipsoidal human heart muscle of mass 218 g, constructed based on physiological parameters [40]. The right ventricle was assumed to take the form of a super ellipsoid of exponent 2.5 and the left ventricle was represented by a simple ellipsoid. Truncation of the ellipsoidal substrate was chosen so that the mass of the tissue corresponded to a reasonable physiological value given morphological data for ventricle thickness. (ii) Blood flow through each of the terminal segments of the tree was assumed to be constant, with each arteriole supplying an equal volume of tissue and homogeneous perfusion throughout the tissue parenchyma [41]. These assumptions greatly simplify fluid dynamical calculations for estimating the total power needed to pump blood through the tree. (iii) The metabolic cost of maintaining a given volume of blood was assumed to be 641.3 J s−1 per metre cubed of blood [30]. For convenience, each arteriole supplies a sphere of tissue with a size calculated by assuming a mean blood flow per unit mass for cardiac muscle of 0.8 ml min−1 g−1 [42]. The value taken from the literature was chosen such that it not only lay within the given error, but also conformed reasonably with both the ellipsoidal heart model, input flow and radii. (iv) The larger arteries with diameters greater than 0.01 mm were constrained to avoid penetration of the outer layer of heart tissue. This simplification differs slightly from real coronary vasculature, where progressive intrusion of arteries into the myocardium can be observed [43]. However, as the major arteries modelled by our method are far larger than the intra-myocardial vessels, a sharp cut-off is thought to provide a reasonable approximation. (v) The starting positions of the two root arteries were fixed with a total input flow of 4.16−6 m3 s−1 [44]. Relative radii of the two inputs to the tree were constrained via ; however, the relative sizes of root arteries and division of perfusion territories are determined by the method alone. (vi) The branching exponent varies throughout the coronary arterial tree, but for the larger arteries its value remains in the range 1.8–2.3. A variable branching exponent would greatly increase the computational cost of the approach, so a compromise value of 2.1 was chosen for the entire tree [45].

Coronary arterial trees containing increasing total numbers of vessels grown using the SA-based method are presented in figure 4. In real human coronary trees, there are three identifiable main coronary arteries (e.g. the schematic from [46]): left anterior descending, right cardiac artery (RCA) and left circumflex artery. The positions and relative dimensions of these are similar in most humans, with major variations observed in less than 1% of healthy individuals [47]. Trees grown using SA (figure 4) adhere well to this structure. There is a consistency in the placement of the larger arteries, although the RCA appears slightly lower, and the right marginal artery appears slightly shorter, in our models. Overall, visual inspection of the arterial structure appears extremely promising.

Figure 4.

Images showing arterial trees grown with the approach detailed here. The number of terminal arterioles is increased from 500 to 6000 (the total number of arterial segments is roughly twice this). There is consistency in the positioning of the larger arteries between the numerical method and the typical arrangement of the major arteries, suggesting that the coronary arteries may be the result of a biological process seeking the global minimum in metabolic demand.

To provide a quantitative comparison of our trees with anatomical data, the topological characteristics of the computer generated coronary artery trees were extracted and compared with morphological data characterizing the pig coronary arteries published by Kassab et al. [9]. They used a combination of corrosion casting and optical sectioning to obtain detailed morphometric data, tabulated using the Strahler (or stream) ordering scheme to denote elements of the tree of varying scale. Within this scheme, the lowest Strahler order numbers correspond to the smallest arterioles and the largest numbers refer to major vessels (for details on Strahler ordering see Methods). To directly compare arterial diameters, lengths, and branching properties of our computer-generated arterial tree with real data from pig coronary arteries, averages were obtained over all elements of the same diameter defined Strahler order.

The mean vessel diameters are shown as a function of order number for a tree comprising 6000 arterioles (12000 vessel segments) in figure 5a. Excellent agreement is found between the trees generated in silico and the morphological data. Only slight deviations from the morphological data can be seen for the smallest vessels (lowest order arteries) in the generated tree. This is likely to be due to the combination of integer order numbers and the condition that terminal sites are of constant radius. The result of this constraint is that the terminal radii will only match the anatomical data for a correct choice of the number of arterioles. Figure 5c shows the effects on diameter of increasing the number of arterioles from 500 to 2000. Agreement is generally good, regardless of the number of terminal arteries, and there is a clear trend towards matching the experimental data as simulated tree size increases. Figure 5b compares average vessel length in the model and porcine morphological data as a function of order number. For the largest arteries (high order numbers), the agreement is excellent. Although the lengths of the smaller arteries (Strahler orders <7) in the computer-generated tree tended to be overestimated, this can be easily explained by the fact that the smallest vessels are required to bridge a gap that would normally be filled by inclusion of lower order vessels in a larger simulation. As the number of generated vessels is increased, the agreement with morphological data improves (figure 5d).

Figure 5.

(a) Vessel diameter as a function of order in a tree with 6000 arterioles. Excellent agreement is found for vessels on all length scales. (b) Vessel length as a function of order number. Agreement is excellent for the major vessels (large order). The large variation seen for arterioles (lower order) is a result of early termination. Also shown are the morphological data reproduced from table 2 of [5] for easy comparison. (c,d) Similar to (a,b), but for smaller trees to highlight the trend towards the morphological data as tree size increases. (Error bars show standard errors, both axes are logarithmic.)

Previously, the best methods available for the computer generation of arterial trees struggled to recreate realistic branching asymmetry. Figure 6 shows the ratio of daughter to mother vessel radii for the largest and smallest daughter vessels as a function of order number. This provides a measure of the branching asymmetry of the tree, where small ratios indicate that branching is symmetric, while ratios approaching 1 suggest a large trunk vessel with small branches. For Strahler orders corresponding to microvascular arterioles, both the computer generated and true morphology approach 0.7, which is consistent with perfectly symmetric branching where both daughter vessels are of similar size. Agreement with the morphological data from [45] improves as the size of the computer-generated tree increases. This is not the result of any special input parameters or initial conditions. The trees are topologically and spatially randomized before SA optimization begins, and are allowed to explore the entire parameter space during optimization. The observed asymmetry is purely the result of a balance between pumping power and metabolic maintenance cost, and is a major improvement in predicting the trunk-like structure of major vessels.

Figure 6.

The ratio of daughter vessel diameters (Ds and D are the diameters of the smallest and largest daughter vessels, respectively) to diameters of parent segments D as a function of order number, showing how the tree tends towards more symmetric branching at lower orders. Agreement with morphological data reproduced from tables in the online supplement of [45] is good, if the early termination of the generated trees is taken into account, with the trend towards the morphological data as the tree size increases. Both graphs demonstrate that there are large trunks at high orders with the largest daughter vessel (b) of similar size to the parent vessel and another side artery which is much smaller (a). At smaller orders, the ratio becomes similar showing that the branchings of the smaller arteries are near symmetric. Realistic branching asymmetries are a clear advantage over other methods of generating arterial trees in silico.

Our final figures show the effect of altering the metabolic energy cost of blood per unit volume m. The largest morphological change is found in the lengths of the larger arteries (figure 7). As m increases, bifurcation symmetry is also increased in the larger arteries and as a result there is an increase in the number of Strahler orders present in the tree (figure 7). The explanation for these scaling behaviours is evident when considering the limiting cases. For m=0 the power involved in pumping the blood dominates the optimization, which leads to a large, ‘snaking’ artery with small side branches that supply the tissue. This large artery would cover the entire surface of the heart, and the configuration is equivalent to a completely asymmetric binary tree. For a large m value (or small power cost) there is a huge penalty associated with larger arteries, and so their lengths are contracted. In order to accommodate the reduction in length, the larger arteries must bifurcate more frequently and symmetrically. Additionally, the high volume cost causes the trunk artery to minimize its total length, resulting in a much straighter path across the tissue. Less extreme examples of this behaviour can been seen in figure 8, with meandering arteries for small m and straight arteries for large m.

Figure 7.

(a) Diameter as a function of order number for trees with 1000 vessels. Decreasing m, which describes the relative energy cost of an amount of blood and the power required to pump it, has little effect on the agreement of the diameters with morphological data. (b) For lengths however there is an obvious effect in the larger arteries, with regimes of high pumping cost being more accurate. The primary optimization for high pumping cost then is to increase the length of the largest arteries. (c,d) The main effect is a change in the asymmetry of the branching of the largest arteries—for large m, the branches are more symmetric than for small m. As m becomes very small, the limiting behaviour is broad trunks that wind around all the tissue, with a large number of very small offshoots that supply blood in the direct vicinity of the large vessel.

Figure 8.

Example trees generated with different values of m, which changes the relative weight of the pumping power to cost of maintaining blood in the optimization. For small m (corresponding to small hearts), vessels in the trees wind around—this is because there is little penalty to make a single wide vessel that curves to supply blood, rather than bifurcating. For large m (corresponding to large hearts) the vessels travel as straight as possible.

The change in m can also be interpreted as a change in length scale as follows: once the large vessels have been excluded from the tissue and all tissue is supplied, the remaining cost function that is optimized has the form Now, make the transformations , . Then the cost function becomes As the optimum in the cost function is the same independent of a multiplicative factor that acts on all terms, then we can absorb a factor of 1/A3 into the cost function to obtain Identifying a new the cost function now has the same form. As changing m is equivalent to changing the length scale, these results suggest that there are likely to be structural differences between species of different sizes, as the power required to pump blood becomes relatively more important than the metabolic demand to maintain blood volume in small vessels. In the absence of morphological data, visual comparison of the coronary arteries tentatively indicates that vessels meander around in smaller species [48] and that vessels are straighter in larger species [49].

Discussion

We have developed a powerful and universal method for growing arterial trees in silico, which is capable of identifying the near globally optimal configuration of arteries for arbitrarily shaped tissues with heterogeneous blood supply demands. As input, the method only needs information about the tissue structure and the entry point positions of the largest arteries. From this information, the approach generates morphologically and structurally accurate coronary arterial trees at almost every length scale. This is a significant improvement on previous optimization methods, which failed to reproduce the consistent structure found in the coronary arteries. We have shown that the method improves with the number of vessels modelled, so that, as computing power increases, there is a systematic improvement in the accuracy of the generated trees. To our knowledge, no other method can generate realistic arterial trees that closely match morphological data by taking only the shape of the tissue as input, and claim systematic improvement in the generated trees with increased computational power.

We expect that our method could have several useful applications. Our first application is to use cardiac and partial cerebral vasculature structures (e.g. MCA territory) computed using this method as input to models of embolic stroke and other infarctions, as these models require detailed vasculature structure over a range of length scales which are not available to imaging techniques [50-52]. This will require further validation of the algorithm against cerebral arterial data. Such vasculature would be downstream of the Circle of Willis, a source of major anatomical variation and probably not a structure reproducible by the current algorithm. Computational models of stroke combined with Doppler ultrasound have potential to provide further information regarding embolic burden during major operations [53].

There are several other potential applications. Models of arterial trees generated by our method may help to improve the interpretation of medical images though advanced image segmentation techniques. Vessels identified through automatic segmentation techniques can be connected via algorithms such as the one presented [54]. In addition, segmentation can be limited to the location of a reduced set of bifurcation points, and the algorithm used to fill in any missing vessels which connected them [55].

We also speculate that the algorithm could be of use in designing the structures of vasculature for artificial tissues. Once the very difficult and intricate process of making networks of vessels in artificial tissue has been achieved [56-58], the opportunity to optimize or inform their design will be available. A common problem during the growth of artificial tissue is that regions of cells can die due to lack of nutrients and oxygen. For instance, an artificial skin graft may be optimized for increased healing, by having its vasculature designed such that there is higher overall perfusion with minimal loss of useful tissue. Even more speculatively, artificially grown organs may have a vasculature designed to minimize their impact on the cardiovascular system.