Fractal Properties of Perfusion Heterogeneity in Optimized Arterial Trees : A Model Study

It is now widely accepted that regional blood flows in organs and tissues, including the heart, the lung and skeletal muscles, are remarkably heterogeneous (Yipintsoi et al., 1973; Glenny and Robertson, 1990; Iversen and Nicolaysen, 1995; Sonntag et al., 1996). Regional flows per gram of uniformly sized tissue pieces have been documented to vary 6–10-fold with only small fluctuations over time (King and Bassingthwaighte, 1989), and recent experimental results suggest that perfusion heterogeneity is maintained even down to microscopic levels (Matsumoto et al., 1996; Glenny et al., 2000; Bauer et al., 2001; Kalliokoski et al., 2001). In particular, the magnitude of spatial heterogeneity, usually characterized by the relative dispersion (RD = SD/mean) of regional flows, is dependent on the scale of spatial resolution, i.e., on the size of the sample pieces. Bassingthwaighte et al. (1989) have empirically established a now well-known power-law relationship describing—within certain limits of piece sizes—the increase of RD with decreasing piece size,

where v is the volume of the tissue pieces used to calculate RD and v_ref is an arbitrarily chosen reference volume. On a double logarithmic scale Eq. 1 yields a linear relationship between the logarithm of RD and the logarithm of piece size with slope (1 − D), i.e., RD exhibits self-similarity upon scaling with respect to the size of sample pieces. In this sense, Eq. 1 represents a fractal relationship and the parameter D, which is a global measure of heterogeneity, is identified as a spatial fractal dimension (Mandelbrot, 1983; Bassingthwaighte et al., 1994). D = 1.0 indicates uniform flow and D = 1.5 indicates complete spatial randomness. For myocardial perfusion heterogeneity, an average D of ∼1.2 has been found (Bassingthwaighte et al., 1989), suggesting that heterogeneity is not random: Flows in adjacent tissue samples are positively correlated (i.e., they tend to be alike). By means of Eq. 1, the correlation coefficient r₁ between flows in nearest-neighbor pieces is directly expressible in terms of D (Van Beek et al., 1989),

i.e., for fractal heterogeneity Eq. 1 and 2 are equivalent and (as a consequence of self-similarity) nearest-neighbor correlations are independent of piece size.

There has been much debate whether the anatomic structures of the arterial trees or local metabolic needs of the supplied tissues are causative for the observed heterogeneities in regional flows (Van Beek, 1997; Deussen, 1998; Bassingthwaighte and Li, 1999; Balaban and Arai, 2001). In an attempt to reveal the contribution of vascular anatomy to regional flow distributions, numerous model studies have been reported, ranging from fractal self-similar branching models (Van Beek et al., 1989; Glenny and Robertson, 1991; Parker et al., 1997) and models derived from the statistics of morphometric data (VanBavel and Spaan, 1992), to whole-organ network reconstructions (Beard and Bassingthwaighte, 2000) based on the extensive set of anatomic data from VanBavel and Spaan (1992) and Kassab et al. (1993).

Here, we employ the method of constrained constructive optimization (CCO) to generate three-dimensional arterial model trees of some 10⁴ vessel segments (Schreiner, 1993; Karch et al., 1999, 2000a). Given the simplifying assumptions of this model, CCO trees exhibit a realistic visual appearance and reproduce various statistics of real arterial trees, such as branching angles (Schreiner et al., 1994) and diameter ratios of parent and daughter segments (Karch et al., 2000b). The method of CCO generates end-point positions compatible with the existing tree structure and allows—contrary to the approach of previous modeling efforts—to prescribe terminal flows as an independent boundary condition, thereby coming close to the notion that the primary causes of flow heterogeneity are the metabolic tissue needs and vascular anatomy is only a consequence thereof.

The values of D and r₁ observed in the left ventricular myocardium of animals under physiological conditions range from D = 1.21 (Bassingthwaighte et al., 1989) to D = 1.37 (Iversen and Nicolaysen, 1995), and from r₁ = 0.18 (Matsumoto et al., 1996) to r₁ = 0.92 (Mori et al., 1995), respectively. However, Kleen et al. (1997) reported D = 1.6 for the subepicardial myocardium of pigs at baseline and after the induction of coronary stenosis, and Mori et al. (1995) observed r₁ = −0.16 for aggregated adjacent sample regions in a dog myocardium under reduced coronary perfusion pressure. In this paper, we study spatial dispersions and nearest-neighbor correlations of regional flows in model trees generated by two different sets of boundary conditions for the terminal flows: (a) equal terminal flows (model 1) and (b) terminal flows set proportional to the volumes of Voronoi polyhedra associated with terminal segment end-points (model 2).

