Fractal ¹⁵O-Labeled Water Washout From the Heart

Analytical Methods and Indicator Dilution Theory

Basic conservation theory, as summarized by Zierler,¹⁵ maintains that after an impulse injection at t=0 into the entrance to an organ, the fractional residual content R(t) at the next instant is 100% of the injectate. R(t) remains at 1.0 until there is an appearance of tracer in the outflow. The fraction of the injectate appearing in the outflow per unit time is h(t) (fraction per second). The accumulated outflow is the integral of the outflow response and is the total dose, unity, minus the fraction of dose still retained within the organ, R(t): where λ is a dummy variable used in the integration. Further, in the ideal situation, in the absence of error in the data acquisition, the outflow fraction of the dose appearing per second is the derivative of the organ content: In practice, one records neither R(t) nor h(t) because injecting an impulse input is impossible. Instead, the signal recorded from the heart, Q(t), is a convolution of R(t) with the input: where F̅ is the mean flow, and the asterisk denotes the convolution integration, more formally written as Q(t)=F̅∫_0^tC_in(λ) · R(t−λ)dλ. Likewise, the outflow concentration-time curve, C(t), is as follows: When the input curve, C_in(t), is completed within a period that is short compared with the mean retention time within the organ, then the tail portion of h(t) has a form influenced almost solely by the network characteristics and scarcely at all by the form of C_in. (For these experiments, the input was 99% complete in 18 to 22 seconds, which was short compared with washout, which took 3 to 6 minutes.) Given these conditions, a good approximation for the residue function is as follows: where Q_max is the maximum value of the curve recorded via the coincidence counters. Although the approximation is poor in the initial seconds, it is good thereafter. Likewise, the approximation for the tail of h(t) is the normalized outflow concentration-time curve: where the denominator is the area under the curve C(t). Alternatively, h(t) can be calculated by using the injected dose, q_o, and the measured flow, F̅, so that h(t)=F̅C(t)/q_o, illustrating that h(t) is the fraction of the dose reaching the outflow per unit time.

The fractional escape rate, η(t), is defined as the fraction of residual tracer escaping per unit time, and is therefore given directly by the amount of tracer appearing in the outflow per unit time, F̅C(t), where F̅ is the total flow, divided by the amount of tracer retained within the organ at each time t: This review of standard indicator dilution theory and of the secondary approximations sets the stage for the special approach to the analysis using fractals. Note that there are two assumptions in the normalization of Equation 9 in going from the middle to the right expression: (1) Q_max in Equation 7 is obtained at a time when all of the tracer is in view of the coincidence detectors, with all of the tracer having entered and none having left. (2) The area of the outflow washout is essentially complete. Assumption 1 seems reasonable in view of the fact that the input curves, C_in, are 99% complete by 18 to 22 seconds (not illustrated), although R(t) does have a measurable plateau since the washout begins at about this time. Assumption 2 is valid within <2%, as can be seen from the final values of R(t) at t=300 seconds in Fig 1⇓.

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

Residue curves, R(t), and normalized outflow concentration-time curves, h(t), from a rabbit heart after injection of ¹⁵O-labeled water into the coronary inflow. Both R(t) and h(t) have tails that are power law functions of time; ie, they are straight on log-log plots (right) and clearly not monoexponential, being concave upward on semilog plots (left).

A Fractal Theory To Be Tested

The fundamental feature of a fractal is self-similarity or self-affinity. A fractal washout process is therefore one for which the rate of washout decreases by some exact proportion for some chosen proportional increase in time; the self-affinity requirement is fulfilled whenever the “exact proportion” remains unchanged, independent of the moment or the segment of the data set selected to measure the proportionality constant. The length of time or the portion of the washout curve for which this relation holds is known as the “scaling region” and is necessarily finite, as is the case for all natural, nonmathematical, fractals. Fractal washout behavior can begin after the input to the organ is complete. If complete, then one would expect the late portion of R(t) to have the shape of a power law relation defining the proportionality relation: where α is the power law exponent. (The observation of a power law relation is not an explanation of the phenomenon but provides an incentive for a search for an explanation.) This power law relation defines the log-log slope dlogR/dlogt=−α. The corollary of this is that the outflow tracer concentration-time curve should also be fractal and since h(t)=−dR/dt, then taking the derivative of the right side of Equation 10 yields and the slope on the log-log plot is dlogh(t)/dlogt=−α−1.

