International Journal of Hepatology

Volume 2012 (2012), Article ID 357687, 17 pages

doi:10.1155/2012/357687

Research Article

Analysis and Algorithmic Generation of Hepatic Vascular Systems

Lars Ole Schwen¹ and Tobias Preusser1,2

¹Fraunhofer Institute for Medical Image Computing MEVIS, Universitätsallee 29, 28359 Bremen, Germany

²School of Engineering and Science, Jacobs University, Campus Ring 1, 28759 Bremen, Germany

Received 5 April 2012; Accepted 16 May 2012

Academic Editor: Mathieu Vinken

Copyright © 2012 Lars Ole Schwen and Tobias Preusser. This is an open access article distributed under the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Abstract.

A proper geometric model of the vascular systems in the liver is crucial for modeling blood flow, the connection between the organ and the rest of the organism. In vivo imaging does not provide sufficient details, so an algorithmic concept for extending measured vascular tree data is needed such that geometrically realistic structures can be generated. We develop a quantification of similarity in terms of different geometric features. This involves topological Strahler ordering of the vascular trees, statistical testing, and averaging. Invariant features are identified in human clinical in vivo CT scans. Results of the existing “Constrained Constructive Optimization” algorithm are compared to real vascular tree data. To improve bifurcation angles in the algorithmic results, we implement a postprocessing step calibrated to the measured features. This framework is finally applied to generate realistic additional details in a patient-specific hepatic vascular tree data set.

1. Introduction

The most important link between the liver and the rest of the organism is the blood flow through three vascular systems [1, 2]. These include two supplying systems, the portal vein (PV) providing venous blood drained from the digestive system and hepatic artery (HA) providing arterial blood, as well as one draining system, the hepatic vein (HV). The bile duct is a fourth vascular system which transports the secreted bile from the bile canaliculi out of the liver into the cystic and common bile duct finally leading into the duodenum.

On the one hand, for the functioning of the liver metabolism a proper blood supply and drainage are essential. In fact, the liver receives about of the cardiac blood output, which amounts to about mL/min per g net liver weight [3]. On the other hand a variety of pathological conditions result in impaired blood flow conditions. For example cirrhosis decreases total hepatic perfusion while increasing the fraction of arterial blood [4]; metastases also lead to an increased arterial fraction [5]. Also, it is known that the blood flow plays a major role in the regeneration capabilities of the liver. In summary the blood flow plays a central role in the liver function and understanding blood flow, and its regulation can be seen as a key to understanding liver physiology and pathology.

Biophysical modeling and simulation have become powerful tools in analyzing and understanding the behavior of complex dynamical systems or to predict future states of such systems without actually performing the corresponding experiments. Such modeling and simulation of physiological processes in the human body do not only have impact on the basic science of understanding life but also direct consequences ranging from pharmaceutical developments up to improved quality in surgical interventions.

In order to properly model and simulate the metabolic function of the liver, it is crucial to have an appropriate model of the blood transportation systems. In a multiscale model [6], the range between the whole organ (e.g., pharmacokinetic models such as [7, 8], and the individual lobules [9–11] is covered by the vascular systems).

For modeling and understanding physiological processes, the level of detail needed in the vascular structures depends on the spatial heterogeneity and the scale of the process being considered. For planning liver resection surgery (see, e.g., [12]), one main task is to determine the territories supplied by large vessels. For this purpose, details available from in vivo imaging are sufficient. Certain liver diseases were observed to be spatially inhomogeneous, for example, chronic hepatitis and cirrhosis [13], fibrosis [14, 15], and steatosis [16]. In this case, a multiscale model considers representative volume cells, consisting of groups of lobules, of the liver sufficiently small such that their properties can be assumed to be more or less homogeneous. The vascular structures then need to be sufficiently detailed to properly reflect supply and drainage of such representative volumes. Lobules, their internal sinusoidal network and a potential zonation of hepatocytes, can and should be viewed on a separate scale in a multiscale simulation framework.

However, current imaging and image processing techniques are not capable of resolving the full vascular system in human livers. Thus, studying the vascular structures from coarse to fine scale, from the portal vein to the finest hepatic units, the sinusoids, is not directly possible.

In vivo CT scans of the human liver provide a resolution of a few hundred micrometers. Ex vivo corrosion casting is an established technique requiring huge skills. If not prepared in situ, however, a deformation of the organ needs to be accepted. Corrosion casts are obviously smaller than whole bodies and permit higher doses of radiation and thus provide better image data, but only few of them are available.

Microscopy generates optical images at much higher resolution which are capable of resolving intralobular sinusoids, but the images are essentially only 2D. Reflectance confocal microscopy [17] or fluorescence confocal microscopy allows optical imaging a few  μm below the surface without physical sections. Histological serial slices can be created and scanned optically, but the cutting process introduces physical artifacts that can only be remedied to a limited extent by image processing (registration, inpainting).

To bridge this resolutional gap, algorithms for generating vascular systems can be applied, thus allowing to build full three-dimensional vascular trees that permit modeling and simulation of liver blood flow from the coarsest to the finest scales. However, to assess and improve the quality of such algorithmically generated structures, that is, to determine whether they are sufficiently realistic for modeling purposes, requires better techniques than merely a visual comparison.

The purpose of this paper is twofold. On the one hand, we describe a method to quantify the difference between structural and geometric features of different vascular structures in Section 3. We will compare a significant number of real contrast-enhanced liver CT data sets and several corrosion casts. Thus we will be able to quantify the similarity of hepatic vascular structures over a large group of individuals. On the other hand, we use these findings to improve an existing algorithmic method to more closely resemble hepatic vascular structures found in humans in Section 4. After this calibration, the method is finally applied to extending measured vascular networks by more details than we can obtain from the image data. While we restrict the description and analysis to one of the vascular systems of the liver in this work, the framework is generic and can be applied to other organs as well. Moreover, the extension to more than a single vascular structure is possible and currently being investigated.

2. Review of Related Work

Vascular structures have been studied for a variety of human organs such as lungs [18, 19], hearts [20, 21], retinas [22], and livers [23] as well as in whole rats [24] and animal organs such as cat lungs [25], pig hearts [26–28], rat hearts [29], and rat kidneys [30]. Even though this section provides a relatively detailed overview on existing work, it is by no means meant to be an exhaustive literature review.

Borrowing vocabulary from graph theory (see, e.g., [31]), many vascular systems can be described as trees consisting of nodes and edges. In the tree hierarchy we will use the terms “parent” and “daughter” edge and permit only strictly bifurcative trees; that is, each edge has exactly zero or two successors, with a single edge incident to the root node. In most cases we are not merely interested in the topology but also the geometry of vascular structures, so nodes are typically assigned a geometric position and edges a radius. As the edges play a more important role in our considerations than bifurcations/nodes, we assign geometric quantities to edges in our analysis. To avoid confusion with “Couinaud segments” going back to [32], we avoid the term “segment” when dealing with vascular networks.

Topological Classification

The simplest topological classification for edges in a tree is the generation number [33] (also named “bifurcation level” [20], “level” [24], “Weibel order” [34]) which, however, is incapable of distinguishing the “importance” of edges in terms of their supplied subtree). The most commonly used classification is the “Strahler” order (e.g., in [19, 35, 36] and many more articles) originally used to classify rivers [37, 38], a measure nondecreasing from leaves to the root. More precisely, edges leading to leaves are assigned order ; parent edges are assigned the larger order of its daughters (if they differ) or the two daughters’ order plus ; see also Figure 2. While it may seem problematic at first glance to start at the smallest edges with highest measurement uncertainty, this scheme is actually very robust. Compared to the more intuitive generation counting, Strahler ordering is a “major improvement […] because it takes into account the asymmetric branching pattern” [39]. Two edges near the root have the same generation number even if one is “small” and leads to a leaf node and its sibling supplies a “large” subtree; see also [25].

The Strahler ordering scheme has been refined by also including radius information (e.g., in [34, 39, 40]). The order at a bifurcation is then only changed if there is a sufficient increase of radius. In our analysis we chose not to use this type of levels because it is no longer strictly topological and because the radius data in our data sets is not very accurate [41].

Geometrical Features

