The 3D Spatial Autocorrelation of the Branching Fractal Vasculature

The fractal branching vasculature within soft tissues and the mathematical properties of the branching system influence a wide range of important phenomena from blood velocity to ultrasound backscatter. Among the mathematical descriptors of branching networks, the spatial autocorrelation function plays an important role in statistical measures of the tissue and of wave propagation through the tissue. However, there are open questions about analytic models of the 3D autocorrelation function for the branching vasculature and few experimental validations for soft vascularized tissue. To address this, high resolution computed tomography scans of a highly vascularized placenta perfused with radiopaque contrast through the umbilical artery were examined. The spatial autocorrelation function was found to be consistent with a power law, which then, in theory, predicts the specific power law behavior of other related functions, including the backscatter of ultrasound.


Introduction
In optics, electromagnetics, and acoustics, the behavior of waves within a weakly inhomogeneous medium can be related to properties of the spatial correlation function, a statistical measure of the spatial patterns or fluctuations within the material [1,2]. Both the forward propagating wave and the waves scattered from the inhomogeneities have distributions that have been linked to the spatial patterns of the inhomogeneities. For example, under the Born approximation and assumptions of stationarity, a key integral formula relates the ensemble average measures of backscattered intensity to the material's spatial correlation function as a 3D Fourier transform operation. A specific case emerges if we consider the parenchyma within organs such as the liver, prostate, or placenta to be the reference media, while the long, cylindrical fluid-filled, fractal branching vessels serve as the weak scattering sites.
Under that hypothesis the fractal branching structure of the fluid-filled vasculature will play a major role in determining the propagation of waves in highly vascularized tissues. However, there are few reported experimental results in the literature that can answer the key question: What is the 3D spatial autocorrelation function of a fractal branching vasculature network from normal soft tissue organs? The question has high significance because the 3D spatial autocorrelation function of these scattering sites drives the statistics of wave propagation, including the millions of ultrasound scans each year derived from backscattered echoes. We address this question directly by applying radiopaque contrast to the vasculature system of a freshly delivered, intact placenta, obtaining high resolution computed tomography (CT) scans, thresholding the 3D data set, and determining local autocorrelation functions within interior regions. In theory, these results provide the quantitative link between the properties of the fractal branching vasculature and physical measures such as scattering and fluctuations in the forward propagating waves.

Fractal Structures and Scattering
Since the popular work by Mandelbrot [3], the idea of self-similar structures occurring across a wide range of scales has been applied to natural structures. The concept of fractals in biological systems has seen numerous applications [4,5]. The nature of scattering from fractal structures has also received attention, with the concept of a fractal dimension d f driving power law relations and scattering behavior. For 3D space filling structures, the measured d f would be less than 3. Lin et al. [6] modeled fractal aggregates and essentially argued that the long wavelength (low frequency or Rayleigh scattering) limit would have the standard k 4 dependence of backscatter where the wavenumber k = ω/c, whereas the short wavelength (high frequency) limit would approach k (4−d f ) . Javanaud [7] simply invoked a power law relation for scattering using d f as the scaling power.
Shapiro [8] argued from a number of examples that a fractal structure would have a spectrum (over some range of scales) that follows a power law, and derived a scattering formula for different regimes under that assumption. One conclusion was that over some scale the backscattered intensity would be proportional to k (4+ν) , where ν is a constant related to the fractal power law behavior. Sheppard and Connolly [9] considered the optical scattering from random surfaces. Shear wave scattering and loss were attributed to fractal structures in Lambert et al. [10] and Posnansky et al. [11]. The compatibility of these works requires further study.
To summarize, the literature contains a number of theoretical functions in random media that are related to fractals and are linked to different power law relations. However, we are left with some open questions as to the particular autocorrelation function found in practice, and its implications for wave propagation. These are addressed in the following sections.

