## Abstract

We have addressed 2 major challenges of ^{82}Rb 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 ^{82}Rb 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 (k_{1}) and egress (k_{2}), as well as LV+RV contributions (f_{v}, r_{v}), by orthogonal voxel grouping. In addition, 13 patients were injected with 2.22 ± 0.19 GBq (60 ± 5 mCi) of ^{82}Rb 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 k_{1}, k_{2}, f_{v}, and r_{v} parametric maps; the average error of estimates of k_{1} was 6.8% ± 3.6%. In all patient studies, our approach yielded robust estimates of k_{1}, k_{2}, f_{v}, and r_{v}, which correlated very well with the status of the subject and the catheterization results. **Conclusion:** Quantitative dynamic ^{82}Rb 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 ^{82}Rb 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 ^{13}N-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 ^{82}Rb positron emitter, the major challenges of parametric ^{82}Rb 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 ^{82}Rb 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 ^{82}Rb 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 ^{82}Rb 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 ^{82}Rb 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 ^{82}Rb 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 ** f_{LS}**:
Eq. 2

Nonnegativity constraints were applied by adding another penalty term, *f*_{n}(*C*,*F*), to the objective function (Eq. 3). ** f_{n}**(

*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 *f*_{uni}(**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 *f*_{uni} under nonnegativity constraints imposed by *f*_{n}(**CR,R ^{−1}F**). The total objective function that was used in the second step (

*f*

_{uni}+

*b*

**f**_{n}) 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 ^{82}Rb 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 k_{1} (mL/min/g) and k_{2} (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 TAC_{i}(t) is the value of the voxel i at time t, TAC_{i} is the TAC of voxel i, and I(t) is the input function (i.e., measured LV function). *k*_{1}^{i}, *k*_{2}^{i}, *f*_{v}^{i}, and *r*_{v}^{i} are the kinetic parameters for voxel *i*. Assuming a soft-tissue density similar to that of water, *k*_{1}^{i} (mL/min/g) characterizes myocardial tissue extraction (inflow). *k*_{2}^{i} (min^{−1}) characterizes myocardial tissue egress (outflow), *f*_{v}^{i} (dimensionless) represents the contribution to the total voxel activity from the blood input function I(t), and *r*_{v}^{i} (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 *k*_{1}^{i}, *k*_{2}^{i}, *f*_{v}^{i}, and *r*_{v}^{i} are zero-valued for some voxels and that, unlike the widely used (1 − f_{v}) formulation, our approach allows fitting the model to all image voxels, including voxels where f_{v} = 1 or r_{v} = 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

*TAC*(_{i}*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^{82}Rb 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*TAC*(_{i}*t*)—that is, we obtained 3D vectors*P*(*TAC*(_{i}*t*))_{τ}, where τ = 1… 3 that described the location of a given dynamic voxel*TAC*(_{i}*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*(*TAC*(_{i}*t*))_{τ}. First, Euclidian distances*d*=_{ij}*P*(*TAC*) −_{i}*P*(*TAC*) 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_{j}*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 ^{82}Rb PET Studies

#### Kinetic Modeling of Dynamic ^{82}Rb Studies.

LV and RV TACs were measured in 10 human rest and dipyridamole stress dynamic cardiac ^{82}Rb 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 ^{82}Rb 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 ^{82}Rb 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 k_{1} values of 0.4, 0.7, and 1.2 mL/min/g and for stress studies with k_{1} values of 0.7, 1.4, and 2.1 mL/min/g. Because ^{82}Rb is an analog of potassium that is taken up rapidly and retained for a long time within myocardial cells (10), we used very low k_{2} values of 0.01 and 0.02 min^{−1} in rest and stress studies. The selected values of k_{1} and k_{2} spanned the range of values reported in the literature in rest and stress studies (10). The lowest two k_{1} values were associated with k_{2} = 0.01 min^{−1}, whereas the highest k_{1} value was associated with k_{2} = 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 (k_{1}, k_{2}) 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).

#### 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 ^{82}Rb 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 mm^{3}). 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 *R*_{single,1} and R_{single,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 ^{82}Rb 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 ^{82}Rb (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 ^{82}Rb from the LV and RV cavities.

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 (k_{1}, unlike k_{2}, 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 (*k*_{1}) and egress (k_{2}), as well as RV and LV contributions (f_{v}^{i}, r_{v}^{i}), were computed over the 30 Monte Carlo simulations. The rate constants k_{1}, k_{2}, as well as f_{v}^{i}, and r_{v}^{i}, 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.

The 2-compartment model kinetic analysis also yielded estimates of k_{1}, k_{2}, f_{v}, and r_{v} parametric maps. The error in estimating of the value of k_{1}, averaged over the myocardium, was 6.8% ± 3.6%. This error was slightly greater for k_{2} (8.4% ± 7.6%) as the simulated k_{2} values were very low (i.e., 0.01–0.02) and, hence, more sensitive to noise. Errors in the estimates of f_{v}^{i} and r_{v}^{i} 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 k_{1}, k_{2}, f_{v}^{i}, and r_{v}^{i} , 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 k_{2} 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 k_{1} 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 k_{1} values at stress was 0.74–1.84 mL/min/g. Furthermore, k_{1} 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 (k_{1}) 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). k_{2} values associated with VOIs were slightly greater (5%–10%) than those with GFADS; however, f_{v} (0.39 ± 0.07) and r_{v} (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 (k_{1 stress}/k_{1 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. k_{1} and k_{2} values were not strongly correlated at rest or stress (*r*^{2} < 0.69).

## DISCUSSION

GFADS was used to estimate noninvasively the LV and RV input functions from the dynamic ^{82}Rb 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 *f*_{v}^{i} and *r*_{v}^{i} but not of k_{1} or k_{2}. 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 ^{82}Rb, because the full width at half maximum of the positron range is 2.6 mm as compared with 1 mm for ^{18}F, so the effect of positron range on spatial resolution and image quality is greater for ^{82}Rb than for ^{18}F. Likewise, modeling random coincidences is of paramount importance in the case of ^{82}Rb cardiac PET because of the high counting rates after injection of 2.22 GBq (60 mCi) of ^{82}Rb, 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 k_{1} and k_{2} by orthogonal grouping, as shown earlier, but only the estimation of f_{v} and r_{v}.

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 H_{2}^{15}O (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, k_{1} 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 k_{1} values on a voxel-by-voxel basis.

## CONCLUSION

Quantitative dynamic ^{82}Rb PET without arterial sampling is feasible using GFADS and compartmental modeling of ^{82}Rb 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