The method of CCO for arterial tree representation is derived from the principle of optimal design (Rosen, 1967), which has long been hypothesized for single arterial bifurcations (Murray, 1926a,b; Zamir, 1976a,b). Starting from a root segment, CCO trees are generated by adding new terminal segments that connect randomly chosen points of a geometric model of the perfusion area with segments of the tree grown thus far, until a given number N_term of terminal segments is reached. At each step of growth, the geometry of the newly created bifurcation is optimized under a prescribed set of physiologic boundary conditions. However, no attempt is made to model the growth process of real arterial trees (Hudlicka et al., 1986).

For ease of reference, we restate the main points of CCO: Arterial trees are modeled as binary branching trees and segments are assumed rigid cylindrical tubes. Arterial bifurcations (Zamir et al., 1983; Zamir, 2001) are characterized by the diameters and lengths of the parent segment and its two daughter branches—the specific shape of the vascular wall of bifurcations and its impact on hemodynamics are not considered. Likewise, blood is assumed an incompressible, homogeneous Newtonian fluid at steady-state and laminar flow conditions, and Poiseuille's law is used to estimate the hydrodynamic resistance R of individual segments,

where l and d are segment length and inner diameter, and μ is the viscosity of blood. The position of each new bifurcation is calculated so as to minimize total intravascular volume (Kamiya and Togawa, 1972; Schreiner et al., 1995) under the following boundary conditions: At a given perfusion pressure p_perf, perfusion flow Q_perf, terminal flows Q_term, and at a preset diameter of the root segment, segment diameters are scaled to obey the relationship (Sherman, 1981)

at each bifurcation, where d₀, d₁, and d₂ are the diameters of the parent and daughter segments, respectively, and the constant γ > 0 characterizes the reduction of segment diameters across bifurcations.

After model trees are grown to their final size (i.e., to the prescribed number N_term of terminal segments) we consider the end-points of these terminals. These points represent locations of blood supply for the surrounding tissue regions, and in model 2 we hypothesize that these regions can be approximated by Voronoi polyhedra: Given a set of distinct points in space, the Voronoi polyhedron associated with one of these points is defined as the convex region of space closer to this point than to any other point of the given set (Preparata and Shamos, 1985). The Voronoi tessellation of an entire point set yields a partition of space into convex nonoverlapping polyhedra that completely fill the space (Fig. 1). Voronoi tessellations originally emerged in pure mathematics as a number-theoretical problem (Dirichlet, 1850; Voronoi, 1908) and have since then been successfully applied in biology and medicine to model cellular patterns (Honda, 1978), capillary domains (Hoofd et al., 1985), and to describe the apparent end-point distribution of arteries in the chorioallantoic membrane of chicken eggs (Kurz and Sandau, 1997). To construct the Voronoi tessellation of terminal segment end-points, we have used the Qhull software package (designed by C.B. Barber and H. Huhdanpaa, University of Minnesota, MN). This package implements the Quickhull algorithm (Preparata and Shamos, 1985), which reduces the construction of a n-dimensional Voronoi diagram to the problem of finding a convex hull (i.e., the smallest convex set that contains a given set of points) in (n + 1) dimensions. Qhull provides vertex coordinates of the Voronoi polyhedra faces for all pairs of adjacent input sites and thus allows straight-forward calculation of polyhedra volumes V: From the areas A_k of individual faces (Ruocco et al., 1992),

where n_v^(k) is the number of vertices of the k-th face, r_j^(k) is the position vector of the j-th vertex of this face, and × denotes the cross product of the respective vectors, we have

where h_k is the normal distance of the polyhedron-center from its k-th face and N_f is the total number of polyhedra faces.

In model 1 we assume that all terminal segments carry the same fraction of the total perfusion flow, Q_term = Q_perf/N_term, so that the variability of local flows is solely determined by the terminals' spatial position, i.e., by the variability of the number of terminals within the respective sample pieces.

In model 2 we assume that segments with relatively high flows also perfuse large volumes of tissue (VanBavel and Spaan, 1992) and we model the tissue regions supplied by individual terminal segments by their associated Voronoi polyhedra and set terminal flows proportional to the corresponding polyhedra volumes. Note that the method of CCO assigns terminal flows after each step of growth and rescales segment diameters during optimization to fulfill these new boundary conditions together with Eq. 4. For model 2 this algorithm is computationally demanding, since each addition of a new terminal requires a complete Voronoi tessellation of all terminal segment end-points. However, CCO trees generated in this way are not significantly different (with respect to, for example, the distribution of segment diameters, pressure profile, and total intravascular volume) from trees in which Voronoi tessellation with terminal flow assignment and rescaling is performed only once after the trees are grown to their final size (Karch et al., 2003). Therefore, the trees of model 2 are generated using the latter approach and actually represent rescaled versions of the respective trees of model 1.