The two most obvious geometrical features for tubular edges are their radius and length [42–44]. At bifurcations, different angles can be considered: the angle between the continued parent edge and either of the daughters [28, 45, 46], the angle between the two daughters [29], the inclination of one of the edges with respect to the plane spanned by the other two (cf. [43]). In [47], Voronoï cells of the leaf nodes are considered, and their asphericity is measured. They can be combined to supplied volumes of edges, but in any case the resulting general polyhedra are not easily comparable.

From the previous “natively” geometric features, derived geometric features can be computed. From the radius, we can compute the radius ratio to the parent [28], the radius asymmetry between larger and smaller daughter radius (“branching ratio” in [20], “bifurcation index” in [28, 48]) or vice versa [42]; cross-section areas and analogous cross-section area ratios [20, 48] and asymmetries [49]. Moreover, a “bifurcation exponent” can be computed for where is the parent radius and are the daughter radii; see, for example, [48, 50]. The exponent is an optimal trade-off between power dissipation for moving the blood and metabolic cost for maintaining the blood vessel wall if laminar Poiseuille flow can be assumed [51, 52]. While this is a property of a single vascular edge, flow (mass) conservation at bifurcations allows considering to be a property of a bifurcation. An alternative value is used in [53] as it minimizes vascular wall material [54] and was measured in [55]. The value minimizes the reflection of pulse waves at bifurcations [56].

Based on the purely geometric features, derived features have been investigated, computed according to additional assumptions. These include pressure drops, inertance, compliance [57], perfusion heterogeneity and asymmetry of supplied volume [58], edge flow, exit pressure and transit times [59], additional compliances [60], pressure profiles [61] and perfusion inhomogeneity in surrounding tissue [61, 62], and wall shear, internal pressure, and circumferential tension [63]. Theoretically optimal angles are discussed in [64, 65].

Various geometric features were determined depending on the bifurcation order/generation number [20, 29, 44] or on the Strahler order [18, 19, 30, 66] and on radius-based Strahler orders [18, 34, 58]. Also, the dependency of various features on edge radii, which are a noninteger quantity was investigated [33, 59, 67–69]. Also, relations between different features were discussed in [21, 44, 50], to name just a few. Software for analyzing vascular trees is presented in [29]. Self-similarity and other fractal properties were frequently observed and investigated; see, for example, [70–75].

Our geometric analysis will be based on the strictly topological Strahler order scheme and only take into account purely geometric features so that there is no dependency on additional physical or physiological assumptions.

Constructive Optimization

Different types of algorithms for generating vascular structures have been presented in the literature. One class of algorithms is based on constructive optimization, that is, determining vascular networks satisfying physiological optimality criteria (e.g., minimal intravascular volume) for a given supply requirement. This requirement is typically given by a set of pseudorandomly distributed leaf nodes interpreted as connections to lobules on the next finest geometric scale. The most prominent approach here is “Constrained Constructive Optimization” (CCO) [61, 69, 76, 77]. These were originally developed for supplying a convex domain (organ); an extension to the nonconvex case is presented in [78]. In contrast, “Global Constructive Optimization” [79, 80] performs a multiscale optimization finding an optimal tree for all leaf nodes at the same time.

The basic principle of CCO is to construct a strictly binary tree by adding one leaf node at a time to an initial tree, each time introducing an optimal bifurcation. Thus, CCO can be seen to be driven by the assumption of equal in- or outflow at all leaf nodes representing constant supply/drainage for each lobulus. Moreover, at bifurcations the radii are balanced such that the flow resistance according to the Hagen-Poiseuille law [81] is equal for both subtrees. This finally results in equal exit pressures at the leaf nodes.

Deterministic Geometric Construction

Other authors use a deterministic geometric construction using Cartesian [82] or hexagonal [53] cell-based models or constructions inspired by self-similar (fractal) or area/space filling objects [43, 83–86]. Such approaches benefit (in terms of algorithmic complexity) from avoiding optimization problems. The resulting vascular networks, however, suffer from an artificial overall structure, which is even visually perceptible, so that they are not appropriate for general models.

Angiogenesis-Based Construction

Yet another approach is modeling angiogenesis, the actual process by which vascular structures grow. This type of approach requires more involved models and algorithms than basic optimality conditions. Earlier results in this area [46, 53, 82, 87] exhibit a visually artificial structure or “somewhat stylized appearance” [82, Figure 1]. More recent work includes [88, 89], the latter combining introducing new vascular edges due to angiogenesis with subsequent geometric optimization of the vascular tree. Grid-based methods typically produce visual geometric artifacts reflecting the grid used. To us there seems no easy way of introducing parameters in these algorithms for being able to calibrate them to better match geometric features measured in the vasculature of human livers. One could combine angiogenesis models with geometric parameters by changing the way how new vascular edges form. Instead of only considering gradients of angiogenetic factors generated by ischemic cells, also properties of the existing vascular edge for which a new bifurcation is to be introduced, on the branching angles of that bifurcation or on other properties of the existing vascular structure, could be taken into account. This would, however, involve additional assumptions for the model and parameters in the algorithmic implementation which are not easily observed experimentally.

In this paper we base the construction of vascular structures on CCO as this is an established technique and avoids the drawbacks of deterministic or angiogenesis-based methods.

Figure 1: A bifurcation with the radii, lengths, and angles considered for our geometric analysis.

(a)

(b)

(c)

Figure 2: The (a) sketch shows the definition of Strahler* orders: the order decreases towards the root if two edges of the same order meet, else the minimum is used, and the values are shifted such that the root edge has order 0. Images (b) and (c) show Strahler* orders (HSV color coded from blue = 0 to red = 5) for two of the clinical PV data sets.

3. Morphometry of Hepatic Vascular Systems

In this section, we aim at quantifying the similarity between different vascular networks. We first describe the specimens used for this study. A similarity measure is introduced and results are presented.

3.1. Data Acquisition

For the geometric analysis, two types of data were used. Corrosion casts of 6 human PVs from [86, 90] and described therein plus one additional PV data set were used. Only 3 corrosion cast HV data sets were available from these studies, not enough for a robust statistical analysis, so we ignored them for our purposes. CT scan image data was segmented and skeletonized by a semiautomatic procedure described in [23]. This yields a graph representation of the vascular structures with curved edges of nonconstant radius. Moreover, 167 PV and 165 HV vascular trees were obtained from clinical in vivo CT scans using contrast agent, for which image preprocessing the same procedure was applied. Edges of the HA scans are thinner, and the extracted data is less reliable than for the other structures, so they were excluded from our consideration. As the HA is essentially parallel to the PV (and bile ducts) and is responsible for roughly one-third of the PV flow (albeit of different composition), a separate analysis of these was considered of minor importance. The bile duct is not considered here.

Depending on the quality of the image data and the preprocessing steps the graph representations were further modified because different properties of the tree may violate the strictly bifurcative property needed for our analysis. Isolated nodes, monofurcations, or multifurcations (edges with more than two daughters) may be present. Loops may exist in two different flavors, either as two distinct paths between two nodes in edge-parent direction (improper loops) or as nontrivial paths in edge-parent direction from a node to itself (proper loops). Improper loops contain at least one node with two incident parent edges. The root node may have more than one incident daughter edge. Thus a correction step is needed to make the extracted graphs suitable for our further analysis.

In this correction step, existing trivial edges (initial node equals terminal node) are removed. Nodes are then tested for the number of incident parent edges. If multiple parent edges coming from the same initial node are found, their radii are combined such that the cross-section area is preserved. For multiple parent edges with different initial node, one with maximum radius is preserved, and the others are removed. Next, monofurcations are removed and multifurcations are split by again introducing trivial edges. Removing monofurcations requires averaging the radii of the two edges involved. As these radii are already estimated radii from the image data, this is only a slight modification comparable to inaccuracies introduced by the skeletonization. Introducing trivial edges does not affect the geometry of the vascular structures at all.

Next, we test for the existence of proper loops which, however, did not occur in our data sets. Any isolated nodes are removed and we finally copy the subtrees for all root nodes found, introducing trivial edges if necessary such that root nodes only have one incident daughter edge.