General Theory
When considering scattering from random media, it can be shown [12][13][14] that the differential scattering cross section per unit volume σ d (k) and the spatial correlation b(r) function of the inhomogeneities are related by: where k is the wavenumber,r is the vector separation between two points within the ensemble average, and A is a constant. Assuming the correlation function is isotropic and simply dependent on the separation distance r, the volume integral reduces to: similar to the integral found in Equations (4), (9) and (11) from Parker [15], covering both acoustic and electromagnetic scattering under the Born approximation. By comparing this to the form of the 3D Fourier transform in spherical coordinates with spherical symmetry 3DS {} [16,17], and shown in Appendix A, we may write: where 3DS {b(r), k} = (2π/k) ∞ 0 r · b(r) · sin(2kr)dr. Thus, the calculation of backscattering reduces to the question of the 3D Fourier transform of the isotropic and spherically symmetric ensemble-averaged correlation function of r. This spatial correlation function of cylindrical fluid-filled vessels is considered next.

The Long Cylindrical Shape
We considered fractal branching networks of the fluid-filled vasculature and other parallel fluid channels to be comprised of long cylindrical sections, branching out in successive generations from largest to finest.
Assuming an isotropic spatial and angular distribution of each generation of fractal branching structures, the goal was to average the autocorrelation function of a basic element across all angles of incidence with respect to the propagating wave and across all size scales from very small micro-channels of fluid to the largest arteries and veins that can exist within the organ. Specifically, we will examine the fluid-filled cylinder of radius a: where κ 0 is the fractional variation in density plus incompressibility, assumed to be 1 consistent with the Born formulation, F(ρ) represents the Hankel transform, which is the two-dimensional Fourier transform of a radially symmetric function, and ρ is the spatial frequency.
First, assuming the fluid-filled cylinder is long in the z-axis, then the shape is one-dimensional and its autocorrelation can be obtained from the inverse Hankel transform of the square of the shape's Hankel transform: where ρ is the spatial frequency. Second, we assumed that cylinders of this kind exist within a fractal geometry in an isotropic pattern. Therefore, copies of this were interrogated by the forward wave over all possible angles. Thus, we formed a 3D isotropic correlation function for the shape by conversion to spherical coordinates. Consider first, one infinitely long cylinder with a material property f (r)-symmetric-as shown in Figure 1a and tilted at some arbitrary angle in a spherical coordinate system. It has a 3D transform that is a thin disk (delta function k z but shown with finite thickness to make the graphic easier to draw and visualize). Its transform is symmetric and given by the Hankel transform of order zero F(ρ) = H{ f (r), ρ}. A particular radius of value of q 0 is shown for reference. Now, if we have an ensemble average of F 2 (ρ) across many cylinders at random angles (orientations) as shown in Figure 1a,b, but consider them as independent and uncorrelated (a cloud of sparsely spaced scattering cylinders), then the ensemble average 3D transform is formed from the disk function adding up over all realized angles. In the limit, any cylindrical radius q 0 in Figure 1b forms a spherical shell in the ensemble average ( Figure 1c). Now, if we have an ensemble average of ( ) 2 F ρ across many cylinders at random angles (orientations) as shown in Figure 1a,b, but consider them as independent and uncorrelated (a cloud of sparsely spaced scattering cylinders), then the ensemble average 3D transform is formed from the disk function adding up over all realized angles. In the limit, any cylindrical radius 0 q in Figure 1b is the Hankel transform of order zero [16], and ( ) δ ⋅ is the Dirac delta function. Next, the thin disk (delta function) in cylindrical coordinates represented as ( ) ( ) where the sifting property of the Dirac delta function in curvilinear coordinates is used in the z Δ and θ Δ directions, respectively. In spherical coordinates, this result is independent of θ (see Equation (17) of Baddour [17]). By setting these equal and thus independent of the coordinate system, and setting , q ρ Δ = Δ we find: Thus, the ensemble average simply takes the Hankel transform function F 2 (ρ) found in Figure 1b and populates the spherically symmetric 3D transform B s (q), shown in Figure 1c. However, there is a scaling of the cylindrical transform function over the spherical ensemble average. Specifically, the 3D spatial Fourier transform of the long cylinder oriented along the z-axis can be expressed as: where [16], and δ(·) is the Dirac delta function. Next, the thin disk (delta function) in cylindrical coordinates represented as F 2 (ρ 0 )δ(k z ) must be converted to a thin disk in spherical coordinates B s (q 0 )δ(θ − π/2) as considered in Figure 2, where the common radius ρ 0 = q 0 and ∆ρ = ∆q. Of importance is the relative scaling of the functions within the conversion to spherical coordinates.
Thus, with reference to Equation (3), we conclude that: Finally, the cross-correlation of an ensemble of these elements and all other (larger and smaller) elements within a fractal structure needs to be derived. The simplest assumption that can be made is that each scattering element has an autocorrelation function with itself that has been determined (above), and that within the ensemble average the cross terms with all other branches (larger and smaller within the fractal structure) are simply a small constant that is nearly invariant with the position and therefore can be neglected except for spatial frequencies nearing zero. Under that very simplistic assumption, the overall autocorrelation function is given by the sum (or integral in the continuous limit) of the different sizes' correlation functions over all generations of branches, weighted by their relative numbers (number density in the continuous limit). Fractal structures in 3D and volume filling in 3D may be characterized by number density functions N is a global constant, and b is the power law coefficient [18]. Within this In the limit of a small element size, the integral of each function, or Riemann middle sum, around a point q 0 on a plane orthogonal to the k z axis will be: where the sifting property of the Dirac delta function in curvilinear coordinates is used in the ∆z and ∆θ directions, respectively. In spherical coordinates, this result is independent of θ (see Equation (17) of Baddour [17]). By setting these equal and thus independent of the coordinate system, and setting ∆ρ = ∆q, we find: Thus, with reference to Equation (3), we conclude that: where Finally, the cross-correlation of an ensemble of these elements and all other (larger and smaller) elements within a fractal structure needs to be derived. The simplest assumption that can be made is that each scattering element has an autocorrelation function with itself that has been determined (above), and that within the ensemble average the cross terms with all other branches (larger and smaller within the fractal structure) are simply a small constant that is nearly invariant with the position and therefore can be neglected except for spatial frequencies nearing zero. Under that very simplistic assumption, the overall autocorrelation function is given by the sum (or integral in the continuous limit) of the different sizes' correlation functions over all generations of branches, weighted by their relative numbers (number density in the continuous limit). Fractal structures in 3D and volume filling in 3D may be characterized by number density functions N(a) that are represented by N 0 /a b , where a is the characteristic radius of the canonical element, N 0 is a global constant, and b is the power law coefficient [18]. Within this framework, the overall correlation function is: (10) where N(a) is the number density of elements as a function of size a. Alternatively, in the 3D transform domain: For the specific case of the fluid cylinder, this becomes: where f 1 (b) is a function of b, and with reference to Equation (9), we find that the predicted backscatter is: The autocorrelation function B(∆r) is also found from the inverse Fourier transform of Equation (12) to be B(∆r) = C · f (b)/r γ and where γ = b − 3 for the convergence of the inverse transform integral, 5 > b > 3, and where C is a constant. Note that the integral operators (transform, average over angles, and average over radii) are all linear so the order of these operations can be transposed for computational ease.
Our general approach is shown in Figure 3.
The autocorrelation function ( ) B r Δ is also found from the inverse Fourier transform of Equation > > , and where C is a constant. Note that the integral operators (transform, average over angles, and average over radii) are all linear so the order of these operations can be transposed for computational ease.
Our general approach is shown in Figure 3.