According to the hypothesis, expressed in Equations 10 and 11, η(t) should be a power law function with an exponent equal to −1: Three equations pose the tests of the fractal washout hypothesis: Over some substantial period of time late in the washout, (1) does Equation 10 provide a good fit to R(t), (2) does Equation 11 provide a good fit to h(t), and (3) can one value for α be used to fit Equations 10 and 11 to the data simultaneously obtained for R(t) and h(t)? If so, then does Equation 12 provide a good fit to the data for the fractional escape rate η(t), with the specific power law slope of −1; ie, does η(t) decay at a rate proportional to 1/time?

Unconstrained Fitting of R, h, and η

The tails of the residue R(t), outflow h(t), and the escape rate η(t) for two experiments are plotted on log-log scales in Fig 2⇓. Each curve was fit with the log-log regression to give best estimates of the power law exponent. The values from 18 experiments for α from the residue functions, fitting Equation 10, were 2.12±0.33 (mean±1 SD); the values for α from the outflow curves, fitting Equation 11, were 2.13±0.34. These values were obtained independently from residue and outflow data. The coefficients of variation for the fits to the R(t)s averaged 0.073±0.042 (N=18); the coefficients of variation for the fits to the outflow curves averaged 0.074±0.070 (N=18). Thus, the independently estimated values of α from residue and outflow were close to each other (2.12 and 2.13). The paired differences averaged 0.005, with a standard deviation of 0.43 for the 18 pairs. Thus, the first two tests listed below Equation 12 are satisfied and do not cause us to reject the fractal relation, since the two independent estimates of α are not statistically different, as tested by either unpaired or paired differences by Student’s t test.

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Fitting power law curves of general form kt^−β to residue and outflow curves, R(t) and h(t), respectively, and to the fractional escape rate, η(t), for two separate experiments. These are unconstrained best power law regressions fitted to the individual curves for R(t), h(t), and η(t).

The curves for the fractional escape rate η(t) are necessarily noisier than those for either h(t) or R(t), since η(t) is calculated by dividing one curve by the other (Equation 9). Because of the noise, the coefficients of variation for the power law regression fitting are larger (0.21±0.16, N=18). The power law exponents obtained for unconstrained fitting of the equation η(t)=kt^−β gave estimates of β=1.12±0.47 (N=18), so that the average of the unconstrained estimates of β did not differ statistically from 1.0, the value theoretically anticipated in Equation 12. Again, the fractal hypothesis cannot be rejected.

Comparing Power Law With Monoexponential Fits

Fig 3⇓ shows optimized best fits using both the power law expressions, R̂(t)=a₁t^−α₁ and hˆ(t)=a₂α₂t^{−α₂₋₁}, and the exponential expressions, R̂(t)=A₁e^−k_1t and hˆ(t)=A₂e^−k_2t. The test was designed as a test of appropriateness of extrapolation to predict the shape of the last parts of the tails of R(t) and h(t); to do this, the model functions were fit to only a limited segment of the data: for R(t), from 0.2 to 0.1 only; for h(t), over exactly the same time period as for R(t). For this test, the power law fits were not constrained to use the same α for R(t) and h(t) but were best-fitting regressions, log-log for the power law and log-linear for the exponential. The rate constants for the exponential fits represent best fits over the same segments of R(t) and h(t) as were used for the power law fitting.

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Extrapolation of residue and outflow curves, R(t) and h(t), respectively, with power law versus exponential fitting. Fitting was done over the time period where 0.2≤R(t)≤0.1 for both R(t) and h(t). The power law and exponential curves beyond R(t)=0.1 represent the extrapolated predictions beyond 60 seconds (top) and 50 seconds (bottom).

The two sets of data shown in Fig 3⇑ on semilogarithmic, not log-log, plots are representative of the predictive capacity of the two extrapolation methods. The power law extrapolation beyond the time where R(t)=0.1 always lies above the exponential best fit for both R(t) and h(t). For 15 of 18 R(t)s, the power law extrapolation predicted and fit the data for R(t)≤0.1 better by far than did the exponential extrapolation; in two sets, the data sets were not long enough beyond R(t)=0.1 to make the distinction, and in one set, R(t) was better fit by the exponential, but h(t) for the same set was much better fitted by the power law. For h(t) the results are less secure, because the choice of the region fitted, being defined by R(t) and not by h(t), means that the curves are relatively noisy. For 12 curves, the power law fit was better by far; for 4 curves, the exponential fit was as good; and for 2 curves, the exponential fit was distinctly better. Of the 4 intermediate results, the data in 2 cases lay clearly in between the 2 extrapolations (power law and exponential) and did not fit either, and in 2 cases, the data points were too scattered to make the distinction.

The areas under the tail of R(t) were underestimated by the monoexponential extrapolation, with the ratio to that estimated by the power law being 0.64±0.12 (N=18). For h(t), the ratio of areas by exponential fit over that for power law was 0.63±0.08 (N=18).

