|Home | About | Journals | Submit | Contact Us | Français|
The purpose of the present study was to determine the accuracy and the sources of error in estimating regional myocardial blood flow and vascular volume from experimental residue functions obtained by external imaging of an intravascular indicator. For the analysis, a spatially distributed mathematical model was used that describes transport through a multiple-pathway vascular system. Reliability of the parameter estimates was tested by using sensitivity function analysis and by analyzing “pseudodata”: realistic model solutions to which random noise was added. Increased uncertainty in the estimates of flow in the pseudodata was observed when flow was near maximal physiological values, when dispersion of the vascular input was more than twice the dispersion of the microvascular system for an impulse input, and when the sampling frequency was <2 samples/s. Estimates of regional blood volume were more reliable than estimates of flow. Failure to account for normal flow heterogeneity caused systematic underestimates of flow. To illustrate the method used for estimating regional flow, magnetic resonance imaging was used to obtain myocardial residue functions after left atrial injections of polylysine-Gd-diethylenetriaminepentaacetic acid, an intravascular contrast agent, in anesthetized chronically instrumented dogs. To test the increase in dispersion of the vascular input after central venous injections, magnetic resonance imaging data obtained in human subjects were compared with left ventricular blood pool curves obtained in dogs. It is concluded that if coronary flow is in the normal range, when the vascular input is a short bolus, and the heart is imaged at least once per cardiac cycle, then regional myocardial blood flow and vascular volume may be reliably estimated by analyzing residue functions of an intravascular indicator, providing a noninvasive approach with potential clinical application.
Noninvasive assessment of regional myocardial blood flow would be of great use for planning the therapy and monitoring the response in the treatment of coronary heart disease. A general approach involves imaging regional time course data on blood-borne indicators as they pass through myocardial tissue after intravascular injection. Imaging modalities with adequate spatial and temporal resolution for the analysis of tissue residue function kinetics include single-photon emission computed tomography (22), positron emission tomography (24, 27, 31), X-ray computed tomography (CT) (10), and magnetic resonance (MR) imaging (MRI) (23, 34). To obtain quantitative information on myocardial blood flow, tissue indicator curves must be analyzed using an appropriate mathematical model.
The central volume principle states that the mean transit time of an indicator in a vascular system is equal to the volume of distribution of the indicator divided by the flow. Estimating regional blood flow (F), independent of vascular volume, by use of the central volume principle to analyze externally detected residue function data is difficult, because it is necessary to obtain measures of regional mean transit time () and regional blood volume (V) to calculate F = V/. With use of the area under the tissue curve divided by the height of the curve to calculate mean transit time (35), dispersion of the vascular input and indicator recirculation cause error in estimates of regional flow (6). Because of normal variations in regional tissue volumes (15), it is necessary to estimate volume in every tissue region in which flow is estimated. Regional transit times cannot be directly calculated from imaging data, because the first moment of a tissue time-activity curve is not equal to the mean transit time through the organ (32). As a result of these complications, estimates obtained using the central volume principle, which might be considered model independent, may show systematic error (6, 13).
Highly permeable indicators offer the highest sensitivity for assessing regional blood flow, because their exchange is flow limited in a large volume of distribution, and the positron emitting agent, [15O]water, has been exploited for this reason with use of the positron emission tomography technique (9). However, highly permeable indicators are not available for CT and MRI, and, instead, low-molecular-weight extracellular indicators have been used for CT (11) and MRI (34). Analysis of extracellular indicators is complicated, because exchange between vascular and interstitial regions must be accounted for to obtain quantitative flow estimates. Therefore, intravascular indicators may prove advantageous, because their restriction to the vascular region simplifies the analysis, even though it also results in reduced imaging sensitivity.
The purpose of the present study was to evaluate the possibility for quantitative estimation of regional myocardial blood flow and vascular volume with use of noninvasive imaging of an intravascular indicator by determination of fundamental sources of error inherent in the approach independent of technical limitations associated with a specific imaging modality. Residue function kinetics of an intravascular indicator were analyzed using a spatially distributed mathematical model, that accounts for flow heterogeneity and vascular architecture. To test the accuracy of the analysis in the presence of signal noise, flow heterogeneity, increased dispersion of the vascular input, and decreased sampling frequency, the model was used to analyze “pseudodata”: realistic solutions generated by the model to which random noise was added. The model was also used to estimate regional myocardial blood flow by analyzing residue functions obtained using MRI during bolus injections of the intravascular contrast agent polylysine-Gd-diethylenetriaminepentaacetic acid (DTPA) in anesthetized dogs (33).
The model of the vascular system describes dispersion and delay of the externally detected input in the left ventricular blood pool [Cin(t)] due to transit through a conduit large vessel (LV) upstream of the myocardial tissue region of interest. All the flow (F) entering the region of interest passes through the single large vessel. The large vessel has the transport function hLV(t) (Fig. 1). Therefore, the input to the myocardial region of interest is given by
where * denotes convolution and τ is a dummy variable for convolution integration. Within the microvascular system representing the region of interest, there are N parallel microvascular flow pathways. Each pathway is composed of a small vessel (SV) unit in series with an axially distributed capillary (cap) unit. The transport function for the ith pathway (where i = 1,N) is hSVi(t) * hcapi(t). The flow in the ith pathway is fiF, where F is the average myocardial blood flow in milliliters per minute per gram and fi is the local flow relative to F. The fraction of the organ mass having this local flow is wiΔfi, and the probability density function of relative local flows in the heart is w(f). Δfi, is the width of the ith flow interval in the discretized flow distribution, such that their weighted sum is unity, i.e.
The weighted mean of the local flows is also unity, i.e.
The input to each parallel pathway is CLV(t), and the fraction of the average flow that goes through the ith pathway is wifiΔfi. Except for the local flow, each pathway is identical. With flow in the large and small vessels independent, the model can be said to have “random coupling” of the single large vessel with the small vessels. Because flow in the small vessel and the capillary of each pathway is the same, the model can be said to have flow coupling between small vessels and capillary-tissue units.
The output [Cout(t)] is given by
The total quantity of indicator within the model region of interest [q(t)] corresponds to the residue curve, the externally detected regional content time curve, and includes the small vessels and the capillary-tissue units, whereas it excludes the large vessel. The quantity of indicator within the ith pathway [qi(t)] is
where Ri(t), the residue function for the ith pathway, is given by
Ri(t) represents the fraction of the indicator retained in the system at time t, after an impulse input to the ith pathway. Convoluting the actual input to the pathway [CLV(t)] with Ri(t) and weighting the relative flow (fi) by the fractional mass (wi) gives the content of the ith pathway [qi(t)]. The content of the entire microvascular system [q(t)] is the sum of the N individual pathways, excluding the large vessel
This definition of the microvascular system assumes that the indicator signal is detected only within the region of interest from which it originates, i.e., there is perfect spatial resolution in the imaging system, and that there are no contributions to the signal from large vessels or left ventricular cavity. Assumptions of the model are that the system is linear and stationary. The Ri(t) values are the integrals of the input minus the integral of the output, so Eq. 6 for q(t) provides for exact mass balance.
Vascular operators (18) were used to model indicator dispersion and delay in the large conduit vessel (VLV) upstream from the region of interest and in small vessels (VSV) within the region. The vascular operators are dispersive delay lines, having a relative dispersion [RD, standard deviation (SD) divided by ] of 0.48. Each flow pathway was modeled using a small vessel operator and capillary tissue unit in series (Fig. 1). The capillaries were assumed to have a volume (Vcap) of 0.035 ml/g to approximate the capillary plasma volume. The capillary model is an axially distributed blood-tissue exchange unit (2) in which capillary permeabilities, axial diffusion, and consumption were set to zero to model an intravascular indicator. Although the model normally provides for intravascular dispersion, when the dispersion (diffusion) coefficient is set to zero, the impulse response becomes a pure delay without dispersion. The outflow mean transit time () of the large conduit vessel is given by = VLV/F, whereas for the individual pathways within the region of interest itself, i = (VSV + Vcap)/(Ffi).
To model the effects of regional heterogeneity of myocardial blood flow, flow was apportioned among the multiple parallel pathways according to a slightly right-skewed lagged normal density distribution (17). For analyzing outflow indicator-dilution curves from the whole heart, it is appropriate to use a RD of regional flows of 0.55 (20), based on measured tracer microsphere distributions (17), extrapolated according to a fractal relation between RD and tissue piece size (7).
The parameters of the model that were adjusted to fit the MRI residue function data were F, VLV, and VSV. Constraints applied for the optimizations were that the lower and upper limits of VLV were 0.01 and 0.03 ml/g and the lower limit of VSV was 0.05 ml/g. The parameter starting values used for the optimizations were F = 1.3 ml·min−1·g−1, VSV = 0.15 ml/g, and VLV = 0.01 ml/g. A nonlinear least squares optimizer routine based on sensitivity functions (12) was used for the fitting.
To evaluate the confidence in the parameter estimates, the model was used to generate solutions to which random uniform noise (mean = zero, SD = 1) was added, representative of that observed in the MRI data. The model was used to fit these pseudodata, for which the “true” values of the parameters were known, in a fashion similar to that used for fitting the MRI data. For fitting the pseudodata, the same three model parameters that were used for fitting the MRI data were optimized; however, the upper limit for VLV was increased to 0.2 ml/g. The parameter starting values for the pseudodata analysis in Figs. 9 and 10 were the same as those for fitting the MRI data, whereas for Figs. 11 and 12 the starting values were randomized about their true values.
Three chronically instrumented adult mongrel dogs, prepared as described previously (25, 33), were studied under anesthetized closed-chest conditions in a 1.5-T MRI magnet, in accordance with guidelines in “Position of the American Heart Association on Research Animal Use.” Briefly, animals were anesthetized using pentobarbital sodium (30–35 mg/kg), and after a left thoracotomy, heparin-filled catheters were inserted into the aorta for measuring arterial blood pressure, into the left ventricular cavity for measuring left ventricular pressure and for deriving heart rate and first derivative of left ventricular pressure, and into the left atrium for injections of polylysine-Gd-DTPA. In two dogs, a 1.5-cm segment of the proximal left anterior descending coronary artery was dissected and an adjustable hydraulic occluder was wrapped around the artery to produce transient underperfusion. In these two dogs, a silicone Elastomer catheter was placed within the left anterior descending coronary artery distal to the occluder to measure mean coronary arterial pressure. The instrumentation leads were exteriorized through the chest wall, the chest was closed, and the animals were allowed to recover for 5 days. For the imaging experiments, the dogs were anesthetized, ventilated, and placed into the 1.5-T magnet, while hemodynamic variables were continuously recorded. The intravascular contrast agent polylysine-Gd-DTPA was injected as a bolus into the left atrium, while a series of consecutive cardiac images were acquired at intervals of one or two cardiac cycles per image, gated to the cardiac cycle using the left ventricular pressure signal, as previously reported (33). Ten baseline images were acquired immediately before the bolus injection, and imaging was continued for 40 s. A total of six bolus injections were performed for first-pass imaging: three in one animal, two in another, and one in another. In the two animals with coronary occluders, one of the bolus injections was performed after inflation of the occluder to produce a decrease in mean intracoronary pressure to 45–55 mm Hg. The stenosis was released immediately after completion of the imaging sequence.
Nine or 10 regions of interest were drawn in myocardial tissue in the MR images, whereas 1 region was drawn in the left ventricular cavity to provide the vascular input function. In addition, a region was drawn outside the image of the animal’s body for background subtraction. The intensity-time curves for the regions of interest were prepared for analysis by subtracting the background intensity in each image from the intensity of each region of interest and by subtracting the average intensity level in each region of interest in the 10 baseline images before each bolus from each intensity-time curve. A total of 56 tissue intensity-time curves were obtained from the myocardial regions of interest during the six bolus injections. The 56 intensity-time curves were analyzed individually.
To correct for nonlinearity in the relation between MRI signal intensity and contrast agent concentration, an in vivo calibration curve was obtained in each experiment following the first-pass imaging described above. Ten to 15 min after a vascular injection of polylysine-Gd-DTPA, when blood concentrations had reached a steady state, blood samples were drawn for measurement of the concentration of polylysine-Gd-DTPA (inductively coupled plasma atomic emission spectrometry) simultaneously with measurement of steady-state MRI signal intensities in the regions of interest (33). This procedure was repeated 7–10 times after a series of increasing contrast agent injections to ensure that the calibration curve included points at least as high as the peaks of the blood pool curves during the first-pass measurements. The calibration curves were used to transform the background- and baseline-corrected intensity-time curves to content-time curves for analysis with use of the model.
The MRI perfusion images were acquired using a 1.5-T Magnetom SP imaging system (Siemens, Erlangen, Germany) employing a Turbo-FLASH sequence with use of a preparatory nonselective radiofrequency inversion pulse and a Helmholtz surface coil positioned on the chest near the heart, as previously reported (33). The Turbo-FLASH gradient-echo pulse sequence was heavily weighted for T1 (time constant for spin-lattice relaxation) effects, such that signal changes due to T2* (time constant for transverse magnetization decay) effects were minimized. The delay between the inversion pulse and image acquisition was adjusted to minimize the myocardial signal and provide optimal sensitivity to the signal intensity enhancement caused by T1 reduction due to polylysine-Gd-DTPA. The matrix size was 60–90 × 128 with a field of view of 250 mm, resulting in a voxel size of 2–3 × 2 × 10 mm on the basis of a 10-mm slice thickness. The number of encoding steps was adjusted between 60 and 90 to ensure that the total acquisition time (354–531 ms) lasted less than one cardiac cycle. The images were gated to the cardiac cycle, and in two experiments one image was acquired every heartbeat, whereas in one experiment, where the heart rate was elevated, one image was acquired every two heartbeats. Cardiac motion appears frozen, because the images are primarily determined by the 10–30 central Fourier lines, acquired in 60–180 ms, that might be considered the “effective” image acquisition time.
Poly-l-lysine-gadolinium diethylenetriaminopentaacetic acid (polylysine-Gd-DTPA) is a contrast medium consisting of poly-l-lysine covalently linked to moieties of gadopentetate (Gd-DTPA). Polylysine-Gd-DTPA has a mean molecular mass of 52,300 Da, but because of dispersion, ~15% of the preparation exists in forms with molecular mass >75,000 and <30,000 Da (26). After intravenous injection of polylysine-Gd-DTPA, ~3% of the material was distributed into the interstitial space after 1 min (29). In the present study, there were no detectable hemodynamic effects due to the bolus injections of the isosmolar solutions of polylysine-Gd-DTPA.
The human studies were conducted for clinical diagnostic purposes in a 1.5-T Magnetom SP system similar to that described above with use of a similar gradient-echo pulse sequence to obtain one scan every heartbeat. The extracellular paramagnetic contrast agent gadolinium diethylenetriaminepentaacetic acid (Gd-DTPA, Schering, Berlin, Germany) was used in the human studies, rather than the intravascular agent polylysine Gd-DTPA, which is not approved for human use.
The patient protocol was approved by the Ethics Committee of the hospital, and all study subjects (n = 8) gave written informed consent for participation in the study. An intravenous catheter was positioned under local anesthetic into the superior vena cava, and electrocardiogram (ECG) leads were connected for gating the imaging to the cardiac cycle. Arterial blood pressure was measured using an inflatable cuff. The subjects were positioned supine in the magnet with the left ventricle at the center of an 18-cm Helmholtz surface coil. The Gd-DTPA contrast agent (0.02 mmol/kg) was preloaded into the intravenous catheter, recordings of arterial blood pressure and the ECG were started, the Turbo-FLASH sequence was initiated, and after 10 precontrast images an automatic power injector rapidly flushed (<1.5 s) the Gd-DTPA with 10 ml of saline solution into the superior vena cava while imaging continued. A total of 64 systolic short-axis images (matrix = 90 × 128, voxel size = 1.8 × 2.7 × 15 mm) were acquired at a rate of one image per heartbeat. In seven patients, a second injection of Gd-DTPA was performed after an interval of ≥20 min. Bolus injection of Gd-DTPA caused no detectable changes in arterial blood pressure, heart rate, or ECG and produced no clinical evidence of side effects.
Dispersion of the vascular input for the model was expressed using the SD of the vascular input (SDinput) of the left ventricular blood pool curves in the human and the canine data. The first-pass intensity-time curves of the background- and baseline-subtracted blood pool regions were fitted using a lagged normal density function by simultaneously adjusting the mean time (), area, RD, and skewness with use of the optimizer routine. The recirculation portion of the curves was excluded. SDinput was calculated as RD × .
There were no detectable hemodynamic changes in the dogs due to the left atrial bolus injections of polylysine-Gd-DTPA during acquisition of the MR image sequences (Fig. 2, Table 1). This finding provides evidence of stationarity, which is required for quantitative analysis.
MRI after left atrial bolus injection of polylysine-Gd-DTPA was used to determine the time course of the intravascular indicator during its transit through the left ventricular blood pool and myocardial tissue regions (Fig. 3). The Turbo-FLASH imaging sequence and cardiac cycle gating resulted in the freezing of cardiac motion in the images. Intensity-time curves were obtained from multiple regions of interest from the complete image series and, after background and baseline correction, were transformed to content-time curves by use of in vivo calibration measurements in the same animals (see METHODS).
Residue curves for the left ventricular blood pool showed an initial peak that nearly returned to the baseline during the first pass of the contrast agent followed by a smaller broader peak due to recirculation (Fig. 4). Peaks of the content-time curves for myocardial tissue were ~10% as high as the peaks of blood pool curves. This is appropriate for an intravascular indicator, for which the area under the tissue curve divided by the area under the blood pool curve equals the blood content of the tissue region (6). The tissue curves showed that the blood pool input waveform was delayed and dispersed because of indicator transport through the coronary vascular system. The peak of the curve from a region with decreased blood flow due to the coronary stenosis (Fig. 4, anterior wall) was reduced and delayed compared with the curve from a normally perfused region (Fig. 4, posterior wall). The reduction in height indicates a reduced blood volume, whereas the delay indicates that, despite the reduction in blood volume, there was a large reduction in flow.
An optimized model fit to a representative MRI tissue indicator content-time curve (Fig. 5A) provided estimated values of 0.77 ml·min−1·g−1 for regional blood flow, 0.18 ml/g for VSV, and 0.01 ml/g for VLV. To illustrate the effect of flow on the model solution, flow was increased (dashed curve) or decreased (dotted curve) by 60% from the optimized value (solid curve), with all other parameters kept constant. When flow was increased, the tissue curve increased above baseline earlier and more rapidly than in the optimized solution because of earlier arrival of the indicator. In addition, the curve with increased flow had a higher peak because of a greater inflow of indicator material, and the curve decreased more rapidly because of more rapid washout. Converse effects were observed in the solution with decreased flow.
The effect of VSV on the model solution is depicted in a similar fashion (Fig. 5B). When VSV was increased (dotted curve) above the optimized value, with flow and VLV held constant, the model solution showed the upslope identical to the optimized fit, because the arrival time of the indicator was unchanged. However, the curve had a higher peak because of a greater amount of indicator in the tissue, and the curve descended later as a result of delayed clearance. This behavior illustrates the general principle that the area under the tissue curve is proportional to the regional blood volume (6, 32). When blood volume was decreased, opposite effects were observed, whereas the upslope remained unchanged.
The effect of the third parameter that was optimized simultaneously, VLV, was that of a delay with dispersion when flow and VSV were held constant (Fig. 5C). There was no effect on the height of the curve, because VLV was not included in the residue function (Fig. 1). Although the effects of all three parameters show distinct features on the waveform of the solution, flow and VLV affect the upslope of the residue curves, indicating that these two parameters are not completely independent. Interdependence of flow and VLV, which causes increased variability in their estimates, is considered in more detail below.
Effects of blood flow and vascular volumes on the residue curve can be systematically compared using the sensitivity functions of the model parameters (Fig. 6), defined as the percent change in the model solution due to a percent change in the parameter value. Each parameter has its own sensitivity function, expressing the relative influence of that parameter on the model solution at every point in time (12). Nonzero sensitivity values (positive or negative) are necessary for parameter optimization. Because sensitivity functions are not invariant but depend on the parameter values, sensitivity functions were evaluated at flow rates of 0.2, 1.1, and 3.2 ml·min−1·g−1, with other parameters kept constant. For simplicity, the model input was a lagged normal density function that approximated the waveform of measured blood pool curves but was free from noise and recirculation.
The sensitivity functions for flow (Fig. 6B) rose abruptly to a peak at the earliest arrival of indicator in the region of interest (Fig. 6A), then began to decrease and became negative during the downslope of the residue curve. Thus an increase in flow caused an increase in the model solution during the early portion of the residue curve and a decrease during later portions (Fig. 5A). The sensitivity functions for VSV (Fig. 6C) within the region of interest remained near zero until the peak of the residue curve, because VSV has no effect on the indicator arrival time, which is influenced only by the form of the input, by flow, and by VLV, which is not a component of the residue function. An increase in VSV caused an increase in the height of the residue function. The sensitivity for VLV (Fig. 6D) upstream from the region of interest was initially negative, because it delayed indicator arrival, and then rose toward zero, having little effect after the peak of the residue curve.
The sensitivity to the level of flow was reduced when the flow rate was high, i.e., 3.2 ml·min−1·g−1, and increased when the flow rate was lowered to 0.2 ml·min−1·g−1 (Fig. 6B). The reason for the effect is that the higher the flow rate, the more nearly the waveform of the residue curve approximated that of the input function. If the waveforms of the input and solution are identical, then flow cannot be estimated independent of vascular volume. Increased flow also abbreviated the sensitivity function for flow, increasing the influence of random signal noise. Both effects indicate that flow estimates may be more reliable when flow is in the normal resting range than when it is elevated.
If one of the sensitivity functions is linearly related to one of the others or to a linear combination of the other two, then the three parameters cannot be independently determined (8). The estimation of flow is susceptible to error, because in the early portions of the residue curves there is a reciprocal relation between the sensitivity functions of VLV and flow. At later times the sensitivities for VLV and flow diverged, providing the necessary degree of independence.
Additional information on parameter identifiability was obtained from the covariance matrix during parameter optimization in the fitting of the MRI tissue curves. Correlations between the model parameters in fitting all the MRI curves averaged 0.33 between flow and VLV, −0.083 between flow and VSV, and −0.15 between VLV and VSV. These results confirm the outcome of the sensitivity function analysis and show that, in fitting the experimental MRI data, the main limitation in parameter identification was the correlation between flow and VLV.
The RDs of the large and small vessel operators and the volume of the capillary region (Vcap) were held constant because they had little effect on the model solution. For example, when the RD of the small vessel operator was decreased from the usual value of 0.48 to 0.15 in a representative solution, the differences between the two solutions over time averaged 1.4%, with a maximum difference of 6.1%. When the RD of the large vessel operator was decreased from 0.48 to 0.15, the differences over time averaged 1.7% with a maximum of 7.9%. When Vcap was increased by 50%, there was at most a difference of 7.9% in the two solutions, whereas the difference averaged 3.6%. When the two RDs were constrained to 0.15 instead of 0.48 during optimization, the myocardial flow estimates decreased by an average of 8%.
To assess the importance of regional flow heterogeneity on the residue curves, an optimized solution for a MRI tissue content-time curve was altered by decreasing the flow heterogeneity by narrowing the distribution of flows in the parallel model pathways (Fig. 7). Figure 7, inset, shows the flow distribution used for normal heterogeneity and the narrow distribution used for nearly uniform flow. The effect of decreasing flow heterogeneity was that the peak of the tissue residue curve was higher and the downslope was more rapid than for normal heterogeneity (Fig. 7). The reason for the effect is that, with heterogeneous flow (RD = 0.55), a portion of the indicator carried in the highest flow pathways traverses the system and exits before all the tracer has entered. With a more uniform flow distribution (RD = 0.055), less indicator exits the system before all of it has entered, resulting in a higher sharper peak. This finding illustrates a problem described by other investigators in using the height-over-area method for estimating 1/ to calculate F = V/ with use of the central volume principle (6, 13), because it is necessary to assume that all the indicator passing through a region is present within it for at least one sampling interval. Flow heterogeneity makes it unlikely that this condition can be satisfied for an intravascular indicator.
To investigate the effect of signal noise on the reliability of model estimates of flow and vascular volumes, a lagged normal density function without recirculation was used as a vascular input for the model to generate pseudodata: model residue functions to which uniform random noise with a mean of zero was added. Model solutions were obtained at a frequency of 1.25 samples/s to simulate the MRI frequency, and the level of added noise was at least as great as that encountered in most of the MRI content-time curves. Model fits to the noisy pseudodata were obtained using optimization procedures that were identical to those used for fitting the experimental MRI data, except the upper limit for VLV was increased from 0.03 to 0.2 ml/g to obtain more general results (see METHODS). Figure 8 shows one set of pseudodata (noise-free solutions), noisy pseudodata, and two starting solutions representing the starting values of the parameters used by the optimizer to fit the noisy pseudodata. Both sets of starting values converged to the identical final solution, which differs somewhat from the true values used to generate the pseudodata (Table 2) because of effects of the added noise. Similar results were obtained when the pseudodata were generated using an experimentally measured left ventricular blood pool signal including recirculation as the model input (results not shown).
To determine the effects of noise more systematically, pseudodata residue functions were generated with the model by use of 15 different combinations of blood flow and vascular volume spanning the physiological range with use of normal flow heterogeneity in the model (Fig. 9). The true values of flow and VSV used to generate the pseudodata are represented by the 15 intersections of the lines in Fig. 9. The true value of VLV was 0.028 ml/g for all cases. Ten different realizations of random uniform noise were added to each set of pseudodata, and the 150 sets of noisy pseudodata were independently fitted by optimizing three parameters: flow, VSV, and VLV (see METHODS for starting values). Figure 9 shows the final individual estimates of flow and VSV when the model used for the fitting included normal flow heterogeneity (heterogeneous flow model).
Parameter estimates were most reliable in a large central region of Fig. 9, where values of flow and volume were within the normal physiological range. Variability was greatest when the flow rate was high and the volume was low, as predicted by the sensitivity function analysis above. There was no evidence of systematic error in the mean estimates under the conditions examined, except for the physiologically improbable case of the lowest flow and highest volume, where there appeared to be a bimodal distribution of the volume estimates. This case also provided the worst estimate of VLV (0.0175 ml/g, true value = 0.028 ml/g). In the other 14 cases, the mean estimate of VLV was 0.0271 ml/g (range 0.0257–0.0296 ml/g).
Selected sets of the identical noisy pseudodata were also fit using a model describing completely homogeneous flow (Fig. 9) to assess the error if flow heterogeneity is neglected. In all cases examined, neglecting heterogeneity caused systematic underestimates of the correct values of flow by as much as 25% and, in most cases, also caused underestimates of VSV by lesser amounts.
The model estimates of flow and volumes were also tested for accuracy as a function of increasing dispersion of the vascular input. Pseudodata were generated while the SD of the lagged normal density function used as the vascular input (SDinput) was increased from 0.20 to 20 s (Fig. 10A) while the determinants of the microvascular system were held constant (flow, volumes, heterogeneity). Uniform random noise was added to the solutions as above, and the noisy pseudodata with their respective dispersive inputs were used for optimizing flow, VSV, and VLV. Each set of pseudodata was independently fit 10 times by use of 10 realizations of uniform random noise.
The estimates of flow showed increasing variability as the input dispersion (SDinput) was increased (Fig. 10B). The abscissa shows SDinput normalized to the SD of the residue function of the microvascular system for an impulse input (SDsystem), which was constant because flow was constant. Despite increased variability in the estimates, which were only minor for VSV, there was no systematic error in estimating the true values. The true value of VLV was 0.028 ml/g, whereas the estimates ranged from 0.027 ± 0.0022 to 0.033 ± 0.046 ml/g for the narrowest and broadest inputs. The relatively large variability was not due to negative values (estimates were constrained between 0.01 and 0.2 ml/g), but rather to right-skewed distributions of the estimates of VLV.
In the MRI dog experiments, indicator bolus injections were made into the left atrium, and the value of SDinput averaged 3.06 ± 0.528 s (SD, n = 6). Values of SDsystem were dependent on the flow, vascular volumes, and flow heterogeneity. Table 3 reports values of SDsystem and SDinput/SDsystem (left atrial injection) at a variety of flows for representative values of vascular volumes and flow heterogeneity with SDinput of 3.06 s. Figure 10B shows that the expected coefficient of variation (SD/mean) of the flow estimates was >10% only when SDinput/SDsystem was >2. The results in Table 3 show that SDinput/SDsystem was <2 after left atrial injections if the flow was <5 ml·min−1·g−1. This finding indicates that dispersion of the input for left atrial injections in the dog is not a major source of error in the estimates of regional flow. Results for intravenous injections in human subjects (intravenous injection) are discussed below.
The accuracy of the model estimates was also tested as a function of the sampling frequency of the pseudodata in the presence of a constant realistic level of random uniform noise. Pseudodata residue function curves were generated with a sampling frequency, or temporal resolution, ranging from 0.25 to 10 samples/s. Two true values of flow were used to generate noisy pseudodata (1.0 and 2.0 ml·min−1·g−1) while true values of VSV (0.11 ml/g) and VLV (0.028 ml/g) were held constant. Flow, VSV, and VLV were simultaneously estimated by fitting the pseudodata by use of 75 realizations of noise for each sampling frequency, as described above (Fig. 11). For both flow rates, there was a marked increase in the variability of the flow estimates for sampling frequencies of <2 samples/s and evidence of systematic overestimates of flow for frequencies below ~1 sample/s (Fig. 11A). Estimates of VSV showed less variability, because of decreased sampling frequency, and no evidence of systematic error (Fig. 11B). Estimates of VSV are shown for only one flow, because the effect of flow was negligible. There was also increased variability in estimates of VLV at sampling frequencies of <2 samples/s (Fig. 11C).
Imaging frequency of the MRI data ranged between 1.14 to 1.64 images/s, depending on the heart rate of the animals. Sampling frequency in this range is suboptimal and may be expected to introduce undesired variability in the individual flow estimates, although there should be no systematic error.
To investigate more systematically the interaction of sampling frequency and noise level, pseudodata were created with a range of signal-to-noise ratios for sampling frequencies of 0.75 and 2.0 samples/s (Fig. 12). Flow estimates for the pseudodata were expressed as the ratio of the estimated to the true value, so 1.0 is a perfect estimate. For a given signal-to-noise ratio (defined as the peak of the vascular input divided by the SD of the baseline noise), uncertainty in the flow estimates was less with the higher sampling frequency, indicating that deleterious effects of random signal noise can be overcome by increased sampling frequency. In practice, the value of this approach may be limited, if the increased sampling frequency leads to decreased signal-to-noise ratio. It was also observed that low signal-to-noise ratios caused systematic over-estimates of flow, although the error was less pronounced for the higher sampling frequency. For the MRI data, the signal-to-noise ratio as defined above was ~120, indicating that if the noise is uniform random, then the accuracy of the flow estimates should be adequate for many clinical and experimental purposes.
The SDinput of the left ventricular blood pool curve after central venous injections of Gd-DTPA in human subjects averaged 4.78 ± 1.55 s (n = 15), which was 56% greater than the SDinput of 3.06 s after left atrial injections in dogs. Error (coefficient of variation) in flow estimates is expected to remain below ~10% if SDinput/SDsystem is <2, as noted above. If it is assumed that SDsystem was the same for human subjects and dogs, then it is expected that error in flow estimates would exceed 10% if flow was >2.7 ml·min−1·g−1 in humans after central venous injections (Table 3). An error of ≤10% may be achieved for flows ≤5 ml·min−1·g−1 in dogs after left atrial injections because of the narrower vascular inputs.
Accuracy of the estimates of regional myocardial blood flow and vascular volume has fundamental limits because of sources of uncertainty inherent in the analysis of residue functions of an intravascular indicator independent of technical limitations associated with specific imaging modalities. In the presence of realistic values of random measurement error, sampling frequency, input dispersion, and perfect spatial resolution, the present results indicate that the fundamental limits readily permit flow estimates with acceptable accuracy for many clinical and experimental applications over most of the physiological range with use of a spatially distributed model describing vascular architecture to analyze indicator kinetics.
Fundamental limitations to the accuracy of flow estimates were observed under three general conditions: when blood flow was near the upper limit of the physiological range (>3 ml·min−1·g−1), when dispersion of the vascular input was more than two to three times the dispersion of the microvascular system for an impulse input, and when the temporal resolution of the data was <1–2 images/s. It was assumed that there were no important limitations to the spatial resolution of the residue function data. The most important limitation in parameter identification was that flow and VLV were not completely independent.
The reason for decreased accuracy of flow estimates due to increased blood flow is that the higher the flow, the more closely the waveform of the residue function approaches that of the vascular input. If the waveforms of the residue function and vascular input are alike, then the sensitivity function for flow is reduced, as shown in Fig. 6, making it difficult to obtain accurate flow estimates. An additional result of increased blood flow is that the residue curve is compressed in time, reducing the number of data points on the upslope of the residue curve, which have the highest sensitivity to flow. Because of a decreased number of data points, random error in the signal is not averaged out, increasing the uncertainty. This problem is not expected to limit accuracy of flow estimates in the normal range of coronary flow in the unstressed heart.
Dispersion of the vascular input does not appear to be an important problem for bolus injections performed in the left atrium (Fig. 10) or for central venous injections (Table 3), indicating that the present approach has potential for clinical application. However, increased input dispersion may seriously limit accuracy for peripheral venous injections of indicators. The error due to increased input dispersion is partially offset by obtaining a larger number of data points along the upslope of the curve, averaging out signal noise in the portion of the data most sensitive to flow. There does not appear to be significant error due to indicator recirculation, because this appears late in the curve and is modeled as a delayed component of the vascular input.
An important source of error in flow estimates is the limited number of data points on the rising portion of the residue curve, the portion most sensitive to flow (Fig. 11). With few points, any random signal error will reduce the accuracy of the fitted model solution. Overly long image acquisition intervals will be particularly troublesome with high blood flows, when short arrival times also decrease the number of data points on the upslope of the curve. These results indicate that flow estimates may be improved by increasing the MRI frequency to ≥2 images/s centered on the times of least cardiac motion near end systole and end diastole. Whereas there would be difficulty in positioning the regions of interest in the pair of image sequences over the same myocardial tissue regions, the paired residue curves could be analyzed simultaneously for flow and volumes to account for changes in myocardial blood volume (VSV) during the cardiac cycle.
Estimates of VSV showed less variability than flow under all conditions. The probable explanation is that the sensitivity functions for VSV were high throughout the residue curves and were distinct from the sensitivity functions for flow and VLV (Fig. 6). Thus it appears that vascular volume within a tissue region of interest may be more accurately estimated than flow (M. Jerosch-Herold and N. Wilke, unpublished observations).
The importance of including vascular architecture in the model was underscored by the finding that neglecting flow heterogeneity in the analysis caused systematic underestimates of regional flow by as much as 25% (Fig. 9). The representation of microvascular architecture in the model was based in part on measurements of regional myocardial flow heterogeneity by use of tracer microspheres in baboons (17) and sheep (5). Indicator kinetics reflect dispersion at every level of the vascular system, including capillaries, but measurements of microsphere deposition are limited to macroscopic tissue pieces. Regional flow heterogeneity was observed to increase as the size of the tissue pieces examined was decreased, suggesting a fractal relation (3). The flow distribution in the multiple pathways of the model is important in general for avoiding bias in parameter estimation, and for the present analysis the fractal relation was used to extrapolate the degree of heterogeneity to that expected at the capillary level (1, 7).
There is some uncertainty as to whether a RD of flow of 0.55 represents the appropriate heterogeneity to be used for modeling residue function kinetics. It was proposed that the degree of heterogeneity can be calculated from the ratio of the weight of the tissue region of interest to the weight of the whole heart (4). If it is assumed that a RD of 0.55 is correct for a whole heart weighing 150 g and that the fractal dimension (D) is 1.2, then a better estimate for the RD of a 2-g region of interest may be 0.55 × (2/150)(D−1), or 0.23. However, the RD of hamster hearts weighing 0.5 g averaged 0.31, on the basis of quantitative autoradiography (28), from which it may be calculated that a more appropriate value of RD would be 0.31 × (2/0.5)(D−1), or 0.41. Therefore, it is possible that the RD of 0.55 used in the present study exaggerated somewhat the heterogeneity within realistic imaging regions of interest. In particular, it is possible that under conditions of increased blood flow the RD should be reduced to values <0.55, because there is evidence that flow heterogeneity is reduced during increased coronary flow (14). The effects of altered microvascular architecture on flow estimates on the basis of residue functions have been considered elsewhere (13).
Distributed modeling used in the present study overcomes limitations involved in the use of the central volume principle to estimate flow by analyzing mean transit times and area over height of the tissue curves (6, 13, 32). The major limitation in the analysis of residue functions with use of the central volume principle is that it requires that all the indicator that will pass through a region of interest is present within the region for at least one sampling interval. Because of the effects of heterogeneous blood flow and input dispersion, this condition will rarely be achieved in vivo. An advantage of distributed modeling is that this requirement is avoided. Another advantage of distributed modeling of myocardial residue functions is that indicator kinetics measured in the left ventricular blood pool are used as the input function for the model (30).
A recent advance in the analysis of indicator kinetics is to deconvolve the vascular input curve with the tissue residue curve, from which the regional mean transit time for an impulse input may be calculated (13). The approach may provide accurate estimates of regional blood flow if blood volume is also determined, although systematic error in the transit time estimates was observed if there was dispersion between the site at which the vascular input was measured and the region of interest.
The results in the present study may underestimate the error in regional flow estimates obtained by analyzing residue functions obtained using real imaging systems. Error will be greater if the image noise has a structure different from the random uniform additive noise with a mean of zero used in the pseudodata. For example, MR image-triggering errors due to cardiac arrhythmias would cause nonrandom noise. Another possible cause of increased error in flow estimates is that whereas the model is a perfect representation of the pseudodata, it is only an approximation for the physiological processes of microvascular transport in vivo. For example, the vascular architecture described in the model is a simplification of the real microvascular anatomy, although it has been successfully used for modeling indicator-dilution curves in the heart (20). In addition, the model does not describe the effects of conduit vessels passing through the region of interest under observation supplying capillary beds in neighboring regions. Such nonnutritive flow may affect residue function kinetics, particularly if large epicardial coronary vessels are included in a region of interest.
MRI may be particularly well suited for studying regional flow with use of an intravascular indicator, because the high spatial and temporal resolution of MRI makes it possible to analyze first-pass kinetics. Another advantage is that the low concentrations of the contrast agent do not cause transient disturbances in hemodynamic conditions that may be observed using hyperosmolar ionic contrast media (19). A disadvantage is the nonlinearity of the relation between MR image intensity and indicator concentration.
Improved MRI systems under development, with T1-weighted image acquisition times as short as 100–200 ms, may allow imaging frequencies three times faster than those used in the present experiments. Signal-to-noise characteristics may be improved by the use of phase-array radiofrequency coils, instead of Helmholtz coils (21). In addition, MRI pulse sequences with driven equilibrium preparation (16) to decrease signal variations due to changes in heart rate and cardiac arrhythmias may reduce image artifacts and noise. If such technical innovations make it possible to increase the MRI frequency to 2 samples/s without sacrificing (and perhaps improving) the image signal-to-noise ratio, then it will be possible to increase the accuracy of regional flow estimates.
Results of the present study show that accurate estimates of myocardial blood flow and vascular volume may be obtained using a spatially distributed model describing vascular architecture to analyze noninvasively detected residue function first-pass kinetics after bolus injection of an intravascular indicator into the central circulation. Increased variability in regional flow estimates is expected when flow is near maximal physiological values, when dispersion of the vascular input is more than twice the dispersion of the microvascular system, and when the sampling frequency of the imaging system is <1–2 samples/s.
We are indebted to Kamil Ugurbil (Center for Magnetic Resonance Research, University of Minnesota) for critical suggestions and encouraging support during frequent discussions. We also thank Schering (Berlin, Germany) for providing the contrast agent and for carrying out the concentration measurements.
This study was supported by National Institutes of Health Grants HL-50238 and RR-01243.