Methods
The fetal side placental perfusion was adapted from dual lobular placental perfusion methods previously published [19][20][21]. Briefly, human placentae from normal term deliveries were obtained within 5-10 min of delivery and examined for tears and gross lesions. Within 25 min, the umbilical artery and vein were cannulated near the insertion of the cord on the surface of the chorionic plate with five French umbilical catheters for perfusion at 10-15 mL min −1 with a finger pump. Flow rate was adjusted to maintain fetal vessel pressure at ~60 mmHg. The fetal perfusate consisted of M199 media without

Methods
The fetal side placental perfusion was adapted from dual lobular placental perfusion methods previously published [19][20][21]. Briefly, human placentae from normal term deliveries were obtained within 5-10 min of delivery and examined for tears and gross lesions. Within 25 min, the umbilical artery and vein were cannulated near the insertion of the cord on the surface of the chorionic plate with five French umbilical catheters for perfusion at 10-15 mL min −1 with a finger pump. Flow rate was adjusted to maintain fetal vessel pressure at~60 mmHg. The fetal perfusate consisted of M199 media without phenol red (Gibco) modified by the addition of dextran (35-45 kDa; 30 mg/mL fetal), D-Glucose (2 mg/mL), sulfamethoxazole (80 µg/mL), trimethoprim (16 µg/mL) gentamicin (52 µg/mL), and heparin (20 USP IU/mL). The fetal perfusate was gassed with 20% O 2 /75% N 2 /5% CO 2 in a 250 mL vessel and bubbles trapped before delivery to the placenta. The cannulated placenta was placed in a plastic bag and immersed in a 37 • C water bath for an hour followed by Doppler and ultrasound elastography experiments as described in McAleavey, et al. [22]. At the conclusion of these experiments, the placenta was perfused with a 37 • C suspension of 30% barium sulfate in 1% agarose prepared from a 60% emulsion oral contrast suspension (Barium-Liquid E-Z-Pague; Bracco) diluted with a 2% agarose in water solution with a gelling temperature of 35 • C. Perfusion continued until no further change was apparent. The placenta was then immersed in 10% neutral buffered formalin for fixation before imaging with a Philips Brilliance 64 computerized axial tomography system. The slice dimension was 768 × 768 pixels, each 0.25 × 0.25 mm; slice thickness: 0.67 mm, spacing between slices: 0.33 mm.