From the point of view of estimation of total areas of curves, this difference is not great, a 3% to 4% underestimation for R(t) and less than that for h(t). The errors are of course greater for mean transit times, where the tails have more influence on the estimate.

Using One α for Fitting R, h, and η

The most stringent test of the hypothesis is test 3. In this approach to obtaining the best estimate of α for each experiment, the data for R(t) and h(t) were fit-ted simultaneously by using R̂(t)=a₁t^−α and ĥ(t)=a₂αt^−α−1. By this approach, the hypothetical fractional escape rate must have the form η(t)=a₃(a₂/a₁)αt⁻¹; the exponent is forced to be −1. The deviations of a₃ from 1.0 represent the difference in the sensitivities of the outflow detector (for positron emissions) and the residue detector (for 511-keV coincident gamma emissions). The best estimates of the a_i values and α were found by nonlinear optimization using sensop, a sensitivity function–based modified Levenberg-Marquardt type routine (Chan et al¹⁶ ). In this approach, the fractal hypothesis of a difference of 1.0 between the power law exponents for R(t) and h(t) is assumed, and with this assumption having been made, the test lies in the goodness of fit to the data for the set of three curves. Two examples are shown in Fig 4⇓. The values of α so estimated are reported in the Table⇓. The values of α range widely around 2.0, averaging 2.2±0.27 (mean±SD, N=18). The constrained theoretical lines appear to be close enough to the three data sets so that the hypothesis cannot be rejected. The inference from both the slopes of R(t) and of h(t) and of the relation between their slopes is that the tails of the transit time distributions are approximately fractal in these isolated perfused rabbit hearts.

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Optimized best fit of combined fractal power law relation for residue R(t)=a₁t^−α and outflow h(t)=a₂αt^−α−1 for the same value of a, forcing the fractional escape rate η(t)=a₃(a₂/a₁)αt⁻¹. Top, The exponent α was 2.22, and the parameters a₁, a₂, and a₃ were 13485, 9396, and 1.01, respectively. The 95% confidence range on α was 2.22±0.01. Bottom, α was 2.41±0.03 (95% confidence limits), and a₁, a₂, and a₃ were 1342, 1178, and 1.15, respectively.

Generally, myocardial residue curves exhibit similarity upon scaling time by dividing by mean transit time.² This would be the same as scaling by a flow-to-volume ratio or simply by multiplying time by F when V is constant. One might intuitively expect a higher value of α at higher flows, and this is suggested by the positive slope of α versus F in Fig 5⇓, but it is not definitive.

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Fractal exponent (α) versus myocardial blood flow (F).

In Fig 6⇓, a more specific test is performed with respect to the question of whether α is related to flow. Three sets of curves were obtained from one heart at flows of 0.76, 1.34, and 3.44 mL · g⁻¹ · min⁻¹. The three residue curves have been superimposed on each other by scaling time with respect to mean transit time (Fig 6⇓, upper curves). The three outflow curves were treated likewise (Fig 6⇓, lower curves). The close juxtaposition of the curves upon each other shows that their shapes are similar, meaning that their higher moments (variance and the skewness and kurtosis obtained from the third and fourth moments) are close to being the same. In Appendix A, this time scaling is applied to the function h(t)=Aαt^−α−1, and it is shown that the exponent α is the same on the time-scaled h(t/t) as on the original h(t). Since scaling by flow should be equivalent to scaling by time, the theory would predict that in one heart at three different flows, the residue and outflow curves should exhibit the same power law exponent α. In this experiment (the one with the widest range of flows in the present study), the α values were 1.82, 1.91, and 2.03; ie, there is a 10% difference between the lowest and the highest. The trend for each individual heart in which there is a range of flows is like that for the whole data set in Fig 5⇑. Thus, although the idea that α is related to flow is neither refuted nor supported by such data, there is room for suspicion despite the theory in Appendix A. Both Fig 5⇑ and the analysis of the similarity-scaled data sets of Fig 6⇓ leave open to further study the question of whether or not α increases with flow.

• Download figure
• Open in new tab
• Download powerpoint

Figure 6.

Residue (upper three) and outflow (lower three) curves, R(t) and h(t), respectively, for ¹⁵O-labeled water from one heart superimposed by scaling the time axis by individual mean transit time, t̅. There is close similarity in the shapes of the curves obtained at 0.76, 1.34, and 3.44 mL · g⁻¹ · min⁻¹.