In fractal dispersion analysis we are interested in the scaling behavior of the relative dispersion of regional flow as a function of spatial resolution: The perfusion domain is successively subdivided into 8 = 2 × 2 × 2, 27 = 3 × 3 × 3, … , 4096 = 16 × 16 × 16 cubic volume elements of equal size and within each individual element the regional flow q is calculated as the sum of flow through all terminal segments within the respective element, divided by the element volume. For each resolution (i.e., element volume v), we calculate the overall mean <q> (which is constant and equal to the total perfusion flow Q_perf divided by the volume V_perf of the perfusion area), the standard deviation SD(v), and relative dispersion RD(v) = SD(v)/<q> of regional flows and obtain RD as a function of element size v. On a double logarithmic scale, Eq. 1 yields a linear relation between RD and resolution,

and the spatial fractal dimension D is given by the slope (1 − D) of this linear relation over a range of resolutions, where the slope is constant (Glenny et al., 1991).

Fractal dimensions observed for regional flows in the lung (Glenny and Robertson, 1990) and the myocardium (Bassingthwaighte et al., 1989) demonstrate that flow heterogeneity is not random—flows in adjacent tissue samples are spatially correlated. Glenny (1992) proposed Pearson's linear correlation coefficient,

as a measure of the linear relationship between nearest-neighbor regional flows. The variables X and Y are paired observations (x_i, y_i) of the magnitudes of flows in adjacent volume elements, <X> and <Y> are their respective means, and n is the number of paired observations. Note that we assume that the tissue is isotropic and the average in Eq. 8 can be performed over all 3 spatial directions. Thus, we effectively conduct a one-dimensional analysis. By virtue of Eq. 2, determining r₁ provides an additional method to estimate D,

i.e., the global measure of heterogeneity, D, is related to the measure of averaged local flow relationships, r₁. Given, the flow distribution obeys a fractal relationship, Eq. 1, D is independent of the resolution used to measure r₁, i.e., the correlation between nearest neighbors at one level of resolution is the same as between neighbors at a higher (or lower) level (Bassingthwaighte and Beyer, 1991).

The fractal branching tree (FBT) model, originally introduced by Van Beek et al. (1989), represents a dichotomously branching tree in which, at each bifurcation, a fraction α (flow asymmetry parameter, 0 ≤ α ≤ 0.5) of the parent blood flow enters one daughter segment and a fraction (1 − α) enters the other. After n bifurcations, i.e., in the n-th generation, there are 2ⁿ terminal segments carrying discrete flows of magnitude

each occurring with a relative frequency

where Q₀ is the flow into the root segment and k = 0, 1, … , n. Thus, this model describes flow dispersion without considering the spatial organization of the arterial tree.

Statistical analysis of model results was performed using the SAS system version 8.1 (SAS Institute). Group data are shown as mean ± SD. Normal distribution of data was tested with the Shapiro Wilks test. All correlation values refer to Pearson's correlation coefficient. Differences in D-values were tested by means of Student's t test for paired data. The level of statistical significance was set to 5%.

We have computed five realizations of model 1 trees by applying different seeds for the random number sequence used to generate new terminal locations. Simulation parameters were set so as to model the large arteries of the coronary arterial tree supplying ∼100 g of myocardial tissue (having a volume of ≈100 cm³) under cardiac arrest and maximum vasodilation (Table I). Trees were grown up to a size of N_term = 16,000 terminal segments in a cube of 4.642 cm ≈ (100 cm³)^1/3 side-length. N_term was chosen to yield RD-values for model 1 similar to experimental data (Glenny and Robertson, 1995). The trees of model 2 were obtained from the previously generated model 1 trees by constructing the Voronoi polyhedra of terminal segment end-points, calculating the respective polyhedra volumes, and setting terminal flows proportional to these volumes. Note that Voronoi tessellations were performed in a way preserving the boundaries of the perfusion area (i.e., original area boundaries are exactly recovered as parts of the polyhedra, see Fig. 1). In a final step, these trees were rescaled, i.e., segment diameters recalculated, so as to meet the modified boundary conditions. Therefore, the trees of model 2 represent rescaled versions of model 1 trees and have identical topological structures. Fig. 2 A shows a specific realization of a model 1 tree and Fig. 2 B illustrates polyhedra resulting from the Voronoi decomposition of end-points of a tree grown to N_term = 400 terminal segments.

Fig. 3 shows the probability density functions of relative regional flows, q* = q/<q>, i.e., regional flows q divided by the mean regional flow <q>, for four specific volume element sizes v. The distributions of regional flows are fairly symmetric about q* = 1 for each element size (ranging from v = 12.5 ml with n = 8 pieces to v = 0.024 ml with n = 4,096 pieces) and each realization of models 1 and 2, so that the RD is an adequate measure of their spread (Bassingthwaighte et al., 1994). With increasing resolution (i.e., decreasing element size) the widths of the probability density functions grow and their relative dispersions RD increase. At each specific level of resolution, RD of model 2 is smaller than the respective value of model 1 (see Table II), indicating that “Voronoi-scaled” terminal flows (model 2) render regional flows more homogeneous than equal terminal flows (model 1). This result is illustrated in Fig. 4, where the clustering of high- and low-flow regions is less pronounced for model 2 (Fig. 4 B) than for model 1 (Fig. 4 A): model 2 looks more random (i.e., less heterogeneous) than model 1.