Results
The raw data set (shown in the maximum intensity projections in the left panel of Figure 4) was thresholded using a constant threshold limit of 180. The threshold was set at the minimum level necessary to zero out the poorly vascularized edges of the placenta, while maintaining well connected branches. Projections of the thresholded binary data can be viewed in the right panel of Figure 4. Also shown in the right panel of Figure 4 is the convex hull of the placenta (solid black line) as well as the boundary of the sampling region (dashed black line). The sampling region contains locations whereby small autocorrelation sample cubes (6 mm on a side, or 27 × 27 × 21 voxels) can be translated along any axis or diagonal by 6 mm without being outside the whole placenta's boundary. For reference, the size of the autocorrelation cubes is shown in the bottom right corner of the right panel of Figure 4, and the surrounding region that contributes to the autocorrelation spectrum is shown by the red circle of radius 9 mm. thresholded using a constant threshold limit of 180. The threshold was set at the minimum level necessary to zero out the poorly vascularized edges of the placenta, while maintaining well connected branches. Projections of the thresholded binary data can be viewed in the right panel of Figure 4. Also shown in the right panel of Figure 4 is the convex hull of the placenta (solid black line) as well as the boundary of the sampling region (dashed black line). The sampling region contains locations whereby small autocorrelation sample cubes (6 mm on a side, or 27 × 27 × 21 voxels) can be translated along any axis or diagonal by 6 mm without being outside the whole placenta's boundary. For reference, the size of the autocorrelation cubes is shown in the bottom right corner of the right panel of Figure 4, and the surrounding region that contributes to the autocorrelation spectrum is shown by the red circle of radius 9 mm.   From within the sampling region, 1000 locations were chosen randomly as the center coordinate of an autocorrelation window. Shifts were applied over all possible translations up to 6 mm, and the autocorrelation function B(∆r) = I(x, y, z) * I(x + ∆x, y + ∆y, z + ∆z) was calculated, where I(x, y, z) is the thresholded computed tomography (CT) intensity, as a binary function {1, 0}, and ∆r = ∆x 2 + ∆y 2 + ∆z 2 . The normalized autocorrelation B(∆r)/B(0) was then calculated and used in subsequent analyses.
Individual locations' autocorrelation functions were found to be dependent on the number of samples included within the 6 mm local cube, and this is related to the structure(s) included. For example, a box centered on a large arterial segment has many pixels above threshold, and has a rather slow decorrelation since the object I(x) is large. At the other extreme, if a box is centered on a region with no obvious vascular architecture and only a few small isolated points above threshold, the decorrelation is steep as a function of ∆r. These extremes point to the inherent difficulties in characterizing a multi-scale process with limited sample volumes and limited resolution. For many sample locations, however, a branching structure was included, and this corresponded to numbers of voxels N above threshold between approximately N = 100 and 1000 (out of 15,309 voxels per cube). These characteristic types are shown in Figure 5.
There is also a theoretical question relating to the nature of the autocorrelation function: Can it reasonably be described as a power law function? To address this, power law functions were fit for autocorrelation lags ∆r from 0 to 6 mm. The results of the power law coefficient for B(∆r) = B 0 · (∆r) γ are plotted in Figure 6. The minimum mean squared error fit to a power law provides an exponent of approximately −1.3. Some indication of a flattening at larger separation distances, greater than 4 mm, was seen however larger multiscale analyses would be required to determine if that was a general result indicating a tendency of the branching vasculature correlation function at longer ranges. decorrelation since the object ( ) I x is large. At the other extreme, if a box is centered on a region with no obvious vascular architecture and only a few small isolated points above threshold, the decorrelation is steep as a function of r Δ . These extremes point to the inherent difficulties in characterizing a multiscale process with limited sample volumes and limited resolution. For many sample locations, however, a branching structure was included, and this corresponded to numbers of voxels N above threshold between approximately N = 100 and 1000 (out of 15,309 voxels per cube). These characteristic types are shown in Figure 5. There is also a theoretical question relating to the nature of the autocorrelation function: Can it reasonably be described as a power law function? To address this, power law functions were fit for autocorrelation lags r Δ from 0 to 6 mm. The results of the power law coefficient for ( ) are plotted in Figure 6. The minimum mean squared error fit to a power law provides an exponent of approximately −1.3. Some indication of a flattening at larger separation distances, greater than 4 mm, was seen however larger multiscale analyses would be required to determine if that was a general result indicating a tendency of the branching vasculature correlation function at longer ranges. We find some variation in the power law fits for the complete range of lags that was dependent on the lower and upper bounds on voxel counts chosen, as shown in Table 1. The bounds were varied to exclude cases that were extreme in either containing too few voxels above thresholds or too large of a structure to be meaningful, however the final average did depend on the particular values chosen, but clustering around for a reasonable range of populated autocorrelation boxes. We find some variation in the power law fits for the complete range of lags that was dependent on the lower and upper bounds on voxel counts chosen, as shown in Table 1. The bounds were varied to exclude cases that were extreme in either containing too few voxels above thresholds or too large of a structure to be meaningful, however the final average did depend on the particular values chosen, but clustering around γ = −1.3 for a reasonable range of populated autocorrelation boxes.