The absolute levels of myocardial blood flow influence the washout curve, h(t), with higher flows giving an earlier peak and a shorter mean transit time. The fractional escape rate of tracer from the organ is higher at high flows. What is interesting about the fractional escape rate is the exponent: the power law slope of −1 is independent of the flow and of the shapes of R(t) and h(t). This universality is striking: whenever the system is fractal, the relation between R(t) and h(t) gives rise to η(t)=at⁻¹, a power law fractal such that in all situations the escape rate diminishes as t⁻¹. This is a strong statement because η(t) is equal to the derivative of R(t) and also to h(t) divided by 1.0 minus its integral. What this means is that when η(t) has an exponent of −1, both R(t) and h(t) must be power law functions if either one is. The provocation provided by this observation is to find the physical explanation, as discussed below, and to ask whether it is necessarily fractal or allows some other descriptor.

Fractal Extrapolation

For the practical purpose of measuring cardiac output, Hamilton et al¹³ proposed using monoexponential extrapolation of the downslope of dilution curves to obtain an estimate of the area under the primary first-pass indicator uncontaminated with recirculated indicator. This was an excellent and powerful suggestion and has been the method used ever since, even by those who recognized that it resulted in a minor underestimation of the area under the curve when the tails deviated (always upward) from monoexponential. Their technique is simple and allows an analytical calculation of the area under the extrapolated tail beyond the section of the downslope used to estimate the exponential rate constant. Now we propose an equally simple method of extrapolation, but we use the power law expression to exclude the recirculation while accounting for all the first-pass indicator and avoiding the systematic underestimation of the exponential extrapolation.

The method is applied to either outflow curves, h(t), or to residue curves, R(t). The power law exponent is the negative of the slope of the regression of log h(t) versus log t for the tail region where the relation is a straight line: log h(t)=log A−β log t. (For the regression between times t₁ and t₂, an estimate of β is log [h(t₁)/h(t₂)]/log [t₁/t₂].) The area, area 1, up to the time t₂ at which the extrapolation begins, can be obtained directly from the data. Area 2, under the extrapolated continuation of the observed data beyond t₂, is calculated analytically: [The condition, β>1, obviously holds for both R(t) and h(t) in the data of these experiments. The integral does not converge when β≤1, which means the area extrapolated would not be finite.] The estimate of cardiac output or flow, F, for a single path system into which a bolus of indicator of amount q_o is injected is F=q_o/(area 1+area 2).

For the residue and outflow curves in the present study, one can calculate the degree of underestimation of the areas that would result from using a monoexponential extrapolation. The calculation is based on using t₂ as the time where R(t₂)/R_peak=0.1, and for the purpose of estimating β, the beginning of the region fitted was t₁, the time at which R(t₁)/R_peak=0.2. The outflow h(t) was fitted over the identical time period. For the 18 data sets, area 2 was calculated two ways: first as Equation 13, the fractal extrapolation, and second by monoexponential extrapolation. The monoexponential rate constant was determined from the best fitting linear regression of log h(t) or log R(t) versus t (linear). The results of the two extrapolation techniques give necessarily smaller values for the monoexponential approach: For h(t), area 2 (monoexponential)/area 2 (fractal)=0.637±0.12 (N=18), and for R(t), area 2 (monoexponential)/area 2 (fractal)=0.630±0.08 (N=18), where ±1 SD is given. Naturally these errors in area due to the use of monoexponential extrapolation are much smaller than are errors in estimated mean transit times.

Why Is Washout Fractal?

The observed fractal washout may be explicable on the basis that regional flows per gram of tissue have fractal spatial distributions, as is well documented.^{4 5} The link between the spatial and temporal events is not yet defined by any theory, but it makes sense that transit times through a network with fractal flow distributions should be fractal, since regional transit times are the local volumes of distribution divided by the local flow.

It is interesting that the observed fractal exponents are possibly dependent on flow. A theory for thinking that α and F should be independent when washout curves can be superimposed on each other by proportional time scaling^{2 20} is given in Appendix A. The upward trend of α versus F in Fig 5⇑ may be a hint that the fractal model is imperfect, which is no doubt the case, and some deviation toward multiexponential form is causing some degree of flow dependency.

The observation that washout is fractal does fit with the fractal paradigm newly recognized to apply to many aspects of biology. Self-similarity over a wide range of scales is found in time-dependent functions of many sorts. As in all physical systems, the summarizing descriptor, “fractal,” applies over only a finite range. Fractals are not forever, except in the mathematics of the ideal, for every real fractal relation fails at both the large and small ends of the scale. The text of Bassingthwaighte et al²¹ presents many examples of these self-similar scaling relations, all of them limited to finite ranges. Fractal power law scaling can emerge from any set of processes repeated multiple times with appropriate scaling between members of the set. See Appendix B for the general theory and an example using multiple exponentials.