When dealing with corrosion casts, we must note that deformations of the corrosion casts not produced in situ particularly affect bifurcation angles for the larger edges. Also we must be aware of the fact that radius information is not very accurate in the graphs [41] and resulting trees. This limitation is mainly due to the fact that the vascular structures are only a few voxels wide in the image data. Thinner vascular segments—if they are present in the corrosion cast or filled by contrast agent at all—are not visible in the image data. The algorithms used in the semiautomatic image processing have a high hardware demand and are currently limited to processing images of about 70 megavoxels. Simpler, fully automatic image processing techniques do not produce satisfying results because the connectivity of the individual networks, but separateness of separate networks, if present in one image data set, cannot be obtained reliably.

3.2. Topological and Geometrical Quantification

As we have different sources of image data (contrast enhanced CT, CT of corrosion casts) with different image quality, we will be comparing vascular structures with potentially different level of detail resolved. Thus, Strahler orders of the respective root edges differ, and edges of a given Strahler order play a different role in different trees. We hence introduce a the notion of , assigning to each edge the difference between the Strahler orders of the root and the current edge; see Figure 2.

For the analysis of the vascular structures we will focus on the following geometric features for each edge in the tree. For an overview of the many possible geometric properties and of what has been considered in the literature we refer to Section 2.

Radii

Each edge has a radius , a ratio , except for the root edge. Moreover, for each nonterminal edge there is a ratio . Moreover, the bifurcation exponent described previously is considered a property of the parent edge.

Lengths

Each edge has a length , and are defined as the obvious analogues of and .

Angles

A bifurcation is fully described by three angles. Without loss of generality, assume the two daughter segments lie in the plane such that their angle bisector is aligned with the axis. Then let be the angle between the two daughter edges, let be the angle between the projection of the parent edge to the plane and the axis, and let be the inclination of the parent edge to the plane; see Figure 1. Note that these angles can be computed as properties of the parent without distinguishing the two daughters.

If, at a given bifurcation, a geometric feature cannot be computed, the value is not considered for the evaluation later on. Hence, the number of bifurcations actually considered for different features may vary.

Fractal properties, despite being interesting for individual, were not considered for a comparison later on. Furthermore, features computed from the purely geometric features were not considered such as to avoid redundancy and/or dependency on physical or physiological assumptions. These could be interesting in later studies.

3.3. Quantification of Similarity

In order to obtain a scalar measure of similarity within one set of vascular trees (e.g., measured trees) or between two sets and of vascular trees (e.g., measured versus generated ones), we perform six steps in which we successively compute

