Abstract
We have addressed 2 major challenges of 82Rb cardiac PET, noninvasive estimation of an accurate input function and absolute quantitation of myocardial perfusion, using a generalized form of least-squares factor analysis of dynamic sequences (GFADS) and a novel compartment analysis approach. Methods: Left and right ventricular (LV+RV) time–activity curves (TACs) were generated from 10 rest/stress studies, and 30 myocardial TACs were modeled to cover a range of clinical values. Two-dimensional PET Monte Carlo simulations of the LV, RV, myocardium, and other organs were generated separately and combined using the above TACs to form 30 realistic dynamic 82Rb studies. LV and RV TACs were estimated by GFADS and used as input to a 2-compartment kinetic analysis that estimates parametric maps of myocardial tissue extraction (k1) and egress (k2), as well as LV+RV contributions (fv, rv), by orthogonal voxel grouping. In addition, 13 patients were injected with 2.22 ± 0.19 GBq (60 ± 5 mCi) of 82Rb and imaged dynamically for 6 min at rest and during dipyridamole stress. Results: In Monte Carlo simulations, GFADS yielded estimates of the 3 factors and corresponding factor images, with average errors of −4.2% ± 6.3%, 3.5% ± 4.3%, and 2.0% ± 5.5% in the LV, RV, and myocardial factor estimates, respectively. The estimates were significantly more accurate and robust to noise than those obtained using TACs based on manually drawn volumes of interest (P < 0.01). The 2-compartment approach yielded accurate k1, k2, fv, and rv parametric maps; the average error of estimates of k1 was 6.8% ± 3.6%. In all patient studies, our approach yielded robust estimates of k1, k2, fv, and rv, which correlated very well with the status of the subject and the catheterization results. Conclusion: Quantitative dynamic 82Rb PET using generalized factor analysis of dynamic sequences and compartmental modeling yields estimates of parameters of absolute myocardial perfusion and kinetics with errors of <9%.
Cardiac PET with 82Rb allows the assessment of absolute myocardial perfusion, as well as coronary flow reserve (CFR), using a column generator (1–7) in clinics that lack a cyclotron or other infrastructural support required to produce radionuclides such as 13N-ammonia. The short half-life of the radiotracer (76 s) makes possible rapid rest/stress paired studies within a very short time (∼15 min), allowing rest and stress imaging under virtually identical conditions and decreasing the total time required to scan each patient. In addition to the poor spatial resolution due to the relatively long positron range of the 82Rb positron emitter, the major challenges of parametric 82Rb imaging are estimation of an accurate input function from noisy data without arterial blood sampling and the absolute quantitation of myocardial perfusion. One approach, proposed by Lin et al. (8,9), to deal with noise is to filter the reconstructed volume before region-of-interest quantitation using wavelet transforms. Compartmental modeling for 82Rb cardiac PET has also been investigated (3,8,10), and 2 models have been proposed by Gould (10): a model that reduces the 2 compartments to one unknown, flow, which is calculated from myocardial uptake and arterial input function, and a compartmental model that fits the model equations to observed myocardial and blood time–activity curves (TACs) to estimate the 2 unknowns. The physiologic model of 82Rb kinetics in the myocardium has also been compared with reduced-order models by Coxson et al. (11).
In this study, we used a generalized form of the least-squares factor analysis of dynamic sequences (FADS) to obtain a robust estimate of left and right ventricular (LV+RV) input functions automatically, without the need to draw volumes of interest (VOIs), and developed a compartment model based on orthogonal grouping to estimate on a voxel-by-voxel basis, rather than within selected VOIs, extraction of 82Rb in the myocardium. We validated our approach by Monte Carlo simulation studies and demonstrated clinical feasibility by patient studies.
MATERIALS AND METHODS
Generalized Factor Analysis of Dynamic Sequences (GFADS)
Background.
FADS is a powerful technique for the analysis of dynamic sequences; its major drawback is that unique solutions are not ensured (12). Several techniques have been developed that address this problem and improve the results of factor analysis (13–17). These techniques are based on the use of a priori physiologic information. Therefore, they are tailored for a particular type of clinical study and cannot be used without modification in different settings. For example, a particular factor analysis approach might yield satisfactory results for studies of healthy people but might not work for studies of patients with different degrees of ischemia, different ejection fractions, and, therefore, large differences in the shape and amplitude of blood TACs. The analysis in this article addresses that possible concern by analyzing dynamic cardiac 82Rb PET studies under a wide range of ejection fractions and blood TACs.
Furthermore, although these techniques increase the range of situations in which unique FADS solutions are achieved, they do not ensure a unique solution in all cases. We have previously developed a technique that is more general than previously reported approaches and can be used in a variety of applications, with a wide range of ejection fractions and blood TACs (18). In this work, we address the nonuniqueness problem by penalizing its main effect, spatial overlap between factor images. In our approach, FADS is performed first by any factor analysis method such as the apex-seeking (19) or least-squares (17) approaches. After FADS, a second step is used to minimize the overlap between factor images and, hence, increase the probability of a unique solution. This second step is described and tested in Monte Carlo simulations of realistic dynamic cardiac 82Rb PET studies, and the feasibility of the approach is demonstrated in patients.
Theory.
A dynamic sequence of medical images can be represented by an N × M matrix A. N is the number of voxels in one dynamic image and M is the number of time frames. The factor model of the dynamic data assumes that the data matrix can be represented by the following equation:
Eq. 1
where F contains P factors and n denotes noise in the data. The factor curves define the time course of a given factor whose spatial definition is contained in matrix C (the factor image). To solve Equation 1, the number of factors P must be known a priori. In this work, a value of 3 was chosen; the factor represents the blood activity in the left and right ventricles and myocardial tissue, as described in detail in the compartment analysis section. Estimation of the factors (F) and factor images (C) was based on minimization of the least squares objective function fLS:
Eq. 2
Nonnegativity constraints were applied by adding another penalty term, fn(C,F), to the objective function (Eq. 3). fn(C,F) is a quadratic function that heavily penalizes negative values of C and F and yields a positive solution (17):
Eq. 3
where:
Eq. 4
Therefore, the total objective function that was used in the first step was:
Eq. 5
where a is a penalty coefficient that can be chosen from a wide range of large values (17). In this work, the entries of A were normalized to the range [0,1] and a was set to 3,000. The total objective function was minimized using the conjugate gradient algorithm (20).
Because the factor model described by Equation 1 is not mathematically unique, the result of the optimization described by Equation 5 cannot be used for quantitative analysis. To address this problem, we generalized our previous approach (18) to allow implementation as a postprocessing method following any factor model technique. In this approach, GFADS, we modify the factors and factor images obtained in the first step by minimizing the nonnegative term funi(R) that penalizes the overlap between images of factor coefficients while keeping CF constant:
Eq. 6
The result of the minimization is C′ = CR, where R is the rotation matrix determined during the optimization, so that the results obtained in the first step hold during optimization in the second step:
Eq. 7
Hence, the problem of nonuniqueness of the FADS solution is addressed in our GFADS approach by minimizing the major effect of nonuniqueness, overlap between the factor images. This is accomplished by optimizing the rotation matrix R by a minimization of funi under nonnegativity constraints imposed by fn(CR,R−1F). The total objective function that was used in the second step (funi + bfn) was minimized by a simplex algorithm (20), and convergence was obtained after 100–500 iterations. The value of the penalty parameter, b, 0.1 in this work, was not found to be critical for the convergence to a solution when chosen in the range [0.05, 0.5].
VOI Analysis
LV, RV, and myocardial TACs were also estimated using a VOI-based approach. The initial VOIs used to define the activity distributions in the LV, RV, or myocardium for the Monte Carlo simulations were also used for estimating the TACs. These volumes were first reduced by one eighth (50% in each dimension) to mimic usual clinical practice to minimize spillover among LV, RV, and myocardial tissues. For the patient studies, VOIs were drawn by a skilled operator over several slices of the 3-dimensional (3D) volume of the LV, RV, or myocardium obtained by summing the first 24 dynamic frames, and counts within these volumes were estimated at each frame to generate the corresponding VOI TAC.
Compartment Analysis
Because 82Rb is only partially extracted by the myocardium, a 2-compartment model is required for accurate estimation of myocardial extraction fraction (10). The 2 compartments of the model are the “free rubidium space” (blood perfusing the myocardium plus interstitial space) and the “trapped rubidium space” (muscle of the myocardium). The main parameters of the model are the kinetic transport constants k1 (mL/min/g) and k2 (min−1), which denote the extraction (forward) and egress (backward) rates of transport between the metabolically trapped space (myocardium) and the freely diffusible space (blood pool), respectively.
Kinetic Model.
The TAC in each voxel was modeled as a combination of 3 contributions: the contribution from myocardial tissue, modeled using a 2-compartment model, and contributions from LV and RV cavities, modeled as fractions of measured LV and RV functions. In cardiac imaging it is necessary to explicitly model the contribution from the RV to obtain stable and robust voxel parameter estimates (10). This model can be expressed as:
Eq. 8
where TACi(t) is the value of the voxel i at time t, TACi is the TAC of voxel i, and I(t) is the input function (i.e., measured LV function). k1i, k2i, fvi, and rvi are the kinetic parameters for voxel i. Assuming a soft-tissue density similar to that of water, k1i (mL/min/g) characterizes myocardial tissue extraction (inflow). k2i (min−1) characterizes myocardial tissue egress (outflow), fvi (dimensionless) represents the contribution to the total voxel activity from the blood input function I(t), and rvi (dimensionless) represents the contribution from the activity in the RV, R(t), which in general differs from the input function. Both I(t) and R(t) were determined by GFADS and assumed to be known. Note that k1i, k2i, fvi, and rvi are zero-valued for some voxels and that, unlike the widely used (1 − fv) formulation, our approach allows fitting the model to all image voxels, including voxels where fv = 1 or rv = 1.
Orthogonal Grouping.
The estimation of parametric images on a voxel-by-voxel basis yields very noisy images due to high levels of noise in TACs derived from single voxels. Using a filtering approach reduces noise but can result in a degradation of spatial resolution of the parametric images. In this work we used a new method, orthogonal grouping, which consisted of a 3D grouping of voxels with similar time course into several groups of size G. The average TACs within the groups were used for calculation of the kinetic parameters. This approach leads to reduced levels of noise with minimal effects on spatial resolution. The orthogonal grouping consisted of 2 steps:
Step 1. An orthogonal analysis was performed on TACi(t). The reduced-dimensional data subspace S that contained most of the data was determined as the space spanned by principal vectors by setting small singular values to zero. In the case of cardiac 82Rb scans, we considered a 3D subspace. The principal vectors spanning subspace S do not correspond to the physiologic tissue curves, and some have negative values. The coordinates of each dynamic voxel in the 3D subspace were calculated as the product of the principal vectors with TACi(t)—that is, we obtained 3D vectors P(TACi(t))τ, where τ = 1… 3 that described the location of a given dynamic voxel TACi(t) in the 3D data subspace S, and where P denotes the projection operator onto the τth base vector.
Step 2. The grouping algorithm was applied to the 3D vectors P(TACi(t))τ. First, Euclidian distances dij =
P(TACi) − P(TACj) were calculated between each 2 dynamic voxels. The dynamic voxel for which the sum of distances to all other voxels was highest was determined. Next, the G − 1 closest dynamic voxels were found. These G dynamic voxels (initial 1 and G − 1) constituted the first group. The procedure was repeated until all dynamic voxels were grouped. G was chosen iteratively to be the smallest value that ensures convergence for all groups (G was 30 for the Monte Carlo simulations and 100 for the patent studies). Finally, the estimated kinetic parameters were remapped back to the original volume.
Monte Carlo Simulations of Dynamic 82Rb PET Studies
Kinetic Modeling of Dynamic 82Rb Studies.
LV and RV TACs were measured in 10 human rest and dipyridamole stress dynamic cardiac 82Rb PET studies using small VOIs positioned over the LV and RV cavities. The patients used to generate the LV and RV input functions were imaged at different time points in the life of the 82Rb generator. There were 6 subjects without evidence of obstructive coronary artery disease (CAD) as determined by the absence of regional perfusion defects on the dipyridamole stress images (summed stress score [SSS] = 0). In addition, there were 4 patients with known obstructive CAD and extensive perfusion abnormalities on the dipyridamole stress images (SSS > 18). This yielded LV and RV TACs that varied greatly in shape and represented a wide range of physiologic conditions.
Because 82Rb is a partially extracted tracer, the myocardial activity detected in a small region is the sum of the activity trapped in the cell and the free circulating activity. When modeling the myocardial uptake, we assumed that the myocardial tissue contains approximately 25% of arterial blood (10). The LV TAC, I(t), was used as input function for a 2-compartment kinetic model to generate a pure myocardial extraction fraction (extravascular space) TAC:
Eq. 9
Myocardial extraction TACs were generated for rest studies with k1 values of 0.4, 0.7, and 1.2 mL/min/g and for stress studies with k1 values of 0.7, 1.4, and 2.1 mL/min/g. Because 82Rb is an analog of potassium that is taken up rapidly and retained for a long time within myocardial cells (10), we used very low k2 values of 0.01 and 0.02 min−1 in rest and stress studies. The selected values of k1 and k2 spanned the range of values reported in the literature in rest and stress studies (10). The lowest two k1 values were associated with k2 = 0.01 min−1, whereas the highest k1 value was associated with k2 = 0.02 min−1. Typical myocardial TACs obtained using our kinetic model were compared with those measured using small myocardial VOIs in patients, and several were found to match patient TACs quite well. Since 3 pairs of (k1, k2) values were used with 10 pairs of (LV, RV) input functions, a set of 30 unique realistic simulated dynamic studies was obtained. The TACs for the liver and soft tissue were assumed to be delayed constant functions, identical to each other. Sample LV and RV input functions, as well as 3 myocardial TACs, are shown in Figure 1 (solid lines).
Simulated RV and LV input functions and associated myocardial TACs corresponding to 3 different sets of values of k1 (mL/min/g) and k2 (min−1) as well as corresponding factors estimated with GFADS (open symbols). MYO = myocardium.
Monte Carlo Simulations.
The anthropomorphic Zubal torso phantom was used (21), with the heart segmented into LV and RV cavities, myocardial muscle, and blood pool, which were assumed to have uniform activity distributions. 2-Dimensional (2D) PET Monte Carlo simulations of the LV and RV, as well as the myocardium and other organs (liver, soft tissue, etc.), were performed separately using the SimSET software (22,23) and combined to generate 30 realistic dynamic 82Rb studies using the TACs obtained in the previous section.
We modeled a 2D PET scanner similar to the GE DST Discovery PET/CT scanner (GE Healthcare) with a uniform 3-cm-thick bismuth germanate (BGO) detector (energy resolution, 20% at 511 keV; axial extent, 15 cm; ring diameter, 88 cm). A stack of annular septa (septal thickness, 0.1 cm) was explicitly modeled in the Monte Carlo simulation in the axial direction, yielding twenty-eight 3.3-mm-thick slices (voxel size, 3.3 × 3.3 × 3.3 mm3). Events reaching the detector were binned into 128 projections using only the direct planes. Nine orders of Compton scatter were modeled. Positron range was also modeled using the analytic approach developed by Palmer and Bronwell (24), which uses β-decay energy spectra and empiric range formulae.
Random coincidences, not presently modeled in SimSET, were also modeled in the simulation software that we used in this work. Random coincidences were estimated from singles using:
Eq. 10
where CTW is the coincidence timing window (12 ns), and Rsingle,1 and Rsingle,2 are the single-channel counting rates in the 2 detectors of the pair. Singles were simulated in a separate Monte Carlo simulation by tracking single 511-keV photons (rather than photon pairs) using the same activity and attenuation distributions as those used in the actual PET simulation. Twice as many decays were generated for the singles simulation than for the PET simulation because 1 decay yields 2 annihilation photons. Previously validated variance reduction techniques (stratification and forced detection) (25) were used when tracking photons.
Nine billion photons were generated when simulating each structure (LV, RV, myocardium, soft tissues). This yielded essentially noise-free sinograms that were used as the basis for generating noisy sinograms. Scatter, known from simulation, as well as random coincidences, estimated using the singles-rate formula (Eq. 10), was subtracted. The resulting sinograms were precorrected for attenuation by multiplying each ray-sum by the inverse of the corresponding attenuation factor, and reconstructed using ordered-subsets expectation maximization (OSEM, 5 iterations, 4 subsets) (26). Stationary resolution was assumed in the forward projector and no postfiltering was performed. This yielded 30 studies, each consisting of twenty-four 5-s frames (128 × 128 × 28 voxels each).
Dynamic 82Rb Patient Studies
Thirteen patients were included in this analysis. Five patients had evidence of obstructive CAD, and 8 had no evidence of CAD (Table 1). Rest and stress dynamic protocols were followed on 8 subjects, whereas rest dynamic data were available on 4 subjects and stress dynamic data were available on 1 subject. Patients undergoing the rest/stress protocol were injected with a bolus of 2.22 ± 0.19 GBq (60 ± 5 mCi) of 82Rb (Bracco Diagnostics) in 14 ± 6 mL of saline and imaged dynamically for 6 min at rest and, 10 min later, during dipyridamole stress (intravenous infusion of 0.14 mg/kg/min for 4 min). The dynamic imaging protocol consisted of twenty-four 5-s frames followed by eight 30-s frames. The initial fast frames were used to capture the rapid washin and washout of 82Rb from the LV and RV cavities.
Clinical Data and Calculated Rest/Stress Myocardial Blood Flows for 13 Patients Included in Study
Sinograms were acquired in 2D mode on a DST Discovery PET/CT scanner with BGO detectors. A CT scan was also performed with each PET study (70 mAs, 140 kVp) using an 8-slice helical scanner and shallow breathing (helical thickness, 3.75 mm; pitch, 1.35:1; 27 mm/rotation). The attenuation map used for correction of the 511-keV photon attenuation was derived from the CT scan using a continuous conversion scale with a range of slopes dependent on the CT kilovolts and the CT number (27). Randoms correction was performed by direct subtraction of delayed events and scatter correction was performed using the scatter modeling approach proposed by Bergstrom et al. (28). All dynamic sinograms were reconstructed using the attenuation-weighted OSEM (21 subsets, 2 iterations, as recommended by the manufacturer) into twenty-four 5-s and eight 30-s frames, each being a 128 × 128 × 47 volume. No postfiltering was performed. Next, LV and RV input functions were computed using GFADS, and the 2-compartment analysis was performed, yielding voxel-by-voxel parametric maps. No decay correction was performed (k1, unlike k2, is not affected by radioactive decay) to avoid introducing inaccuracies arising from the short half-life for a 5-s frame.
Data Analysis
Estimation performance was assessed for the Monte Carlo–simulated data by computing known true values. The LV and RV TACs estimated with GFADS or by VOI analysis were compared with the LV and RV input functions used to generate the dynamic studies. Likewise, the myocardial TACs estimated with GFADS and VOI were compared with the true myocardial TAC. For a given Monte Carlo simulation, the total GFADS (or VOI) estimation error εGFADS of the LV, RV, or myocardial TAC was computed over the 24 time frames, and the average error over the 30 Monte Carlo studies was expressed as:
Eq. 11
The SD of the average error was also calculated over the 30 Monte Carlo studies after averaging over the time frames. Likewise, the average errors associated with myocardial tissue extraction (k1) and egress (k2), as well as RV and LV contributions (fvi, rvi), were computed over the 30 Monte Carlo simulations. The rate constants k1, k2, as well as fvi, and rvi, were averaged over each tissue compartment using the known simulated activity distribution as a mask, and errors were calculated using the known parameters’values used to generate the dynamic studies. Means and SDs of the errors were calculated over the 30 studies.
RESULTS
Figure 1 shows examples of the LV and RV input functions and myocardial TACs that were assumed for the simulations, as well as the corresponding curves estimated using GFADS. GFADS yielded accurate estimates of all 3 factors modeled in the Monte Carlo simulations, as well as the corresponding factor images, with average errors of −4.2% ± 6.3%, 3.5% ± 4.3%, and 2.0% ± 5.5%, respectively, in the LV, RV, and myocardial factor estimates. Errors were calculated with respect to the true factors used in the simulation over all time frames. Figure 2 shows the simulated LV and RV input functions and the myocardial TACs, as well as the corresponding factors, estimated by GFADS and by the VOI TACs. Average errors were significantly higher with the VOI TACs than with GFADS (P < 0.001, paired t test) with average errors of LV, RV, and myocardial factor estimates of 30.7% ± 8.3%, 26.5% ± 6.2%, and 7.3% ± 7.1%, respectively. A single transverse image from the factor image corresponding to each factor is also shown. Note the underestimation of the LV and RV input functions at early time points and the overestimation of these input functions at later time points. Also, note the overestimation of the myocardial TAC by the VOI TAC technique at the early time points; this is a result of spillover of counts from the LV compartment into the myocardial compartment. This spillover is not observed with GFADS.
Estimated factors with VOI TACs and GFADS and GFADS-derived factor images. The 3-factor images were well separated with GFADS, whereas the LV+RV input functions were underestimated and the myocardial (MYO) TAC was overestimated with VOI TACs, due to spillover in early time points.
The 2-compartment model kinetic analysis also yielded estimates of k1, k2, fv, and rv parametric maps. The error in estimating of the value of k1, averaged over the myocardium, was 6.8% ± 3.6%. This error was slightly greater for k2 (8.4% ± 7.6%) as the simulated k2 values were very low (i.e., 0.01–0.02) and, hence, more sensitive to noise. Errors in the estimates of fvi and rvi were 6.5% ± 2.8% and 7.3% ± 2.7%, respectively. The corresponding errors when using TACs estimated using VOIs as input functions to the kinetic analysis were 28.6% ± 7.9%, 11.3% ± 6.3%, 45.7% ± 7.8%, and 29.7% ± 6.0% for k1, k2, fvi, and rvi , respectively.
Patient clinical data at rest and stress, as well as the SSS and CFR, are given in Table 1. In all patient studies,GFADS yielded LV and RV input functions, myocardial TACs, and corresponding factor images of the expected form (Fig. 3). Furthermore, the 2-compartment kinetic model yielded parametric maps of myocardial tissue extraction and egress as well as spillover fractions from the blood pool into the myocardial tissue (Fig. 4). Note that these images were obtained despite the short acquisition time of the dynamic frames (5 s each) and without postreconstruction filtering; this implies that the new approach is robust to noise. As expected, the k2 values measured in patients using the LV and RV input functions estimated with GFADS were very low and ranged between 0.012 and 0.028 min−1. The corresponding k1 values ranged between 0.52 mL/min/g in patients with suspected CAD and 1.03 mL/min/g in patients without suspected CAD at rest. The corresponding range of k1 values at stress was 0.74–1.84 mL/min/g. Furthermore, k1 values estimated in anterior, inferior, lateral, and septal walls did not differ significantly (NS, paired t test) in patients with normal myocardial blood perfusion (SSS = 0). The flow values (k1) estimated using the LV and RV input functions estimated by VOI analysis were systematically greater (8.3%–27.6%) than those obtained using GFADS LV and RV input functions (P < 0.05). k2 values associated with VOIs were slightly greater (5%–10%) than those with GFADS; however, fv (0.39 ± 0.07) and rv (0.31 ± 0.02) parameters were significantly greater (13.4%–31.7%) with VOIs as compared with GFADS (P < 0.05). The CFR was also calculated as the ratio of flow at peak stress and flow at rest (k1 stress/k1 rest). CFR values ranged from 1.58 to 2.3 in patients with no prior known CAD and from 1.39 to 1.51 in patients with known prior CAD. k1 and k2 values were not strongly correlated at rest or stress (r2 < 0.69).
Estimated factors and factor images (3 transverse slices) with GFADS for a typical dynamic 82Rb PET study in 93-y-old female with known CAD but with normal myocardial perfusion (SSS = 0). Radioactive decay was not compensated. LV, RV, and myocardium (MYO) were well separated using our uniqueness constraint.
Parametric maps of myocardial tissue extraction (k1) and egress (k2) as well as RV and LV contributions (fvi, rvi) from the patient dynamic study shown in Figure 3. Note absence of aberrant points due to orthogonal grouping.
DISCUSSION
GFADS was used to estimate noninvasively the LV and RV input functions from the dynamic 82Rb PET studies. After fitting the time-varying factor model to the dynamic data using a least-squares objective function, a different objective function that penalized spatial overlap between factor images was minimized. Note that this optimization preserves the least-squares fit obtained in the first step and that both steps incorporated nonnegativity constraints on the factors, as well as on the factor images. This approach does not require a priori knowledge of the kinetics and can be used in other dynamic imaging applications. Furthermore, our approach does not require drawing VOIs to obtain the LV and RV TAC input functions. Only a cube that represents a cropped volume containing the heart is needed. This is a major advantage as it obviates the need for manual intervention in the quantitation scheme and makes it reproducible. Furthermore, the estimates obtained with FADS were significantly more accurate than those with VOI TACs, suggesting that the spillover and tissue overlap are the sources of significant errors with the latter approach. In fact, the spillover from the LV compartment to the myocardium was greatest for VOI TACs at the lowest values of simulated blood flow, suggesting that the largest estimation errors are made when blood flow is reduced by disease. Furthermore, ideally, GFADS provides independent factor images of the LV, RV, and myocardium. Therefore, the kinetic model fitting to Equation 8 using the blood TACs derived from GFADS should, essentially, subtract the influence of LV and RV blood from the myocardial contribution, making a partial-volume correction for spillover unnecessary with our approach. Therefore, a spillover due to myocardial contamination by the input function or to cardiac motion translates into an overestimation of fvi and rvi but not of k1 or k2. Finally, the GFADS estimates were more robust to noise than those obtained using VOI quantitation. This is consistent with the fact that the VOI quantitation is based on fewer voxels compared with the factor model, for which each time point of the curve estimated with GFADS results from the fitting of the entire image of corresponding factor coefficients to the data at that time point. If GFADS were applied to the VOI data only, there would be no advantage in terms of noise, and results similar to the VOI quantitation would be obtained. Although the number of factors must be defined before performing GFADS, we have found that P = 3 always yielded robust estimates of the LV and RV input functions in the Monte Carlo simulations and patient studies. This is consistent with the fact that the first 3 eigenvectors, obtained by principal component analysis, were consistently several times greater than the other eigenvectors in both Monte Carlo simulations and patient studies. This is also consistent with previous results performed with P = 3 in teboroxime cardiac canine studies and ischemic patients (17,29). In the future, we plan to investigate the potential of using a greater numbers of factors. We also plan to determine the optimal cluster size G, although the fact that a fixed value of G = 30 was successfully used for all Monte Carlo simulations and that G = 100 was successfully used for all patients suggests that the choice of G is not crucial.
In the Monte Carlo simulations, positron range was modeled using an analytic approach (24). Modeling positron range is important in PET Monte Carlo simulations of 82Rb, because the full width at half maximum of the positron range is 2.6 mm as compared with 1 mm for 18F, so the effect of positron range on spatial resolution and image quality is greater for 82Rb than for 18F. Likewise, modeling random coincidences is of paramount importance in the case of 82Rb cardiac PET because of the high counting rates after injection of 2.22 GBq (60 mCi) of 82Rb, and the torso attenuation. Dead time is another important factor that affects singles and coincidences differently; although this effect was not modeled in this study, we are currently working on including it in the Monte Carlo simulation.
In the Monte Carlo simulations, patient-derived TACs were used to generate spatially uniform LV and RV compartments, and GFADS yielded accurate estimates of the original simulated TACs with very little overlap between LV, RV, and myocardium. In patient studies, one factor might not be sufficient to describe the dynamic activity distribution in the LV or RV when the activity does not vary uniformly within these structures. In this case, more factors would be needed to fully describe the dynamic ventricular activities. We are presently investigating the impact of patient hemodynamics on the shape of LV and RV input functions. The kinetic approach presented in this work was designed to model such situations; the grouping algorithm does not necessarily result in spatially contiguous clusters as the distance metric is not a spatial distance vector. This allows different locations within a given structure to have different kinetic parameters, a desirable property when nonuniform dynamic activity distributions, such as those in diseased patients, are to be modeled. However, in the presence of respiratory motion, the lack of spatial constraints in the orthogonal grouping can lead to similar time curves that represent different tracer kinetics. One potential solution to this problem is to use spatial constraints in the orthogonal analysis. Note that differences in tissue count recovery due to partial-volume effect or cardiac motion do not affect the estimation of k1 and k2 by orthogonal grouping, as shown earlier, but only the estimation of fv and rv.
In all patient studies, GFADS yielded LV and RV input functions of the expected form and the 2-compartment kinetic model yielded parametric maps of myocardial tissue extraction and egress as well as spillover fractions from the blood pool into the myocardial tissue. The flows estimated with our approach were well within the range reported previously using an approach based on wavelet transforms, and validated in H215O (9), and were consistent with the clinical status and catheterization results in the limited number of patients considered. CFR was systematically higher in subjects with no prior known CAD (2.1 ± 0.3) than in patients with known prior CAD (1.4 ± 0.1). The estimated CFR values in subjects without obstructive CAD (Table 1) are somewhat lower than those usually found in the literature (∼2.5). This may be related to the fact that these subjects had multiple coronary risk factors (e.g., hypertension, dyslipidemia, diabetes) and were not healthy volunteers. Furthermore, as flow increases, myocardial extraction may no longer be proportional to flow (as in the low-flow regime) and peak flow may be underestimated in dipyridamole studies. This flow dependency of the first-pass extraction fraction can be modeled in our kinetic modeling approach after estimating the LV and RV input functions with GFADS if it is known a priori and, hence, would not significantly affect the accuracy of the flow estimates. It can also be included as a postprocessing step.
In subjects with known CAD who underwent dipyridamole stress, k1 was 30%–50% lower when the VOI drawn on the parametric image corresponded to a region irrigated by a stenotic vessel than when the VOI was drawn over a large myocardial area. Therefore, our approach, based on factor and compartment analyses, allowed discrimination between normal and diseased myocardial blood flow in different subjects, as well as in the same subject. Although GFADS taken alone does not allow the discrimination between different myocardial regions in the same patient (because only one factor is currently used to model the myocardial TAC), the combination of the generalized factor analysis and compartment analysis approaches allows this aim to be achieved by using GFADS to estimate the LV and RV input functions and then by estimating k1 values on a voxel-by-voxel basis.
CONCLUSION
Quantitative dynamic 82Rb PET without arterial sampling is feasible using GFADS and compartmental modeling of 82Rb kinetics; it yields accurate estimates of parameters of absolute myocardial perfusion and kinetics with errors less than 9%.
Acknowledgments
The authors thank Ellen Schlossberg for her assistance with data processing and Dr. Sharmila Dorbala for her help collecting the patient data. This work was supported in part by the National Institutes of Health (NIH) under grants RO1-EB000802 and RO1-EB001989. The contents of this article are the responsibility of the authors and do not represent the views of the NIH.
Footnotes
Received Sep. 2, 2004; revision accepted Apr. 5, 2005.
For correspondence or reprints contact: Georges El Fakhri, PhD, Division of Nuclear Medicine, Department of Radiology, Brigham and Women’s Hospital, 75 Francis St., Boston, MA 02115.
E-mail: elfakhri{at}bwh.harvard.edu
REFERENCES
In this issue
Jump to section
Related Articles
Cited By...
- Myocardial Scar But Not Ischemia Is Associated With Defibrillator Shocks and Sudden Cardiac Death in Stable Patients With Reduced Left Ventricular Ejection Fraction
- Reduced Myocardial Flow Reserve by Positron Emission Tomography Predicts Cardiovascular Events After Cardiac Transplantation
- Integrated Noninvasive Physiological Assessment of Coronary Circulatory Function and Impact on Cardiovascular Mortality in Patients With Stable Coronary Artery Disease
- Integration of Quantitative Positron Emission Tomography Absolute Myocardial Blood Flow Measurements in the Clinical Management of Coronary Artery Disease
- Coronary Microvascular Function and Cardiovascular Risk Factors in Women With Angina Pectoris and No Obstructive Coronary Artery Disease: The iPOWER Study
- Comparison and Prognostic Validation of Multiple Methods of Quantification of Myocardial Blood Flow with 82Rb PET
- How Close Is Close Enough?
- Prognostic Interplay of Coronary Artery Calcification and Underlying Vascular Dysfunction in Patients With Suspected Coronary Artery Disease
- Association Between Coronary Vascular Dysfunction and Cardiac Mortality in Patients With and Without Diabetes Mellitus
- Coronary Vascular Dysfunction and Prognosis in Patients With Chronic Kidney Disease
- Habitual dietary sodium intake is inversely associated with coronary flow reserve in middle-aged male twins
- Improved Cardiac Risk Assessment With Noninvasive Measures of Coronary Flow Reserve
- Quantitative Relationship Between the Extent and Morphology of Coronary Atherosclerotic Plaque and Downstream Myocardial Perfusion
- Inflammation is Related to Coronary Flow Reserve Detected by Positron Emission Tomography in Asymptomatic Male Twins
- Cardiac PET Imaging for the Detection and Monitoring of Coronary Artery Disease and Microvascular Health
- Reproducibility and Accuracy of Quantitative Myocardial Blood Flow Assessment with 82Rb PET: Comparison with 13N-Ammonia PET
- Cardiac Positron Emission Tomography
- Comparison of Positron Emission Tomography Measurement of Adenosine-Stimulated Absolute Myocardial Blood Flow Versus Relative Myocardial Tracer Content for Physiological Assessment of Coronary Artery Stenosis Severity and Location
- Clinical Myocardial Perfusion PET/CT