It is intriguing that the coronary vascular network has so many fractal characteristics. The anatomy itself is fractal in terms of segment diameters and lengths,²² the number of branches is a fractal function of the diameters of the vessel from which the branches derive,²³ the ratios of diameters of parent to daughter branches follow log-log scaling, and the flow distributions are spatial fractals.^{4 24} The degree of heterogeneity of regional flow distributions has been approximated by those calculated from artificial fractal branching patterns by using an overly simple dichotomous branching.²⁵ Both Gan et al²⁶ and VanBavel and Spaan⁸ estimate spatial distributions of flows from anatomic information; such calculations performed in our laboratory on artificial coronary systems give similar results and exhibit the same fractal correlation structure in regional flows in the heart as was observed in nature.²⁷ The present study adds to the story but uses the completely independent evidence on washout kinetics in vivo, as opposed to structural considerations or spatial distributions of flow at a given moment of observation. Washout kinetics, quite independently of spatial patterns or branching network structures, drive us to the same conclusion that the system is fractal.

A: Similarity of Time Scaling Means That the Fractal Exponent Should Not Be Influenced by Flow

The question of why the power law exponent α of the transport functions obtained at different flows is independent of flow can be explained mathematically on the basis of previous observations. Similarity on scaling of the time axis by the mean transit time is usually observed for the coronary system.^{2 20 28} The test of “similarity” is a statistical one; namely, can the shape of an impulse response h₁(t) obtained at flow F₁ be considered similar to the response h₂(t) obtained at flow F₂? Because there may be volume changes between two different physiological states, rather than scaling by F₁/F₂, one uses the more general scaling factor, the mean transit time or the ratio of mean transit times, t̅₁/t̅₂, where t̅_i=V_i/F_i and where the i indicates a particular state or condition or time of day when the observation of h_i(t) was made. When similarity holds, then all time-normalized impulse responses have the same shape and are superimposed on each other on a plot of t̅_ih(t/t̅_i) versus t/t̅_i.

For the special situation where h(t) is a power law function, as in Equation 11 in the text, the following equation applies: Then we can use t̅_i=V_i/F_i for clarity: Thus, when similarity scaling holds and h(t) is a fractal power law function, then the similarity scaling by mean transit time results merely in a scaling of the observed h(t) by (t̅)^−α, and the power law exponent α is unaffected by transformation. Since the scaling transformation can be performed either by using t̅_i generally or by using F_i when V_i is constant, then it follows that α should be unaffected by flow, as observed in Fig 5⇑.

The same logic holds for R(t) and η(t), which have absolute values influenced by flow but have power law exponents that are not influenced by flow.

B: Sums of Scaled Functions Can Give Power Law Behavior

A power law function can be represented as the sum of a finite number of fractal-scaled basis functions. Consider approximating of the power law function of Equation 16 with the weighted sum of basis functions f(t) in Equation 17: where a_i is the amplitude scalar and k_i is the time scalar for the ith member. Since the basis functions are not necessarily orthogonal, a finite sum of N scaled basis function is considered.

The minimum mean-squared error between F(t) and a particular f(k_it) over the interval from t=0 to t=∞ is found by calculating a_i: From this, one can solve the relation between a_i and k_i by using a dummy variable, τ=k_it, substituted into Equation 18: or where C is a constant that does not depend on k_i.

A power law function can therefore be represented by a finite sum of the scaled basis functions, where the weight of each basis function is determined by the scale factor raised to the power law exponent:

In general, the k_i can be chosen on the basis of the interval over which the power law slope is fit. If the interval is defined by t=t_a to t=t_b, then k₁ can be chosen by k₁=1/t_a or a conveniently chosen value. In order to evenly distribute all of the k_i in the log-time domain, the rest of the k_i can be calculated over the range chosen:

An example using exponentials as the basis function is demonstrated in Fig 7⇓. F and f are given as follows:

• Download figure
• Open in new tab
• Download powerpoint

Figure 7.

Log-log plot showing multiexponential fits to F=t⁻², where F is flow and t is transit time, using two, three, and four exponentials with t_a=1, t_b=100, k₂/k₁=100 for N=2, k₂/k₁=10 for N=3, and k₂/k₁=4.642 for N=4, where k is the time scalar, and N is the number of exponentials.

The finite-sum approximation is shown for N=2, 3, and 4 exponentials. An approximate fit is achieved using only four exponentials over the interval of t_a=1 to t_b=100. Making t_a and t_b outside of the desired region to be fitted and increasing N allows one to approach exact power law behavior arbitrarily closely.