Consequence of the Model and Power Law
For soft tissue scattering, a number of studies have characterized the optical and acoustical measurements of tissue specimens. Specifically, most carefully calibrated ultrasound backscatter results from soft tissues report increasing backscatter with frequency. Generally this can be fit to a power law with an exponent of 1 < γ < 2. For example, Campbell and Waag found liver backscatter to increase as f 1.4 over the medical imaging band of 3-7 MHz [23] consistent with other pioneering reports from the early 1980s [24][25][26][27][28]. Other examples of increasing backscatter vs. frequency in livers are demonstrated in [29] extending up to 25 MHz, and in Lu et al. [30], also in livers. Similarly, for other tissues, including sarcoma and carcinoma models, the increase in backscatter vs. frequency was reported in Oelze and Zachary [31] and for thyroids in Rouyer et al. [32]. Thus, in assessing our theoretical model, Equation (3), we are seeking a prediction of a backscatter power law fit of an exponent slightly larger than one, given a relatively isotropic fractal fluid network scattering power law of b − 2. In that case, it is clear that the power law b governing the number density of vessels must be greater than 3. This is plausible, although reported data has to be considered in light of the methodologies used; this was considered in Section 5.2. We also note for completeness that for F[q] 2 /q in spherical radial frequency being proportional to q (b−6) , also corresponds in the spatial domain to an autocorrelation function proportional to 1/r b−3 s , where 3 < b < 5. For our placenta measurements, γ = −1.3 = 3 − b implies b = 4.3, also implying backscatter would increase as k 2.3 from Equation (13), over some frequency range.
Although our resolution was limited by the CT system and contrast enhancement, optical studies can obtain finer resolution. Schmitt and Kumar [33] studied phase contrast images of mouse liver histology slides and found a power law power spectrum behavior for spatial frequencies down to at least 1/10 microns.
Similarly, Rouyer et al. [32] studied thyroid tissue using broadband ultrasound systems and by replotting their data we can estimate a power law behavior between 6 MHz and 15 MHz, corresponding to wavelengths down to 100 microns.