(1)Strahler* orders for each segment of each tree and the applicable geometric features,(2)histograms of the geometric features described previously per Strahler* order for each tree,(3)a binary decision of similarity for each geometric feature and each Strahler* order between each pair of trees (within either  or ;  we are then interested in the ratio (number of pairs classified similar) over (number of all pairs), still for each geometric feature and each Strahler* order,(4)a weighted average of these ratios over the Strahler* orders for each geometric feature, (5)weighted averages for the radius, length, and angle features; (6)a weighted overall average.

Statistical Interpretation

The computation in Step 1 is straightforward and described in Section 3.2. For the remaining steps, we consider the values as samples of an unknown underlying probability distribution, for which we compute the (empirical) cumulative distribution function (CDF) in Step 2. Note that we do not use probability densities because the CDFs are directly used in subsequent statistical testing and because they avoid binning when plotting the data.

Testing for Similarity

For two CDFs corresponding to the same geometric features of edges of fixed Strahler* order of two different vascular trees we use the two-sample Kolmogorov-Smirnov test [91, 92] (KS test), to check the null hypothesis that the two CDFs have the same underlying distribution. In Step 3 we classify two trees as “similar” (in terms of this feature and Strahler* order, written as ) unless the test rejects the null hypothesis at significance . Note that there is no particular reason for choosing this value here; other than that is is very common [93]. In comparison to other statistical tests, for example, the Cramér-von Mises test [91], the KS test has the advantage of being independent of scaling in the values of the underlying distributions.

Let us first consider the case of one set of vascular trees, where we are interested in the similarity among, for example, the clinical PVs. This will also be referred to as the single-population case. A similarity ratio for a fixed geometric feature and Strahler* order is then computed as

Where  denotes classification as “similar” in terms of and in the sense defined above.

If we are interested in the similarity between, for example, a set of measured vascular trees and a set of generated trees, we define a similarity ratio as

and refer to this as the two-population case.

Note that due to the nature of statistical testing, we cannot expect these similarity values to be larger than . This also holds for the weighted averages in the following.

Averaging over Strahler* Orders

In order to average over the Strahler orders in Step 4, but still considering a fixed geometric feature , we wish to take into account the number of edges present in the tree at different Strahler* orders. Since different vascular trees in the population considered may have different level of detail, we restrict the averaging to those Strahler* orders () below which at least half of the trees actually have edges. The weighted average is obtained as

In the two-population case, the maximal Strahler* order is the maximal order at which at least half of the trees in both populations actually have edges and the averaging weights are the geometric means of the same values as before. The weighted average for a fixed geometric feature is obtained as

Averaging over Features

For averaging over features in Step 5 in the one-population case, we use the arithmetic mean to quantify an average similarity. In the two-population case, we use the values as averaging weights so that invariant features in the measured trees are weighted stronger than noninvariant ones, such that we can quantify “similarity where it is expected.” The two numbers should hence not be compared directly. Averages over a set of features are defined as

in the one- and two-population cases, respectively.

Finally, in Step 6 we average over different sets of features (related to radii, lengths, and angles, resp.) in such a way that the different numbers of features per set are taken into account. In the one- and two-population cases we average as follows:

Implementation of the KS Test

The two-sample two-sided KS test [91, 92] involves three steps, computing the maximal difference between the two CDFs, rescaling it by the factor depending on the two sample sizes and evaluating the KS distribution to check whether it is smaller than the significance level. Evaluating the KS distribution requires a brief discussion.

Tabulated vales [94, 95] (e.g., [96] based on Monte Carlo simulations) for fixed significance levels were used traditionally, which lacks the possibility to use arbitrary significance levels if this is desired. The limiting KS distribution [97, 98] for large sample sizes

is the same in the one- and two-sample cases [99] with appropriate scaling factors for the sample size.

Evaluation of the distribution in the one-sample case for small samples is discussed in detail in [98]. A review of current software capable of doing two-sample KS tests reveals that octave [100, 101] and root [102, 103] use the limiting KS distribution regardless of sample size. For small samples, R [104] uses [105] for the one-sample test and an undocumented method for the two-sample test, which seems to be a generalization of [106] to samples of different cardinality.

We performed own Monte Carlo Simulations creating many small samples with (uniformly distributed or the sum of two uniformly distributed) random numbers, obtained using the Mersenne Twister [107], and computing the KS values for them. These should be uniformly distributed in , regardless of the distribution of the sample if all samples are drawn from the same distribution. For the limiting distribution, we confirmed the known “overestimating effect” [108]. Using the one-sample implementation of [98], available as a supplement to that paper, in the two-sample case showed only a slight overestimation of the distribution. The code used in R, however, produced accurate results. The numerical computation of the CDFs and their differences should be performed with at least IEEE double accuracy.

3.4. Results

Results of the topological analysis are shown in Tables 1 and 2 show that how many of the vascular trees considered had edges at different Strahler* orders and how many edges there were on average, as well as a histogram of the connectivity between different Strahler* orders. The maximal Strahler* order in the clinical CT scans and the corrosion cast data sets is even though the corrosion casts contain more edges, in particular of high Strahler* order, which is also reflected in the connectivity histogram. For a strictly binary tree with leaf nodes, these histograms would show only connections between subsequent levels, approximately for the connection between levels and .

Table 1: Average number of human PV and HV trees considered for statistical analysis that have edges of Strahler* orders to , along with the average number of edges per tree where such edges are present.

Strahler* order

PV
No. of trees

avg. no. of edges

1

HV
No. of trees

avg. no. of edges

corr. PV
No. of trees

avg. no. of edges

Table 2: Histogram (in percent) of transitions at bifurcations between different Strahler* orders present in the human clinical PV and HV ((a) and (b)) and corrosion cast PV (c) considered.

(a) Clinical PV

↗

(b) Clinical HV

↗

(c) Corrosion cast PV

↗

The cumulative distributions of the bifurcation exponents for all the data sets considered are shown in Figure 3. These turn out to be very similar for all Strahler* orders but differ between the clinical PV and HV as well as between clinical and corrosion cast PV data sets. Radius data, in particular for thin radii at high Strahler* orders, is not very accurate, so part of these distributions may be due to imaging and image processing inaccuracies. However, larger radii are more robust and exhibit the same distribution.

(a)

(b)

(c)

Figure 3: Overall histograms of the bifurcation exponents for our clinical PV (a), HV(b), and corrosion cast PV (c) data sets for different Strahler* orders.

The results of the similarity analysis for our data sets are listed in Table 3. We can observe that the similarity among HV trees is slightly smaller than the similarity among PV. The similarity among the corrosion cast PV data sets is much smaller; this is probably due to the small number of specimens considered and the different conditions under which they were manufactured. Comparing clinical and corrosion cast PV, radii and absolute lengths are relatively different whereas angles are only slightly less similar than among the clinical PV. Moreover, Table 3 shows that the similarity between PV and HV is comparable to the similarity within each of the two populations, except for the angle features.

Table 3: For different geometric features, the table lists the one-population similarity within the human clinical PV, HV, and corrosion cast PV data sets. Moreover, the two-population similarity between the clinical and corrosion cast PV as well as the clinical PV and HV data sets is listed. Unlike for the individual features, the numerical values of the averages in the one- and two-population cases cannot be compared due to the different averaging formulas (7) and (8).

Feature
Similarity
Similarity

PV
HV
corr. PV
PV/corr. PV
PV/HV

Radius

Radius decrease

Radius asymmetry

Bif. exponent

Lengths

Length decrease

Length asymmetry

Angle

Angle

Angle

Radius average

Length average

Angle average

Total average

The averaging procedure over Strahler* orders for two fixed geometric features in the one- and two-population cases is shown in Table 4. The plots of the CDFs also illustrate that not only the average difference but also number of samples (visible as vertical steps in the graph) influences the similarity percentages and —which is what the statistical test is meant to take into account. The decrease of the average radius with increasing Strahler* order, see, for example, [86], and other dependencies of geometric quantities on the topological orders can be confirmed in the detailed results but were not within the scope of the statistical comparison in this study.

Table 4: The tables show the invariance analysis for two example features, edge radii and the bifurcation exponent . In (a), we consider the one-population case of clinical PV, and (b) treats the two-population case clinical versus corrosion cast PV. The average number of trees as well as the average number of edges per tree is those considered for the respective feature and may thus vary due to geometric degeneracies. We moreover list the averaging weights from (3) and (5), respectively, to determine the average invariance from the individual invariances for the different Strahler* orders . The plots show the empirical cumulative distribution functions of the features, red for the clinical and green for the corrosion cast data sets.

(a) Detailed similarity analysis for clinical PV

Strahler* order

No. of trees

avg. no. of edges

No. of trees

avg.  no. of edges

n/a

(b) Detailed similarity analysis for clinical PV/corrosion cast PV

Strahler* order

 No. of trees

avg.  no. of edges

 No. of trees

avg.  no. of edges

n/a

4. Algorithmic Generation of Vascular Structures

In this section, we first discuss our implementation and application of standard CCO. We then propose improvements and quantify their effects finally showing applications of the improved scheme.

4.1. Constrained Constructive Optimization

First, we consider CCO as described in the literature [61, 69, 76, 77]. Since livers are of nonconvex shape, we use a penalty approach [109, Section 4.2] inspired by [78]. We penalize the location of bifurcation nodes lying outside the organ domain by the squared distance from the organ scaled by a constant such that our objective function becomes

with if the unit of lengths is mm. The node-based penalty approach clearly only considers the two endpoints of a cylinder and ignores its radius, but we consider it appropriate as these cylinders are only meant to be an approximation of vascular edges which may be curved in reality. We choose but investigated different values of lambda (see also [110]).

The main task in CCO is to add one new leaf node to an existing vascular tree, finding an optimal new bifurcation. The topology of the optimal connection is not clear a priori; thus different new connections need to be tried and optimized geometrically. The geometric location of the newly introduced bifurcation is optimized using a gradient descent method with Armijo line search [109, 111]. The gradient is not evaluated analytically but using a centered difference quotient. This is because the radii of edges change, and not only locally so, as soon as a bifurcation is moved, as radii are rebalanced as described in Section 2. Leaving this out or only rebalancing radii after optimizing the position for fixed radii produces considerably different results.

It is known that a constant viscosity of blood can only be assumed for radii greater than  m due to the Fåhræus-Lindqvist effect [112]; see also [60] in the lung context. This is relevant because radii in a full vascular network range down to  m [50]. Changes in the effective viscosity can easily be integrated in the CCO procedure by evaluating the viscosity, which is needed when computing flow resistances for determining radii [63]; refer to [113] for a formula:

with Pa s,  m in the range  m. In order not to have a discontinuity in , we add a linear transition to in the range .

Computational Workload

To limit the computational workload for doing so, only a fixed number of nearby existing edges are tried (we heuristically determined that the 40 edges with closest midpoint to the new leaf node practically always contained the optimal new connection). Moreover, we first optimize with a rough stopping criterion (square root of the final tolerance), select the 20 best candidate topologies, then optimize these fully and select the optimum.

Without considering the Fåhræus-Lindqvist effect, the computational complexity for adding one new leaf node to an existing vascular tree with nodes is , where factors and higher logarithms are neglected here. This assumes that the number of gradient descent steps and Armijo steps in each gradient descent step is limited (we observed this in practice but cannot provide a rigorous proof for this claim) so that the number of radius rebalancing loops and cost function evaluations is bounded by a constant. The latter can be performed in if appropriate elaborate caching of unchanged values is used. So the total workload for generating a vascular tree with leaf nodes is . If, in addition, viscosity depends on radius due to the Fåhræus-Lindqvist effect, the radius rebalancing becomes more expensive because caching can no longer make use of relative radii of subtrees.

Applying CCO in a Calibration Procedure

Our goal is not to create vascular trees from scratch but to generate a higher level of detail in patient-specific given data sets so that in particular the coarse anatomic details are fixed and need not be generated algorithmically. To check how good generated vascular trees resemble measured ones, we thus proceed as follows.

(1)A given data set is loaded and preprocessed (made strictly bifurcative) as described in Section 3.1. (2)The vascular tree is pruned to containing only coarse anatomy by preserving the root and recursively its daughters as long as the initial node of the edge is outside the organ mask or the radius is larger than times the maximum radius of the tree (see also Figure 4). The first condition makes sense because in the given data sets, the vascular structures start from outside the segmented organ mask. (3)Next, leaf nodes at pseudorandom positions (based on [107]) with a minimal distance between them are generated. Due to a pruning step later on, we here generate times the actually desired number of leaf nodes (which in turn is chosen uniformly distributed in , thus being in a similar range as the data sets we wish to compare to). (4)One by one, these leaf nodes are added to the tree using CCO as described previously. (5)Finally, the leaf nodes of the “coarse anatomy” tree and additionally all generated leaf nodes are pruned from the tree, reducing the bifurcation level of the resulting data set. This is meant to reflect the fact that the clinical CT images were obtained at limited resolution, leaving out a lot of fine anatomy, and in particular not having leaf nodes near the boundary of the organ. The factor above was determined heuristically to compensate for this pruning. Clearly more than one bifurcation level of nodes is missing in the clinical CT data sets (a human liver has 1 to 1.5 million lobules [114]), but in order to limit the computational workload for the calibration, we consider this one additional set of leaf nodes sufficient. (6)For each clinical CT data set, we generate different CCO trees (different random seed leads to different sets of leaf nodes).

Figure 4: For the “dense” example in Figure 2, the (a, d) shows the preprocessed clinical PV tree before (a, b, c) and after (d, e, f) pruning to the coarse anatomic structures. Two resulting CCO trees before and after angle postprocessing are shown in the (b, e) and (c, f), respectively.

Postprocessing for Geometric Analysis

In order to avoid numerical artifacts due to very short segments in the subsequent geometric analysis and comparison to the measured vascular trees, very small edges shorter than a given threshold were contracted to length , partially resulting in Inf and NaN values being ignored in the analysis.

4.2. Results and Shortcomings of Standard CCO

In order to verify that minimizing intravascular volume in the CCO cost function indeed produces the best overall results, we apply the CCO procedure described previously for to generate different CCO trees for each of of the clinical PV data sets. The effects of different on the appearance of the vascular trees generated by CCO were discussed in [47, 115]. As we can see from Table 5, produces better overall results than the smaller or larger values of considered here. Radii, for which we consider the measurements inaccurate, actually are least similar, but lengths are most similar for . Angles are more similar for smaller , but these are postprocessed later on. Let us point out that the averages obtained for differ slightly from the values presented in Table 6 because only a subset of the clinical PV trees was considered here.

Table 5: The best overall CCO results before postprocessing are achieved for a radius exponent in the cost function (12).

Cost exponent

Radius average

Length average

Angle average

Total average

Table 6: In this table we compare the similarity of vascular trees obtained by standard CCO (sCCO) and improved CCO (iCCO) to the clinical CT data sets. The similarity values for individual features, but not the averages, can be compared to the HV and PV single-population similarities in Table 3.

Feature
PV/sCCO
HV/sCCO
PV/iCCO
HV/iCCO

Radius

Radius decrease

Radius asymmetry

Bif. exponent

Lengths

Length decrease

Length asymmetry

Angle

Angle

Angle

Radius average

Length average

Angle average

Total average

Results of the Geometric Analysis

Results of the comparison of generated and measured clinical PV and HV trees are shown in Table 6. From these data and the underlying cumulative histograms (not shown), a number of observations can be made.

Absolute radii are in a comparable but larger range than in the measured data sets, maybe indicating that the Strahler* ordering scheme is not the optimal topological classification. The radius decrease factor is bounded by in the generated data sets which is not the case in the measured data. Radius asymmetry compares poorly, generated bifurcations are generally more asymmetric than the measured ones. The bifurcation exponent is fixed to in the algorithm, which is clearly an artifact compared to the measured distributions of (see Figure 3). However, the radius data from the measurements is not particularly accurate [41], so a high similarity between generated and measured radius-related features is not to be expected.

Absolute lengths are in a comparable but larger range than the measured data sets, similar to the absolute radii. Length decrease factors and length asymmetries are nicely similar between generated and measured trees.

Angles are nicely similar; in the generated data sets are significantly smaller than in the measured data. The same is true for the , which is not surprising because the optimality criterion causes bifurcations to be flat. Note, however, that not all bifurcations are flat because CCO is not a global optimization procedure.

4.3. Improved CCO

As for further improving the standard CCO as described previously, we first tried to obtain a higher similarity of the bifurcation exponent , fixing which to is clearly an artifact of the algorithm, as suggested by [20, 115]. This did not turn out to be useful, neither by fixing values according to the observed distribution (see Figure 3) already in the optimization procedure nor as a postprocessing step. In both cases, we obtained very unnaturally small absolute radii throughout the generated vascular trees. In case of postprocessing, the leaf edge radii are changed to a wide range, conflicting with the idea behind CCO that leaf nodes should provide homogeneous supply.

The poor similarity of the angles between daughters (which are too small on average) is improved by shifting the bifurcation point in the direction of the angle bisector as a postprocessing step. A value of of the mean of the daughter edges was chosen for this purpose to fit the mean values of the generated vascular trees to the measured data.

Moreover, the poor similarity of the inclination angles is improved in a second postprocessing step. For this purpose, we consider one triangle per bifurcation, its vertices being the initial node of the parent and the terminal nodes of the daughter edges. We then compute the distance of the bifurcation point from this triangle relative to the longest edge of the triangle, also giving us a measure of nonflatness of the bifurcation for which the distribution can be computed separately for the clinical PV and HV data. As postprocessing, we modify the bifurcations in the generated vascular trees by moving the bifurcation points in normal direction to the triangles described before according to the nonflatness distribution, unless the triangle is degenerate. Note that this should not be done on a per-Strahler* order basis because data for the ranges we later want to apply the algorithm to is not available.

Calibration Results

The similarity results of the improved CCO output are listed in Table 6 next to the standard CCO results. Radii are not affected by the postprocessing. In particular, we do not rebalance radii throughout the tree afterwards; this turns out to further decrease similarity. Similarity in terms of length features becomes slightly worse by the postprocessing. Angle features, however, become much more similar, even the third angle for which an improvement was not aimed at by our postprocessing procedure.

The visual macroscopic difference introduced by the postprocessing is marginal; see Figure 4. Also notice that postprocessing the output of an optimization procedure decreases the optimality of the results. This cannot be seen in the cost function; actually the volume decreases by on average (PV and HV) due to shorter daughter segments at bifurcations, but at the same time the flow resistance increases. This indicates that future work should also include “functional” features rather than merely geometric ones.

As we could already observe for the standard CCO case, generated HV is less similar to measured ones than it is the case for PV.

Application of Improved CCO

Finally let us show an application of the improved construction algorithm. Starting from the dense PV example shown in Figure 2 (without pruning), we construct a tree supplying leaf nodes. A human liver has about to million lobules [114], so that in this setting one leaf node corresponds to about lobules. The resulting vascular tree is shown in Figure 5.

Figure 5: This vascular tree was generated applying the CCO algorithm including angle postprocessing to extend the clinical PV data set from Figure 4 without first pruning it. Each of the 10 000 leaf nodes in this setting is meant to supply about 125 lobules in the corresponding human liver.

5. Conclusions and Outlook

Even though the coarsest vascular structures in human livers have a large variability, the overall geometric properties are largely invariant between different individuals. Vascular trees generated by the CCO algorithm, that is, constructed to satisfy physical optimality conditions for homogeneous supply, show a certain similarity to measured vascular trees in reality. The similarity can be improved further by a postprocessing step modifying bifurcation angles.

Limitations

Several limitations of the methods presented here need to be mentioned.

The main limitation for obtaining more robust radius data is the imaging resolution, but more sophisticated image processing methods can possibly also improve the radius data. Higher resolution image data of small subsamples as described in Section 1 could be used for validation of geometric features at a smaller scale. However, a sufficient number of such datasets were not yet available for this study.

The similarity analysis and corresponding assessment of CCO are here limited to purely geometric features of the vascular structures. One could, for example, consider supplied territories by different parts of the vascular trees. In the model, these could be obtained, for example, by determining watersheds or Voronoï cells. Experimentally one could measure territories where an appropriately injected substance is distributed. However, such measurements in humans in vivo are not feasible. Functional properties as described in Section 2 could also be taken into account for validation of the results. Computing these, however, is not a trivial task and requires well-established and validated additional assumptions. Also, measuring appropriate data is difficult, in particular if in vivo measurements are desired.

We work with a simplified geometric representation of vascular trees allowing only straight cylindrical edges. This is only an approximation of real vascular edges which can be curved and of varying, not necessarily circular, cross-section. Varying out/inflow per leaf node (lobulus) can easily be integrated in the CCO implementation if such data is available. The flow model determining flow resistances does not yet take branching angles into account. The effective influence of the microscopic flow at bifurcations could be used to determine correction factors here.

Finally, the computational performance of our CCO implementation is not satisfactory yet for generating highly detailed vascular trees. This can possibly be remedied by using a multiscale approach partitioning the organ in separate supplied territories.

Outlook

Currently, only a single vascular network is generated algorithmically. If PV, HV, and possibly HA are considered simultaneously, penetration of the distinct vascular trees needs to be prevented. This can probably be achieved by considering a joint flow/pressure model also including an effective representation of the organ tissue between the vascular trees following [116] and by adding a nonpenetration constraint in the construction process.

The vascular tree data sets used here are obviously limited to human livers, but the similarity measures and construction algorithm are generic. We plan to apply the methods presented here to in vivo μCT scans of rats and mice and to μCT scans of murine hepatic corrosion casts.

Geometric representation of hepatic vascular trees determined by the methods presented here will be used for flow simulations in the context of multiscale liver modeling.

Acknowledgments

The authors acknowledge funding by the German Ministry of Education and Research (BMBF) via the systems biology network “Virtual Liver.” they also thank Jean Fasel and Horst Hahn for creating the corrosion casts and providing the corresponding vascular tree data sets. Furthermore they thank Uta Dahmen, Olaf Dirsch, Andrea Schenk, Stephan Zidowitz, Niko Komin, and Karolis Blaževičius for fruitful discussions on the topic.

References

M. A. Choti, “The Liver,” in The Johns Hopkins Atlas of Human Functional Anatomy, pp. 137–141, The Johns Hopkins University Press, 4th edition, 1997.

D. E. Malarkey, K. Johnson, L. Ryan, G. Boorman, and R. R. Maronpot, “New insights into functional aspects of liver morphology,” Toxicologic Pathology, vol. 33, no. 1, pp. 27–34, 2005.

C. Eipel, K. Abshagen, and B. Vollmar, “Regulation of hepatic blood flow: the hepatic arterial buffer response revisited,” World Journal of Gastroenterology, vol. 16, no. 48, pp. 6046–6057, 2010.

B. E. van Beers, I. Leconte, R. Materne, A. M. Smith, J. Jamart, and Y. Horsmans, “Hepatic perfusion parameters in chronic liver disease: dynamic ct measurements correlated with disease severity,” American Journal of Roentgenology, vol. 176, no. 3, pp. 667–673, 2001.

E. Leen, J. A. Goldberg, J. R. Anderson et al., “Hepatic perfusion changes in patients with liver metastases: comparison with those patients with cirrhosis,” Gut, vol. 34, no. 4, pp. 554–557, 1993.

H.-G. Holzhütter, D. Drasdo, T. Preusser, J. Lippert, and A. M. Henney, “The virtual liver: a multidisciplinary, multi-level challenge for systems biology,” Wiley Interdisciplinary Reviews, vol. 4, no. 3, pp. 221–235, 2012.

S. Willmann, J. Lippert, M. Sevestre, J. Solodenko, F. Fois, and W. Schmitt, “PK-Sim: a physiologically based pharmacokinetic “whole-body” model,” Drug Discovery Today, vol. 1, no. 4, pp. 121–124, 2003.

T. Eissing, L. Kuepfer, C. Becker, et al., “A computational systems biology software platform for multiscale modeling and simulation: integrating whole-body physiology, disease biology, and molecular reaction networks,” Frontiers in Physiology, vol. 2, no. 4, pp. 1–10, 2011.

R. Saxena, N. D. Theise, and J. M. Crawford, “Microanatomy of the human liver—exploring the hidden interfaces,” Hepatology, vol. 30, no. 6, pp. 1339–1346, 1999.

S. Hoehme, M. Brulport, A. Bauer et al., “Prediction and validation of cell alignment along microvessels as order principle to restore tissue architecture in liver regeneration,” Proceedings of the National Academy of Sciences of the United States of America, vol. 107, no. 23, pp. 10371–10376, 2010.

S. Hoehme and D. Drasdo, “A cell-based simulation software for multi-cellular systems,” Bioinformatics, vol. 26, no. 20, pp. 2641–2642, 2010.

I. Endo, H. Shimada, K. Takeda et al., “Successful duct-to-duct biliary reconstruction after right hemihepatectomy. Operative planning using virtual 3D reconstructed images,” Journal of Gastrointestinal Surgery, vol. 11, no. 5, pp. 666–670, 2007.

B. Hølund, H. Poulsen, and P. Schlichting, “Reproducibility of liver biopsy diagnosis in relation to the size of the specimen,” Scandinavian Journal of Gastroenterology, vol. 15, no. 3, pp. 329–335, 1980.

P. Bedossa, D. Dargère, and V. Paradis, “Sampling variability of liver fibrosis in chronic hepatitis C,” Hepatology, vol. 38, no. 6, pp. 1449–1457, 2003.

N. S. Goldstein, F. Hastah, M. V. Galan, and S. C. Gordon, “Fibrosis heterogeneity in nonalcoholic steatohepatitis and hepatitis C virus needle core biopsy specimens,” American Journal of Clinical Pathology, vol. 123, no. 3, pp. 382–387, 2005.

M. Karcaaltincaba and O. Akhan, “Imaging of hepatic steatosis and fatty sparing,” European Journal of Radiology, vol. 61, no. 1, pp. 33–43, 2007.

V. Campo-Ruiz, G. Y. Lauwers, R. R. Anderson, E. Delgado-Baeza, and S. González, “In vivo and ex vivo virtual biopsy of the liver with near-infrared, reflectance confocal microscopy,” Modern Pathology, vol. 18, no. 2, pp. 290–300, 2005.

S. Singhal, R. Henderson, K. Horsfield, K. Harding, and G. Cumming, “Morphometry of the human pulmonary arterial tree,” Circulation Research, vol. 33, no. 2, pp. 190–197, 1973.

W. Huang, R. T. Yen, M. McLaurine, and G. Bledsoe, “Morphometry of the human pulmonary vasculature,” Journal of Applied Physiology, vol. 81, no. 5, pp. 2123–2133, 1996.

S. Aharinejad, W. Schreiner, and F. Neumann, “Morphometry of human coronary arterial trees,” The Anatomical Record, vol. 251, no. 1, pp. 50–59, 1998.

N. Mittal, Y. Zhou, S. Ung, C. Linares, S. Molloi, and G. S. Kassab, “A computer reconstruction of the entire coronary arterial tree based on detailed morphometric data,” Annals of Biomedical Engineering, vol. 33, no. 8, pp. 1015–1026, 2005.

M. E. Martinez-Perez, A. D. Hughes, A. V. Stanton et al., “Retinal vascular tree morphology: a semi-automatic quantification,” IEEE Transactions on Biomedical Engineering, vol. 49, no. 8, pp. 912–917, 2002.

D. Selle, B. Preim, A. Schenk, and H.-O. Peitgen, “Analysis of vasculature for liver surgical planning,” IEEE Transactions on Medical Imaging, vol. 21, no. 11, pp. 1344–1357, 2002.

M. Zamir and S. Phipps, “Network analysis of an arterial tree,” Journal of Biomechanics, vol. 21, no. 1, pp. 25–34, 1988.

R. T. Yen, F. Y. Zhuang, and Y. C. Fung, “Morphometry of cat pulmonary venous tree,” Journal of Applied Physiology, vol. 55, no. 1, pp. 236–242, 1983.

G. S. Kassab, K. Imoto, F. C. White, C. A. Rider, Y.-C. Fung, and C. M. Bloor, “Coronary arterial tree remodeling in right ventricular hypertrophy,” American Journal of Physiology, vol. 265, no. 1, pp. H366–H375, 1993.

G. S. Kassab, D. H. Lin, and Y.-C. B. Fung, “Morphometry of pig coronary venous system,” American Journal of Physiology, vol. 267, no. 6, pp. H2100–H2113, 1994.

M. Gössl, M. Rosol, N. M. Malyar et al., “Functional anatomy and hemodynamic characteristics of vasa vasorum in the walls of porcine coronary arteries,” The Anatomical Record, vol. 272, no. 2, pp. 526–537, 2003.

S. Y. Wan, E. L. Ritman, and W. E. Higgins, “Multi-generational analysis and visualization of the vascular tree in 3D micro-CT images,” Computers in Biology and Medicine, vol. 32, no. 2, pp. 55–71, 2002.

D. A. Nordsletten, S. Blackett, M. D. Bentley, E. L. Ritman, and N. P. Smith, “Structural morphology of renal vasculature,” American Journal of Physiology, vol. 291, no. 1, pp. H296–H309, 2006.

R. J. Trudeau, Introduction to Graph Theory, Dover Publications, 1993.

C. Couinaud, Le Foie: Études Anatomiques et Chirurgicales, Masson, Paris, France, 1957.

K. Horsfield and G. Cumming, “Morphology of the bronchial tree in man,” Journal of Applied Physiology, vol. 24, no. 3, pp. 373–383, 1968.

G. S. Kassab, C. A. Rider, N. J. Tang, and Y.-C. B. Fung, “Morphometry of pig coronary arterial trees,” American Journal of Physiology, vol. 265, no. 1, pp. H350–H365, 1993.

S. B. Barker, G. Cumming, and K. Horsfield, “Quantitative morphometry of the branching structure of trees,” Journal of Theoretical Biology, vol. 40, no. 1, pp. 33–43, 1973.

K. Horsfield, “Morphometry of the small pulmonary arteries in man,” Circulation Research, vol. 42, no. 5, pp. 593–597, 1978.

R. E. Horton, “Erosional development of streams and their drainage basins; hydrophysical approaches to quantitative morphology,” Bulletin of the Geological Society of America, vol. 56, pp. 275–370, 1945.

A. N. Strahler, “Quantitative analysis of watershed geomorphology,” Transactions of the American Geophysical Union, no. 8, pp. 913–920, 1957.

Z. L. Jiang, G. S. Kassab, and Y. C. Fung, “Diameter-defined Strahler system and connectivity matrix of the pulmonary arterial tree,” Journal of Applied Physiology, vol. 76, no. 2, pp. 882–892, 1994.

R. Z. Gan, Y. Tian, R. T. Yen, and G. S. Kassab, “Morphometry of the dog pulmonary venous tree,” Journal of Applied Physiology, vol. 75, no. 1, pp. 432–440, 1993.

J. Drexl, V. Knappe, H. K. Hahn et al., “Accuracy analysis of vessel segmentation for a LITT dosimetry planning system, Perspective in Image-Guided Surgery,” in Proceedings of the Scientific Workshop on Medical Robotics, Navigation and Visualization, T. M. Buzug and T. C. Lueth, Eds., pp. 204–213, World Scientific, 2004.

E. VanBavel and J. A. E. Spaan, “Branching patterns in the porcine coronary arterial tree: estimation of flow heterogeneity,” Circulation Research, vol. 71, no. 5, pp. 1200–1212, 1992.

J. C. Parker, C. B. Cave, J. L. Ardell, C. R. Hamm, and S. G. Williams, “Vascular tree structure affects lung blood flow heterogeneity simulated in three dimensions,” Journal of Applied Physiology, vol. 83, no. 4, pp. 1370–1382, 1997.

J. N. Maina and P. van Gils, “Morphometric characterization of the airway and vascular systems of the lung of the domestic pig, Sus scrofa: comparison of the airway, arterial and venous systems,” Comparative Biochemistry and Physiology, vol. 130, no. 4, pp. 781–798, 2001.

H. B. M. Uylings, “Optimization of diameters and bifurcation angles in lung and vascular tree structures,” Bulletin of Mathematical Biology, vol. 39, no. 5, pp. 509–520, 1977.

M. F. Kiani and A. G. Hudetz, “Computer simulation of growth of anastomosing microvascular networks,” Journal of Theoretical Biology, vol. 150, no. 4, pp. 547–560, 1991.

R. Karch, F. Neumann, M. Neumann, P. Szawlowski, and W. Schreiner, “Voronoi polyhedra analysis of optimized arterial tree models,” Annals of Biomedical Engineering, vol. 31, no. 5, pp. 548–563, 2003.

M. Zamir, “On fractal properties of arterial trees,” Journal of Theoretical Biology, vol. 197, no. 4, pp. 517–526, 1999.

M. Zamir, “Nonsymmetrical bifurcations in arterial branching,” Journal of General Physiology, vol. 72, no. 6, pp. 837–845, 1978.

C. A. Dawson, G. S. Krenz, K. L. Karau, S. T. Haworth, C. C. Hanger, and J. H. Linehan, “Structure-function relationships in the pulmonary arterial tree,” Journal of Applied Physiology, vol. 86, no. 2, pp. 569–583, 1999.

C. D. Murray, “The physiological principle of minimum work. I. The vascular system and the cost of blood volume,” Proceedings of the National Academy of Sciences of the United States of America, vol. 12, no. 3, pp. 207–214, 1926.

K. Horsfield and M. J. Woldenberg, “Diameters and cross-sectional areas of branches in the human pulmonary arterial tree,” The Anatomical Record, vol. 223, no. 3, pp. 245–251, 1989.

R. Gödde and H. Kurz, “Structural and biophysical simulation of angiogenesis and vascular remodeling,” Developmental Dynamics, vol. 220, no. 4, pp. 387–401, 2001.

H. Kurz and K. Sandau, “Modelling of blood vessel development bifurcation pattern and hemodynamics, optimality and allometry,” Comments on Theoretical Biology, no. 4, pp. 261–291, 2007.

N. Suwa, T. Niwa, H. Fukasawa, and Y. Sasaki, “Estimation of intravascular blood pressure gradient by mathematical analysis of arterial casts,” The Tohoku Journal of Experimental Medicine, vol. 79, pp. 168–189, 1963.

T. Arts, R. T. Kruger, W. van Gerven, J. A. Lambregts, and R. S. Reneman, “Propagation velocity and reflection of pressure waves in the canine coronary artery,” The American Journal of Physiology, vol. 237, no. 4, pp. H469–H474, 1979.

G. S. Krenz, J. H. Linehan, and C. A. Dawson, “A fractal continuum model of the pulmonary arterial tree,” Journal of Applied Physiology, vol. 72, no. 6, pp. 2225–2237, 1992.

M. Marxen and R. M. Henkelman, “Branching tree model with fractal vascular resistance explains fractal perfusion heterogeneity,” American Journal of Physiology, vol. 284, no. 5, pp. H1848–H1857, 2003.

N. Mittal, Y. Zhou, C. Linares et al., “Analysis of blood flow in the entire coronary arterial tree,” American Journal of Physiology, vol. 289, no. 1, pp. H439–H446, 2005.

F. Y. Zhuang, Y. C. Fung, and R. T. Yen, “Analysis of blood flow in cat's lung with detailed anatomical and elasticity data,” Journal of Applied Physiology, vol. 55, no. 4, pp. 1341–1348, 1983.

W. Schreiner and P. F. Buxbaum, “Computer-optimization of vascular trees,” IEEE Transactions on Biomedical Engineering, vol. 40, no. 5, pp. 482–491, 1993.

R. Karch, F. Neumann, B. K. Podesser, M. Neumann, P. Szawlowski, and W. Schreiner, “Fractal properties of perfusion heterogeneity in optimized arterial trees: a model study,” Journal of General Physiology, vol. 122, no. 3, pp. 307–321, 2003.

A. Kamiya and T. Takahashi, “Quantitative assessments of morphological and functional properties of biological trees based on their fractal nature,” Journal of Applied Physiology, vol. 102, no. 6, pp. 2315–2323, 2007.

M. Zamir, “Optimality principles in arterial branching,” Journal of Theoretical Biology, vol. 62, no. 1, pp. 227–251, 1976.

M. Zamir, “The role of shear forces in arterial branching,” Journal of General Physiology, vol. 67, no. 2, pp. 213–222, 1976.

K. Horsfield and W. I. Gordon, “Morphometry of pulmonary veins in man,” Lung, vol. 159, no. 4, pp. 211–218, 1981.

W. Schreiner, F. Neumann, M. Neumann, A. End, S. M. Roedler, and S. Aharinejad, “The influence of optimization target selection on the structure of arterial tree models generated by constrained constructive optimization,” Journal of General Physiology, vol. 106, no. 4, pp. 583–599, 1995.

W. Schreiner, F. Neumann, M. Neumann, A. End, and M. R. Müller, “Structural quantification and bifurcation symmetry in arterial tree models generated by constrained constructive optimization,” Journal of Theoretical Biology, vol. 180, no. 2, pp. 161–174, 1996.

R. Karch, F. Neumann, M. Neumann, and W. Schreiner, “Staged growth of optimized arterial model trees,” Annals of Biomedical Engineering, vol. 28, no. 5, pp. 495–511, 2000.

J. H. G. M. van Beek, S. A. Roger, and J. B. Bassingthwaighte, “Regional myocardial flow heterogeneity explained with fractal networks,” American Journal of Physiology, vol. 26, pp. H1670–H1680, 1989.

F. Family, B. R. Masters, and D. E. Platt, “Fractal pattern formation in human retinal vessels,” Physica D, vol. 38, no. 1–3, pp. 98–103, 1989.

M. Sernetz, J. Wübbeke, and P. Wlczek, “Three-dimensional image analysis and fractal characterization of kidney arterial vessels,” Physica A, vol. 191, no. 1–4, pp. 13–16, 1992.

M. Zamir, “Fractal dimensions and multifractility in vascular branching,” Journal of Theoretical Biology, vol. 212, no. 2, pp. 183–190, 2001.

B. R. Masters, “Fractal analysis of the vascular tree in the human retina,” Annual Review of Biomedical Engineering, vol. 6, pp. 427–452, 2004.

D. Mancardi, G. Varetto, E. Bucci, F. Maniero, and C. Guiot, “Fractal parameters and vascular networks: facts & artifacts,” Theoretical Biology and Medical Modelling, vol. 5, p. 12, 2008.

R. Karch, F. Neumann, M. Neumann, and W. Schreiner, “A three-dimensional model for arterial tree representation, generated by constrained constructive optimization,” Computers in Biology and Medicine, vol. 29, no. 1, pp. 19–38, 1999.

W. Schreiner, “Concepts and features of arterial tree models generated by constrained constructive optimization,” Comments on Theoretical Biology, no. 1, pp. 103–136, 2001.

W. Schreiner, R. Karch, M. Neumann, F. Neumann, P. Szawlowski, and S. Roedler, “Optimized arterial trees supplying hollow organs,” Medical Engineering & Physics, vol. 28, no. 5, pp. 416–429, 2006.

H. Hahn, M. Georg, and H.-O. Peitgen, “Fractal aspects of three-dimensional vascular constructive optimization,” in Fractals in Biology and Medicine, Mathematics and Biosciences in Interaction, pp. 55–66, Springer, 2005.

M. Georg, T. Preusser, and H. K. Hahn, “Global constructive optimization of vascular systems,” Tech. Rep. 2010-11, Washington University in St. Louis, Department of Computer Science & Engineerin, 2010.

S. C. Gupta, Fluid Mechanics and Hydraulic Machines, Pearson Education, India, 2006.

M. E. Gottlieb, “Modelling blood vessels: a deterministic method with fractal structure based on physiological rules,” in Proceedings of the Annual International Conference of the IEEE Engineering in Medicine and Biology Society, vol. 12, 1990.

C. Y. Wang and J. B. Bassingthwaighte, “Area-filling distributive network model,” Mathematical and Computer Modelling, vol. 13, no. 10, pp. 27–33, 1990.

C. Y. Wang, J. B. Bassingthwaighte, and L. J. Weissman, “Bifurcating distributive system using Monte Carlo method,” Mathematical and Computer Modelling, vol. 16, no. 3, pp. 91–98, 1992.

T. R. Nelson and D. K. Manchester, “Modeling of lung morphogenesis using fractal geometries,” IEEE Transactions on Medical Imaging, vol. 7, no. 4, pp. 321–327, 1988.

H. K. Hahn, C. J. G. Evertsz, J. H. D. Fasel, and H.-O. Peitgen, “Fractal properties, segment anatomy, and interdependence of the human portal vein and the hepatic vein in 3D,” Fractals, vol. 11, no. 1, pp. 53–62, 2003.

F. Nekka, S. Kyriacos, C. Kerrigan, and L. Cartilier, “A model of growing vascular structures,” Bulletin of Mathematical Biology, vol. 58, no. 3, pp. 409–424, 1996.

D. Szczerba and G. Székely, “Macroscopic modeling of vascular systems,” in Medical Image Computing and Computer-Assisted Intervention—MICCAI 2002, T. Dohi and R. Kikinis, Eds., vol. 2489 of Lecture Notes in Computer Science, pp. 284–292, Springer, 2002.

M. Kretowski, Y. Rolland, J. Bézy-Wendling, and J.-L. Coatrieux, “Physiologically based modeling of 3-D vascular networks and CT Scan angiography,” IEEE Transactions on Medical Imaging, vol. 22, no. 2, pp. 248–257, 2003.

H. K. Hahn, Makroskopische Morphometrie und Dynamik der lebergefäße des Menschen [Diplom thesis], MeVis, Centrum für medizinische Diagnosesysteme und Visualisierung, 1999.

D. A. Darling, “The Kolmogorov-Smirnov, Cramérvon Mises tests,” The Annals of Mathematical Statistics, vol. 28, no. 4, pp. 823–838, 1957.

E. J. Dudewicz and S. N. Mishra, Modern Mathematical Statistics, Wiley Series in Probability and Mathematical Statistics, Wiley, 1988.

S. Stigler, “Fisher and the 5% level,” Chance, vol. 21, no. 4, p. 12, 2008.

J. Kiefer, “K-sample analogues of the Kolmogorov-Smirnov and Cramér-V. Mises tests,” The Annals of Mathematical Statistics, vol. 30, no. 2, pp. 420–447, 1959.

J. A. Koziol and D. P. Byar, “Percentage points of the asymptotic distributions of one and two sample K-S statistics for truncated or censored data,” Technometrics, vol. 17, no. 4, pp. 507–510, 1975.

M. A. Stephens, “EDF statistics for goodness of fit and some comparisons,” Journal of the American Statistical Association, vol. 69, no. 347, pp. 730–737, 1974.

A. N. Kolmogorov, “Sulla determinatzione empirica di una legge di distribuzione,” Giornale dell'Instituto Italiano degli Attuari, vol. 4, pp. 83–91, 1933.

R. Simard and P. L'Ecuyer, “Computing the two-sided Kolmogorov-Smirnov distribution,” Journal of Statistical Software, vol. 39, no. 11, pp. 1–18, 2011.

N. Smirnov, “Table for estimating the goodness of fit of empirical distributions,” The Annals of Mathematical Statistics, vol. 19, no. 2, pp. 279–281, 1948.

GNU octave http://www.octave.org/

J. W. Eaton, GNU Octave Manual, Network Theory Limited, 2002.

ROOT—an object oriented data analysis framework, http://root.cern.ch/

R. Brun and F. Rademakers, “ROOT—an object oriented data analysis framework,” in Proceedings of the AIHENP Workshop, vol. 389, pp. 81–86, Lausanne, Switzerland, September 1996.

R Development Core Team, R: A Language and Environment for Statistical Computing, 2009.

G. Marsaglia, W. W. Tsang, and J. Wang, “Evaluating kolmogorov's distribution,” Journal of Statistical Software, vol. 8, pp. 1–4, 2003.

F. J. Massey Jr., “The distribution of the maximum deviation between two sample cumulative step functions,” The Annals of Mathematical Statistics, vol. 22, no. 1, pp. 125–128, 1951.

M. Matsumoto and T. Nishimura, “Mersenne twister: a 623-dimensionally equidistributed uniform pseudo-random number generator,” ACM Transactions on Modeling and Computer Simulation, vol. 8, no. 1, pp. 3–30, 1998.

P. J. Kim, “On the exact and approximate sampling distribution of the two sample Kolmogorov-Smirnov criterion dmn, ,” Journal of the American Statistical Association, vol. 64, no. 328, pp. 1625–1637, 1969.

D. P. Bertsekas, Nonlinear Programming, Athena Scientific, Belmont, Mass, USA, 2nd edition, 1999.

W. Schreiner, F. Neumann, M. Neumann, A. End, and S. M. Roedler, “Anatomical variability and functional ability of vascular trees modeled by constrained constructive optimization,” Journal of Theoretical Biology, vol. 187, no. 2, pp. 147–158, 1997.

L. Armijo, “Minimization of functions having Lipschitz continuous first partial derivatives,” Pacific Journal of Mathematics, vol. 16, no. 1, pp. 1–3, 1966.

R. Fåhræus and T. Lindqvist, “The viscosity of the blood in narrow capillary tubes,” American Journal of Physiology, vol. 96, no. 3, pp. 562–568, 1931.

R. H. Haynes, “Physical basis of the dependence of blood viscosity on tube radius,” The American journal of physiology, vol. 198, pp. 1193–1200, 1960.

E. Kuntz and K. Hans-Dieter, Hepatology: Textbook and Atlas, Springer, Heidelberg, Germany, 3rd edition.

F. Neumann, W. Schreiner, and M. Neumann, “Computer simulation of coronary arterial trees,” Advances in Engineering Software, vol. 28, no. 6, pp. 353–357, 1997.

T. Ricken, U. Dahmen, and O. Dirsch, “A biphasic model for sinusoidal liver perfusion remodeling after outflow obstruction,” Biomechanics and Modeling in Mechanobiology, vol. 9, no. 4, pp. 435–450, 2010.