Fig. 5 presents the relative dispersion RD(v) of regional flow as a function of volume element size v in a double-logarithmic plot. The fractal dimension D was obtained from the slope (1 − D) in a linear regression of log RD vs. log v for each individual realization of models 1 and 2 (see Eq. 7). The three largest volume decompositions were excluded from the linear fit, because (a) the samples were small (n = 8, n = 27, and n = 64, respectively), and (b) log RD values did not fit the fractal relationship of Eq. 7. Likewise, the three smallest decompositions were also excluded, because of the small mean number of terminals per sample volume (16,000/4,096 ≈ 3.9 for the smallest sample volume v = 0.024 ml). The linear correlation between the remaining paired values of log RD and log v was high, r² = 0.99 ± 0.005 for model 1 and r² = 0.99 ± 0.007 for model 2. The fractal dimension D for models 1 and 2 was estimated as the mean of D-values obtained from individual linear least-squares fits (see Table II). We find D = 1.405 ± 0.009 for model 1 and D = 1.691 ± 0.016 for model 2, both being significantly different from D = 1.5 for randomly distributed terminal flows (Student's t test, P < 0.0001, n = 5).

Fig. 6 shows the linear correlation coefficient r₁, Eq. 8, between regional flows in adjacent sample pieces as a function of sample volume size v for the trees of models 1 and 2. Fig. 6 illustrates two distinctive properties of the models: (a) Although model 1 always exhibits positive correlations, we observe significant negative correlations throughout model 2. (b) In both models, r₁ is roughly independent of the spatial resolution v. For the nearest-neighbor correlations, averaged over the same range of piece sizes used in the previous section to estimate D, we find <r₁> = 0.122 ± 0.025 for model 1 and <r₁> = −0.175 ± 0.027 for model 2 (see Table III). Since r₁ proved roughly constant over the range of piece sizes v where log RD scales linearly with log v, Eq. 9 can be applied for an additional estimate of the fractal dimension D from nearest-neighbor correlations r₁. From the D-values calculated from r₁ at individual resolutions v, we find D = 1.417 ± 0.016 for model 1 and D = 1.639 ± 0.024 for model 2, respectively (Table III). These values are in good agreement with the results we obtained for D from the method of scaling-analysis of log RD vs. log v (Table II).

In the following, we report some consequences of the different terminal flow settings of models 1 and 2 on the distribution of flow within the whole tree and compare our results with the predictions of the FBT model. Fig. 7 A shows the distribution of relative segmental flows for model 1, model 2, and the FBT model with n = 14 generations, i.e., 2¹⁴ = 16,384 terminal segments (yielding roughly the same number of terminals as the CCO models with N_term = 16,000), and a flow-asymmetry parameter α = 0.46, obtained by a best-fit procedure to the flow distribution of model 2. While the probability density function (pdf) for segmental flows of model 1 closely follows a power law, p(x) ∼ x^−β, with an exponent β = 1.83 ± 0.03, the pdfs of both model 2 and the FBT model exhibit a distinctive maximum and decrease at their lower ends, but otherwise behave very similarly to model 1. Fig. 7 B presents the pdfs of relative terminal flows for model 2 and the FBT model. Both models are well approximated by a lognormal distribution,

where <x>_0.5 is the median and σ_g is the geometric standard deviation. The differences in the respective parameters between the models are very small, <x>_0.5 = (5.93 ± 0.04) × 10⁻⁵, σ_g = 1.39 ± 0.01 (model 2) and <x>_0.5 = (5.74 ± 0.04) × 10⁻⁵, σ_g = 1.34 ± 0.01 (FBT).

During the last decade since its discovery (Bassingthwaighte et al., 1989), a wealth of model studies have been performed to broaden the understanding of the fractal character of spatial heterogeneity of regional blood flow in mammalian organs. Van Beek et al. (1989) modeled the myocardial vascular tree as a dichotomous branching network and concluded that a relatively small deviation from a symmetric flow splitting at bifurcations was sufficient to explain the observed heterogeneity. Glenny and Robertson (1991) applied this model to study pulmonary blood flow heterogeneity and reported an excellent agreement with experimental data. VanBavel and Spaan (1992) constructed computer models of the porcine coronary arterial tree from the statistics of their morphometric data and concluded that the coronary branching pattern was a major determinant of flow heterogeneity. Although these models provided valuable insight on flow heterogeneity, they did not represent anatomic features accurately (Bassingthwaighte et al., 1989), such as the spatial organization of the arterial tree. Glenny and Robertson (1995) reported on a space-filling three-dimensional model with branches along the three orthogonal directions and randomly chosen flow asymmetry at each bifurcation and found that this model explained both heterogeneity and spatial correlation of regional pulmonary perfusion. Parker et al. (1997) extended this model to incorporate a range of branching angles and daughter-parent length ratios and concluded that a flow-symmetric model was sufficient to predict the heterogeneity observed in dog lungs, while negative correlations of flows with distance were obtained only with an asymmetric flow-ratio along bifurcations. Based on the extensive set of anatomic data from VanBavel and Spaan (1992) and Kassab et al. (1993), Beard and Bassingthwaighte (2000) presented a whole-organ model of the coronary arterial network and successfully predicted the degree of heterogeneity and the spatial correlations of regional flows observed in the myocardium of animals. Qian and Bassingthwaighte (2000) presented a general flow bifurcation model based on a dichotomous branching tree with asymmetric branching and random flow variation at each bifurcation and showed that this model gives rise to an asymptotic lognormal flow distribution and fractal scaling in the dispersion of regional flows consistent with experimental data. Kendal (2001) applied the statistical theory of exponential dispersion models to regional organ blood flow and recovered Bassingthwaighte's fractal scaling relation, Eq. 1, assuming that the dispersion model was additive, scale invariant, and represented a compound Poisson distribution. Recently, Marxen and Henkelman (2003) explained fractal perfusion heterogeneity in a branching tree model as a consequence of scale-independent branching asymmetry and fractal vascular resistance.

Fig. 5 shows experimental data of myocardial perfusion heterogeneity (log RD vs. log v) of various animals as reported by four independent research groups, together with simulation results of models 1 and 2. We have converted the experimental values of myocardial tissue mass to sample element volumes v assuming a specific mass of 1 g/ml for heart tissue (Spaan, 1991). Results for models 1 and 2, together with experimental data, are summarized in Table IV.

The value of D = 1.41 ± 0.01 for N_term = 16,000 predicted by model 1 (n = 5 simulations) is in good agreement with the experimental findings of Iversen and Nicolaysen (1995), D = 1.37 ± 0.06 for n = 46 rabbits, and of Kleen et al. (1997), D = 1.34 ± 0.04 (D = 1.39 ± 0.06) for n = 7 (n = 6) pigs, and only moderately above the range of individual D-values (1.18–1.37) reported for n = 6 rabbits by Bassingthwaighte et al. (1989). However, model 1 overestimates the average D-values reported for n = 11 sheep (D = 1.17 ± 0.07), n = 10 baboons (D = 1.21 ± 0.04), and n = 6 rabbits (D = 1.25 ± 0.07) in the original study on the fractal nature of regional myocardial blood flow heterogeneity by Bassingthwaighte et al. (1989), as well as the average D of 1.21 ± 0.08 for n = 4 dogs measured by Mori et al. (1995). As detailed below, the fractal dimension D of model 1 decreases with increasing N_term (for N_term = 32,000 we obtain D = 1.38 ± 0.01).

Model 2 predicts a fractal dimension of D = 1.69 ± 0.02 (n = 5 realizations), independent of N_term. To the best of our knowledge, a value of D > 1.5 has not been reported in previous model studies, although—as has been pointed out by Bassingthwaighte et al. (1990)—it is possible to have a fractal dimension D > 1.5 when inverse correlations between nearest neighbors become evident. Yet, in a recent experimental study, Kleen et al. (1997) have observed D = 1.64 ± 0.04 (1.66 ± 0.03) and D = 1.61 ± 0.03 (1.65 ± 0.04) in the epicardial layers of n = 6 (n = 7) pig hearts at baseline and after the induction of coronary stenosis.

Since the absolute values of perfusion heterogeneity, i.e., the relative dispersions RD of regional flows, also depend on the number N_term of end-points (see below), N_term was chosen so as to yield physiologically reasonable RD values for model 1 (Fig. 5). With the same N_term, model 2 predicts realistic RD values only for sample volumes v < 1 ml. In this range, model 2 is in reasonable agreement with the data of Mori et al. (1995). However, model 2 underestimates RD for large volume elements when compared with the experimental data of Bassingthwaighte et al. (1989) for sheep, rabbits, and baboons.

As to the spatial correlation of blood flow to neighboring sample volume elements, model 1 predicts a mean linear correlation coefficient of <r₁> = 0.12 ± 0.03 (n = 9) over the size range of volume elements considered (Table III and Fig. 6). Kleen et al. (1997) have reported r₁ = 0.26 ± 0.07 (n = 6) and r₁ = 0.26 ± 0.04 (n = 7) for pig myocardial tissue samples of 316 mg ≈ 0.316 ml size. For this specific sample element size, model 1 predicts r₁ = 0.14 ± 0.01 (n = 5), thus underestimating the experimental values. Since in model 1 D depends on N_term, we observe the respective variation also for r₁ (see Eq. 2). For N_term = 32,000, model 1 gives r₁ = 0.22 ± 0.02 (n = 5) for the above-mentioned tissue sample size, in good agreement with the experimental values. Moreover, the r₁ values predicted by model 1 are also within the range of 0.18–0.32 reported by Matsumoto et al. (1996) for the within-layer correlation of adjacent regional flows in rabbit hearts. However, the r₁ values of 0.42–0.92 measured by Mori et al. (1995) in dog hearts are underestimated by model 1.

Model 2 gives a negative linear correlation coefficient of <r₁> = −0.18 ± 0.03 (n = 9), which has not been observed in experimental studies under normal physiological conditions, and thus fails to predict regional flow heterogeneity of the normal heart. However, Mori et al. (1995) reported r₁ = −0.16 for sample aggregates of ∼40 mg ≈ 0.040 ml in a dog heart under reduced coronary perfusion pressure. These results suggest the hypothesis that negative nearest-neighbor correlations as predicted by model 2 may indicate myocardial ischemia. This view is consistent with the tendency of nearest-neighbor correlations being attenuated in the myocardium of dogs under reduced coronary perfusion pressure, as observed by Mori et al. (1995). Kleen et al. (1997) reported D > 1.5 for epicardial layers in pig hearts, yet their corresponding r₁ values of 0.21 ± 0.06 (n = 6) and 0.34 ± 0.05 (n = 7) failed to follow the fractal relation, Eq. 2, which would predict negative r₁ values for D > 1.5. These authors interpreted their D values >1.5 as indicating very heterogeneous, nonself-similar perfusion in the subepicardial myocardium. Glenny (1992) observed negative correlations of regional pulmonary perfusion for distant sample elements (but not for adjacent pieces) and attributed these findings to the effect of positive correlations at small distances in a finite volume under the boundary condition of a fixed total perfusion flow.

It is a specific feature of the method of CCO to prescribe terminal flows during growth and optimization. As such, this method allows to model (adaptive) effects of local substrate needs of the supplied tissue on the distribution of resistance and flow within the model trees. Note that model 1 does not define a specific supply volume (neither size, nor shape) associated with terminals. Only model 2 defines both size and shape of terminal supply regions as well as associated terminal flows.

The influence of different terminal flow distributions quickly diminishes toward the root (Fig. 7 A) through the averaging effect of the resistance pathways within the trees. Since the smallest flows can occur only at the periphery, the different pdfs of models 1 and 2 in Fig. 7 A mainly originate from the terminals.

CCO trees with terminal flows set proportional to Voronoi polyhedra volumes and FBT models give rise to approximate lognormal terminal flow distributions and the flow-asymmetry parameter α can be chosen to render these distributions very similar (Fig. 7 B), a finding not immediately obvious, given the different mechanisms that determine the respective flows. This result is further supported by a recent paper of Qian and Bassingthwaighte (2000), who showed that dichotomous branching tree models with random flow variations at each bifurcation yield an asymptotic lognormal flow distribution.

Finally, we note that Kendal (2001) has shown that Bassingthwaighte's fractal scaling relation, Eq. 1, is a consequence of a compound Poisson-gamma distribution, representing regional flow. This author assumed that the sites of microsphere entrapment (corresponding to the end-points of our model trees) are randomly (Poisson) distributed within the sample pieces and that blood flow at the sites of entrapment obeys a gamma distribution,

where Γ(a) denotes the gamma function and the mean and variance of p(x) are given by ab and ab², respectively. Above we have used a lognormal distribution to characterize terminal flows; however, Vaz and Fortes (1988) have shown that lognormal and two-parameter gamma distributions with suitable parameters are very similar. This is illustrated in Fig. 7 B, displaying a least-squares fit of a two-parameter gamma distribution to the relative terminal flows of model 2. We find a = 10.1 ± 0.3 and b = (6.0 ± 0.2) × 10⁻⁶.

The method of CCO does not directly incorporate anatomic information to reconstruct arterial trees. Yet, the model contains adjustable parameters (see Table I) that determine various properties of the simulated trees.

The number N_term of terminal segments, i.e., the trees' spatial resolution, is of particular relevance in predicting the dispersion of regional flows. Fig. 8 shows the dependence of the relative dispersion RD and fractal dimension D on N_term for N_term = 8,000, 16,000, and 32,000 terminal segments. The absolute values of RD systematically decrease with increasing N_term, both for models 1 and 2. For model 1, where RD measures the variability of the number of end-points in the sample volume elements, this behavior can be qualitatively understood by the properties of N randomly positioned points (drawn from a uniform distribution) in a spatial domain of volume V. For such a distribution, the probability P(n_v = k) of finding exactly n_v = k points in a cubic sample box of volume v is given by a Poisson distribution,

where λ = N/V is the spatial point density (e.g., Stoyan and Stoyan, 1994). Since for a Poisson distribution the mean and variance are <n_v> = Var(n_v) = λv, we have

i.e., RD decreases with increasing N as

⁠. Note that the last equation directly gives

with D = 1.5, in accordance with the result of King et al. (1990) obtained by direct simulation of randomly distributed points. RD values of model 2 also show a

dependence, since respective RD values collapse onto a single line when scaled by

\(\sqrt{N}\)

⁠, see Fig. 8 B. Although the fractal dimension D of model 2 is practically insensitive to the choice of N_term, model 1 exhibits a gradual decrease of D with increasing N_term from D = 1.45 ± 0.01 at N_term = 8,000 to D = 1.38 ± 0.01 at N_term = 32,000. At small N_term, D is close to the value of 1.5 for randomly distributed points, while for larger N_term the spatial distribution of end-points is more and more disturbed by the tree structure, resulting in a decreasing D. Our choice of N_term = 16,000 has been dictated mainly by adjusting model 1 to realistic RD values and by available computational resources, in particular with respect to the memory-demanding calculation of Voronoi-tessellations of trees grown to their final size. However, VanBavel and Spaan (1992) have estimated the total number of end segments (diameters between 5 and 10 μm) in myocardial tissue to be ∼5.5 × 10⁶/100 g, i.e., more than two orders of magnitude larger than in the present models. By extrapolation of D to a realistic number of end-points we estimate D ≈ 1.2 at N_term = 10⁶, a value still consistent with experimental data (Bassingthwaighte et al., 1989). We finally note that, by virtue of Eq. 2, a dependence of D on N_term entails a respective variation of r₁. Accordingly, we find for model 1 <r₁> = 0.09 ± 0.02 (N_term = 8,000), <r₁> = 0.12 ± 0.03 (N_term = 16,000), and <r₁> = 0.18 ± 0.05 (N_term = 32,000).

The bifurcation exponent γ, Eq. 4, defines the amount by which segment diameters decrease along arterial bifurcations when moving from the root toward the periphery. As such, γ influences—besides segment length and blood viscosity—the distribution of segmental resistance and flow within the trees, see Eq. 3. As a consequence, different values of γ cause model trees to follow different paths during growth, resulting in trees of different structure and spatial arrangement of terminals (Schreiner, 2001). Thus, we expect γ to be an important determinant of the amount of heterogeneity predicted by CCO trees. For γ = 2.0 we find D = 1.27 ± 0.01 (n = 5) for model 1 and D = 1.60 ± 0.02 (n = 5) for model 2. γ = 3.0 tends to increase D, yielding D = 1.46 ± 0.01 (n = 5) for model 1 and D = 1.69 ± 0.02 (n = 5) for model 2.

The representation of the coronary arterial tree by the method of CCO is an overly simple description of reality. In the framework of this model, heterogeneity of blood flow per unit volume of tissue depends both on the positioning and on the amount of flow delivered by individual terminal segments (for model 1, only the positions are relevant, since all terminal flows are set to the same value). At each step of tree construction, terminals are positioned such as to be compatible with the tree generated thus far, but are otherwise randomly arranged in space, i.e., drawn from pseudo-random numbers with a uniform distribution. As such, the model represents a purely geometric and mechanistic approach, neither taking into account a realistic spatial distribution of arteriolar end-points nor the complex interactions of regulatory processes that have been hypothesized as a possible cause for the observed flow heterogeneities (Balaban and Arai, 2001).

In model 2 we essentially make two assumptions: (a) terminal flow is proportional to the volume of supplied tissue, and (b) the regions of tissue supplied by individual terminals can be modeled by their respective Voronoi polyhedra. While the first assumption seems reasonable (VanBavel and Spaan, 1992), the second is a hypothesis that has to be validated by experimental data. Seiler et al. (1992) have reported a power-law relationship between luminal cross-sectional areas of human coronary artery segments and perfused regional myocardial mass with exponents in the range between 0.62 and 0.82. In good agreement with these results, we find for model 2 a ∼ v_VP^β with β = 0.77 ± 0.01, r² > 0.95 (n = 5), where a is the luminal cross-sectional area of segments and v_VP is the volume of tissue supplied by the respective distal subtree, Fig. 9. Here, tissue volume v_VP is the total volume of the Voronoi polyhedra attached to the subtrees' terminals. Moreover, model 2 yields a scaling-relation of intravascular volume as a function of the volume of tissue supplied consistent with experimental (Prothero, 1980) and theoretical (West et al., 1997) findings (Karch et al., 2003). Yet, setting terminal flows proportional to polyhedra volumes homogenizes regional flows to an extent that reduces RD to unphysiological values, in particular for large sample pieces (Fig. 5), and results in negative nearest-neighbor correlations of regional flows in contrast to experimental data of the normal heart.

The structure of the model trees, which in turn guides the spatial positioning of terminals, is primarily based on an optimality principle. However, for complex biological systems any such principle is only a hypothesis (Kassab and Fung, 1995) and minimizing intravascular volume might not be adequate for all parts of the cardiovascular system with different functions and modes of operation (Zamir, 1977). Moreover, as has been pointed out by Pries et al. (1995), optimization of vascular beds is possible only within the limits imposed by functional requirements. The usage of a single optimization principle might be the reason why the model fails to predict a connective structure that agrees with experimental data at all levels of a topological ordering parameter, such as the diameter-defined Strahler system (Kassab et al., 1993). Thus, it would be desirable to model the topological structure directly from Kassab's connectivity matrix, similar to the approach taken by Beard and Bassingthwaighte (2000).

CCO models assume that blood flow in arterial trees is steady and that the hydrodynamic resistance of each segment follows Poiseuille's law, Eq. 3. Given the pulsatile nature of flow and the complex geometry of arterial walls, in particular in the neighborhood of vessel junctions, these assumptions are an idealization. Actually, blood flow depends on the dimensionless Reynolds and Womersley numbers (e.g., Fung, 1990): only if both are much smaller than 1, deviations from Poiseuille's law can be neglected. Kassab et al. (1997) have pointed out that for the coronary arteries, Womersley numbers are <1 and flow may therefore be considered approximately quasisteady. Yet, application of Poiseuille's law in the larger coronary arteries underestimates flow resistance due to arterial branching and flow entrance and exit disturbances. Not taking into account the pulsatile nature of blood flow means that regional flows predicted by CCO models represent values averaged over many heart beats. On the other hand, experimental data demonstrate that temporal fluctuations in regional myocardial flows are small in comparison with spatial fluctuations (King and Bassingthwaighte, 1989) and that spatial blood flow heterogeneity is temporally stable over hours (Deussen et al., 1996). These results suggest that neglecting pulsatile flow in CCO trees may be of minor importance within the context of modeling spatial blood flow heterogeneity.

Finally, we note that Wang and Bassingthwaighte (2001) showed for a parallel arrangement of capillaries that capillary supply regions, i.e., regions based on the functional measure of capillary diffusion, cannot—except in very special cases—be approximated by Voronoi polygons. However, terminal segments of CCO models do not represent real terminal arterioles or capillaries, nor do they form a system of parallel pipes. In particular, Ellsworth et al. (1994) have pointed out that capillaries are not the only source of oxygen, and that a complex pattern of oxygen exchange may exist between arterioles, venules, and adjacent capillary networks. Yet, although model 2 yields reasonable distributions of terminal flows and scaling of luminal cross-sectional area with the volume of supplied tissue, it fails to predict fractal dimensions D and nearest-neighbor correlation coefficients r₁ of regional flows under normal physiological conditions, but gives D and r₁ values consistent with data observed under coronary stenosis and reduced coronary perfusion pressure and may thus characterize the associated disturbances of regional flows. In particular, the negative nearest-neighbor correlations predicted by model 2 may represent a vascular steal phenomenon in partially ischemic hearts, similar to coronary steal observed under vasodilator-induced conditions in myocardial ischemia (Becker, 1978).

In recent years, there has been increasing experimental evidence that the most important factor for regional flow heterogeneity in the heart is not the structure of the vascular network, but the metabolic need of the tissue itself (Van Beek, 1997; Deussen, 1998; Bassingthwaighte and Li, 1999; Balaban and Arai, 2001; Bassingthwaighte et al., 2001; Decking, 2002). In the light of these advances, the virtue of the present model lies in the approach of prescribing terminal flows as an independent boundary condition rather than obtaining these flows by means of the dispersive effect of an artificially generated vascular network.

We have studied two models for the distribution of terminal flows: Model 1 employs equal terminal flows and predicts, depending on the specific setting of the number N_term of end-points, fractal dimensions of perfusion heterogeneities in the range from 1.20 to 1.40 and positively correlated nearest-neighbor regional flows in good agreement with experimental data observed under normal physiological conditions. In model 2 we hypothesize that (a) terminal flows are proportional to the volume of supplied tissue regions, and that (b) these supply regions can be modeled as Voronoi polyhedra. Although these boundary conditions yield terminal flows well approximated by a lognormal distribution, model 2 substantially homogenizes regional flows and results in inverse correlations of nearest-neighbor regional flows, thus failing to predict D and r₁ values of the normal heart: We find D = 1.69 ± 0.02 and r₁ = −0.18 ± 0.03, independent of N_term and consistent with experimental data observed under coronary stenosis and under reduction of coronary perfusion pressure.

We conclude that flow heterogeneity can be modeled by generating terminal positions compatible with an existing tree structure and assigning terminal flows without directly resorting to the flow-dispersive effects of a specific branching tree model. As such, the current approach comes closer to the notion that metabolic tissue needs are more likely to be causative for the observed heterogeneity than the flow-dispersive effects of a vascular network. Further research is necessary to improve the model, in particular with regard to realistic settings of the total number and spatial positions of end-points, the topological connectivity of the generated trees, as well as an attempt to include a model of the main regulatory processes and signaling pathways that have been hypothesized as a possible origin of regional blood flow heterogeneity in the heart.

The authors would like to thank the anonymous reviewers of this study for their critical comments and valuable suggestions.

Olaf S. Andersen served as editor.