Refinement of Theory vs. Experiments
A more exact comparison of theory to experiment will require greater precision on two fronts, the backscatter vs. frequency power law, and the fractal branching power law of the fluid channels, for specific organs. For example, with respect to liver, Campbell and Waag [23] reported a backscatter proportional to the frequency to the 1.4 power over a range around 5 MHz for beef liver samples. These measurements and the analytical treatment of all system effects were highly rigorous, yet confirmation of these results over a larger bandwidth and in conjunction with imaging of the samples (ensuring avoidance of any unwanted ligaments or trapped bubbles, for example) would help to refine this estimate.
The fractal branching behavior assumed in Equation (10) is a continuous number density of cylinders represented as proportional to N 0 /a b ; this forms a major assumption of this framework, yet is not known precisely. In reviewing the literature results, it must be kept in mind that there are major distinctions running through different analyses, principally 2D vs. 3D, and also the type of measurement utilized. Fractal dimension is limited by topographical dimension, so 2D analysis of slices or projections will have a lower dimension than corresponding 3D measurements. Full 3D measurements at high resolution are rare and difficult due to the demands on resolution and sheer size of the imaging data. Some high quality and high resolution 3D assessments of the branching vasculature have been rigorously acquired for the brain [34] and the placenta [35], however to our knowledge there have not been any published results from 3D fractal analyses of the liver vasculature at high resolution.
Adding to the uncertainty about the important power law parameter b, which defines the number density of cylindrical vessels as a function of their radius, over the ensemble and in linear dimensions, is the wide variety of measurement styles that are typically reported. This is a major source of possible misinterpretation. Our power law parameter b, which governs the number density of vessels as a function of diameter, is not the same as the fractal dimension measured by sandbox or related box counting techniques. Even in the cases of older literature reports where the number of vessels were painstakingly counted (see Figure 1 of [18] indicating a power law of approximately −2.7 over some older studies, Table 1 of [36] indicating a power law of approximately −1.3 in pulmonary veins, or Table  4 of [37] suggesting a power law of −2.5 in pulmonary arteries), there are important methodological details that can swing the assessment. For example, many counting schemes derive the number of vessel vs. branch generation number 1, 2, 3 . . . as ordinal numbers. The slope of this curve is not the same as our number density vs. radius curve. Other schemes count "greater than or equal to" a varying radius, which corresponds to an integration. Since the integration of 1/a b is proportional to 1/a (b−1) , this scheme inherently shifts and reduces the power law parameter by one. Similarly, a scheme of "binning" together vessels by proportional limits over a log scale will shift the power law. For example, counting all the vessels within ±15% of 0.1 mm, then all those within 15% of 1 mm, then all those within 15% of 10 mm, will effectively integrate the distribution within these bins, again converting the continuous distribution to a discrete set with 1/a (b−1) relative distribution. For all these reasons, it is plausible that some literature assessments of branching vascular parameters, using different methodologies, may not be the same as number densities defined in our manner by the parameter b. Clearly, high resolution data sets in 3D are required to better define this important factor in soft tissues such as liver, thyroid, prostate, and others. The particular 3D shapes and autocorrelation functions from these different organs may differ from the results found in the placenta, since the 3D placenta has a unique relationship with the maternal uterine blood flow, so organ-specific parameters are needed to form baseline values.

Conclusions
High resolution contrast CT studies of a normal placenta's branching vasculature in 3D demonstrated a power law behavior for the 3D spatial autocorrelation function with an exponent close to −1.3. Correspondingly, a theoretical hypothesis about scattering from tissue was generated from a primitive cylindrical shape representing the plausible distribution of extracellular fluids and blood in long fluid-filled channels throughout normal soft tissue. Assuming a wide range of diameters of the cylindrical fluid spaces and a macroscopically isotropic distribution over some region of interest within the organ, the predicted backscatter is of the form of a power law with dimension greater than one. This matches some observations about the nature of vessels and measurements of backscatter from soft tissues and is consistent with the power law autocorrelation function measured experimentally from the placenta.
Pathologies that affect tissue morphology would be expected to modify this model in significant ways, but it is prudent to first establish an accurate working model of normal soft tissue scattering. More precise measurements of the key parameters will be required to fully test this hypothesis and prior models based on spheres.