|Home | About | Journals | Submit | Contact Us | Français|
This work reports a principal component analysis (PCA)-based approach for reducing bias in distribution volume ratio (DVR) estimates from Logan plots in positron emission tomography (PET). Comparison has been made of all existing bias-removal methods with the proposed PCA method, for both single-estimate PET studies and intervention studies where pre- and post-intervention estimates are made. Bias in Logan-based DVR estimates is because of the noise in the PET time-activity curves (TACs) that propagates as correlated errors in dependent and independent variables of the Logan equation. Intervention studies show this same bias but also higher variance in DVR estimates. In this work, noise in the TACs was reduced by fitting the curves to a low-dimension PCA-based linear model, leading to reduced bias and variance in DVR. For validating the approach, TACs with realistic noise were simulated for a 11C-labeled tracer with carfentanil (CFN)-like kinetics for both single-measurement and intervention studies. Principal component analysis and existing methods were applied to the simulated data and their performance was compared by statistical analysis. The results indicated that existing methods either removed only part of the bias or reduced bias at the expense of precision. The proposed method removed ~90% of the bias while also improving precision in both single- and dual-measurement simulations. After validation of the proposed method in simulations, PCA, along with the existing methods, was applied to human [11C]CFN data acquired for both single estimation of DVR and dual-estimation intervention studies. Similar results were observed in human scans as were seen in the simulation studies.
For most positron emission tomography (PET) radiotracers, the tracer concentration time course in brain can be modeled as a nonlinear function of rate parameters of a chosen pharmacokinetic model (Carson, 1986). These kinetic parameters are estimated by fitting the dynamic PET data to the kinetic model, using nonlinear estimation techniques. The distribution volume ratio (DVR), a popular index of receptor binding site density related to the non-displaceable binding potential (BPND), is a function of these rate parameters. Nonlinear estimation techniques are computationally expensive, and under noisy conditions may not converge, thus reducing viability for voxel-by-voxel parametric image generation. Logan plot analysis, a computationally efficient alternative, has been proposed for reversible tracers, where ordinary least squares (OLS) can be used to estimate DVR after transformation of dynamic PET data (Logan et al, 1996). Logan plot analysis has been used extensively because of its simplicity and model independence. However, negative bias exists in Logan plot-based DVR estimates, owing to the noise present in time activity curves (TACs) that propagates as correlated errors in the dependent and independent variables of the Logan plot equation (Slifstein and Laruelle, 2000). We have found that this noise propagation is more pronounced (resulting in higher bias and variance) in PET studies that involve a pharmacological or behavioral intervention, primarily because of the restricted temporal range of data available for DVR estimation. In these studies, two pharmacokinetic measurements are obtained from a single-tracer injection, before and after the intervention. It is important to have a method to reduce bias as well as maintain high precision in the parameter estimates from such experiments.
To reduce the bias in single-measurement PET studies, various methods have been proposed previously for the arterial sampling approach (Logan et al, 2001; Varga and Szabo, 2002; Ichise et al, 2002; Ogden, 2001). The operational Logan plot equations for both arterial sampling and reference region approaches have the same problem of correlated errors in both dependent and independent variables. The references mentioned above address this problem of correlated errors and hence should be applicable to both arterial sampling and reference region approaches. We found that these existing methods either removed the bias at the expense of precision or removed only part of the bias. In this paper, we introduce a new approach for reducing bias using principal component analysis (PCA), which is a feature extraction technique used to simplify a data set by reducing its dimensionality while maintaining its relevant characteristics. Principal component analysis has been used extensively in nuclear imaging, including PET, in the spatial domain (Barber, 1980; Pedersen et al, 1994; Thireou et al, 2003; Razifar et al, 2006), where the investigators have analyzed the images of the coefficients of individual principal components. Principal component analysis has been previously used in temporal domain in PET for clustering the TACs before parameter estimation (Kimura et al, 2002) and separation of TACs in dual-tracer scans (Kadrmas and Rust, 2005). Frequency analysis of dynamic structures, where the principal components are rotated to avoid negative values, has been used in temporal domain in single photon emission computed tomography (Sitek et al, 1999). In this work, we apply PCA to achieve temporal smoothing of PET TACs to reduce the bias in DVR estimates.
The differential equations for a two tissue compartment model for a PET tracer are shown below:
where Cp(t) is the arterial plasma input, CND(t) is the radioligand concentration in the nondisplaceable compartment, CS(t) is the radioligand concentration in the specific compartment, and K1, k2, k3, and k4 are the kinetic parameters of the model. The total concentration of the radioligand in the tissue is given by
where IR(t) is the tissue impulse response and is the sum of two exponential terms as shown below:
with α1, α2, A1, and A2 being combinations of the individual rate parameters, K1 to k4.
The non-displaceable binding potential (BPND), a commonly used index of receptor binding density, is the ratio of the equilibrium concentration of specifically bound to nondisplaceable radioligand in any given region (k3/k4). Direct estimation of BPND from equations (3) and (4) suffers from the practical difficulties of arterial blood sampling and the possible errors associated with it, as well as being computationally expensive because of the nature of nonlinear estimation techniques. Conversely, the reference region approach proposed by Logan et al (1996) avoids arterial sampling by using a region assumed to be devoid of specific binding sites to estimate a similar index of binding termed DVR (= BPND + 1).
In the reference region-based Logan plot approach, DVR is estimated using OLS by transformation of the dynamic PET data as shown below:
where INT is an intercept term, Ti is the midpoint time of the ith frame, and is the chosen k2 parameter for the reference region. Equation (5) becomes linear over multiple time measurements of the target and reference tissues after time T′, and DVR is estimated by calculating the slope of the line. The numerators in both the dependent variable
and independent variable
have integral terms and are relatively noise free as compared with the denominators. It has been reported that the presence of the noisy term (CT(Ti)) in the denominators on both sides of equation (5) is the cause of correlated errors leading to negative bias in DVR estimates (Slifstein and Laruelle, 2000). This bias is more pronounced in regions with higher DVR values. Reduction of noise in the measured TACs will assist in reducing the correlated errors leading to a reduction in the bias. In this work, we propose a PCA-based approach to achieve this.
A pharmacological intervention study involves obtaining two measures of one neuropharmacological parameter from a single PET study. A single tracer is administered, and at some point during the PET data acquisition, an interventional challenge or perturbation of the system is made with the intent of changing in vivo distribution of the radiotracer. Logan plots are used to estimate DVR before and after intervention (Koeppe, 2002). The bias and variance is more pronounced in this case because fewer data points are available for DVR estimation compared to when a single measurement is made. Application of the existing bias-removal techniques to such intervention studies has not been reported in the literature. Except for the temporal range of data used, the procedure for estimating DVR in intervention and single-measurement studies is the same. The Logan plots for dual-measurement intervention studies have the same noise-induced bias associated with them as in single-measurement studies. Thus, although the existing methods have been proposed for single-measurement studies, they should be applicable to the dual-Logan case without any modifications. We applied the proposed PCA-based approach as well as the existing methods to reduce the bias in PET intervention studies (both simulations and human studies).
Principal component analysis is a feature extraction technique used to simplify a data set by reducing its dimensionality while maintaining its important characteristics (Faraway, 2004). In feature extraction, data space is transformed into a ‘feature’ space having the same dimensions as the data set. This transformation is such that the data set can then be represented by a reduced number of dominant ‘features’ while retaining all the important intrinsic characteristics of the original data. Using PCA, each original data vector of dimension ‘p’ can be expressed as a linear combination of ‘p’ orthogonal basis vectors. By limiting the number of basis vectors to ‘q’ dominant vectors (q<p), data can be represented with a reduced dimensionality, thus reducing the noise and preserving all the important features in the data. Using PCA, the noise in the PET TAC values (CT(Ti)) can be reduced, thus reducing the correlated errors in equation (5), which is the cause of bias in DVR estimates.
The classical model for a PET TAC is a nonlinear function of rate parameters (equations (3) and (4)). Here, we represented dynamic PET data using a PCA-based linear model. Thus, a noisy tissue curve vector for the jth voxel, j = [CTj(T1), CTj(T2),, CTj (Tp)]′ Rp×1, from a PET study involving ‘p’ temporal frames can be expressed as
where GRp×q is the system matrix to be constructed using PCA (q<p), jRq×1 is the coefficient vector, and j Rp × 1 is the vector of residuals. For obtaining the system matrix, a training set is required. We used the TACs from all the voxels in the PET image volume as the training set.
Let the training set contain ‘d’ TACs from the dynamic PET data. These ‘d’ curves can be arranged row-wise in the matrix form as follows:
Matrix X0 is obtained by subtracting the column mean vector () from each row of X. The principal components are the eigenvectors of the sample covariance matrix of and are obtained using singular value decomposition. The total number of principal components is equal to p, the number of frames in the PET study. These components can be ordered in the decreasing order of their eigenvalues as F=[1,2, …, p], where iRp×1 and FRp×p. The mean vector and a subset of the principal components (depending on their eigenvalues) are chosen to construct the system matrix G = [12 q] (q<p). The coefficient vector j (equation (6)) for each noisy curve j can be estimated by least squares minimization:
The PCA-based fitted tissue curve can then be obtained as shown below:
The process of obtaining the principal components using PCA from the training set has a closed-form expression and hence is straightforward. Selection of ‘q’ for a particular tracer is made from simulations. Training data and test data are simulated using the literature range of tracer parameter values and the expected shape of arterial input function. The principal components are obtained from the simulated training set and ‘q’ is selected to be the minimum number such that G = [12 … q] is a valid system matrix for simulated test data.
In human studies, the training data and the test data are the same and include all the TACs in the PET image volume. Principal components for each subject are obtained from the training set, and system matrix is formed using the number ‘q’ selected from simulations. This system matrix is then used to fit the test data.
Simulations were performed mimicking the relatively slowly equilibrating radiotracer [11C]carfentanil (CFN), a reversible μ-opioid agonist. A partial bolus–continuous infusion administration protocol was simulated for tracer administration (60% bolus and 40% infusion). This was the same bolus to infusion ratio used in the human scans performed at our institution. The local rate of delivery of the radiotracer and the clearance rate of the radiotracer to the plasma were same for the target region and the reference region , although this assumption is not required by the reference region method. The training set was simulated using the kinetic parameter ranges reported previously (Endres et al, 2003) as shown in Table 1. The training set consisted of 2,496 noisy curves (all with unique kinetic parameter combinations) of a 70 min, 16-frame scan (4 × 0.5 min, 3 × 1 min, 2 × 2.5 min, 2 × 5 min, and 5 × 10 min), which is the same protocol used for human scans at our institution. The principal components were obtained from this training set.
The test data consisted of six hypothetical regions (Table 1 and Figure 1). Four of them with receptor binding, having true DVR values of 5, 4, 3, and 2, and one without receptor binding, having true DVR = 1, were sampled from the training set. It is possible for a small region of the brain to have very different kinetics from those elsewhere. As the proposed PCA method is training set based, its sensitivity to a region of unusual kinetics not part of the training set should be evaluated. To examine this, a sixth region with substantially greater DVR than any other curve present in the training set was added to the test set (DVR = 8). Noiseless TACs for these six regions are shown in Figure 1. Realistic voxel-level noise was added to obtain 1,024 realizations for each of these curves. The noise model used for the simulations related to the counting statistics similar to equation (8) in Logan et al (2001). A scale factor that generated noise level similar to that in voxel-wise PET TACs was chosen. The principal components obtained from the simulated training data were used to construct the system matrix and fit the curves in the simulated test data. DVR estimates were then obtained by replacing the noisy TAC values in the denominators of equation (5) with the values of the PCA-based fitted curves. T′ was chosen to be 20 min based on the known characteristics of CFN. The integral in the numerator of equation (5) was calculated using trapezoidal approximation (Ogden, 2001). The curve simulated for the region with no binding (DVR = 1) was used as the reference region in the Logan analysis.
It must be noted that as no blood samples are measured in reference region approaches, correction for blood volume component is not possible. Not correcting for blood volume component will introduce a systematic bias in the DVR estimates unrelated to the noise-induced bias (Logan et al, 1996). To check the magnitude of this systematic bias in the proposed and existing methods, separate simulations were also performed with blood volume component using the following model:
where C(t) is the total radioligand concentration measured in a voxel, VB is the blood volume component (chosen to be 3.5%), and CB(t) is the blood radioligand concentration.
Intervention studies were simulated with the same input function and rate parameters as for single-measurement studies (Table 1) with a naloxone intervention simulated 40 min after initiation of the scan. The simulated data consisted of 100 min, 19-frame scan (4 × 0.5 min, 3 × 1 min, 2 × 2.5 min, 2 × 5 min, and 8 × 10 min), which is the same protocol used for human intervention scans at our institution. The training set consisted of TACs with a linear decrease in k3 to 40%, 50%, and 60% of the initial k3 value over 10 min after the intervention to model the loss in binding potential after the naloxone intervention. The test set was simulated with a 50% decrease in k3. The DVR values were estimated from Logan plots using ‘early’ data from 20 to 40 min for pre-intervention DVR (DVRpre) and ‘late’ data from 60 min until the end of scan data for post-intervention DVR (DVRpost) as shown in Figure 2 (noiseless case). DVRpre is estimated from 20 to 40 min, when the Logan plot has become linear. The perturbation at 40 min (denoted by the vertical gray line) increases receptor occupancy causing the Logan plot to bend (40 to 60 min) before becoming linear again at around 60 min. The 60 to 100 min period can be used to estimate DVRpost. For both DVRpre and DVRpost estimations, the integral terms in the numerators of equation (5) are calculated from the beginning of the scan.
DVR values for the above simulations were also estimated using existing bias-removal techniques. Thus, DVR values were estimated by the following six methods:
The TLS method (Varga and Szabo, 2002) takes into account errors in both dependent and independent variables of the Logan plot to estimate DVR. This can be classified as a method using an alternate cost function. The LEGA (Ogden, 2001) and MA1 (Ichise et al, 2002) approaches can be classified as methods rearranging the original Logan plot equation. In LEGA, a rearrangement of the operational Logan plot equation makes the error term additive. The problem is then solved by maximum likelihood approach. In MA1, equation (5) is rearranged to bring the noisy term on one side of the equation and DVR is estimated by taking the ratio of two regression coefficients. The GLLS method (Logan et al, 2001) can be classified as a temporal smoothing method. In this method, each noisy TAC is separated into two segments, which are individually fitted to one compartment model using GLLS (Feng et al, 1996). The smoothed segments are put together to get a smooth TAC, which is then used to obtain Logan DVR estimates. Implementation of GLLS requires choosing the frame where the TAC is to be segmented. For single-measurement simulations, the 11th frame (15 to 20 min) was found to be the optimal segmentation point for bias reduction and was used as segmentation point for human data as well. The proposed PCA approach can also be classified as a temporal smoothing method where a PCA-based lower dimension linear model is used to reduce the noise in the tissue curves.
For intervention studies, the intervention is applied at the 14th frame (40 to 50 min) and DVR is estimated from the segments before and after the 14th frame using Logan analysis. As mentioned before, the existing methods, although proposed for single-measurement studies, are also applicable to ‘dual-Logan’ intervention studies. Although this is true for GLLS as well, some practical difficulties exist. To apply the GLLS method exactly as proposed by Logan et al (2001) in the dual-Logan case, temporal segments of the TAC before and after intervention may have to be further cleaved into two parts each. This, although possible, is not practical because of a limited number of points in the TACs (11 pre-intervention and 6 post-intervention). Thus, we implemented GLLS without dividing the segments before and after intervention any further. It should be noted that to fairly compare this method with the other methods, the segments before and after intervention may need to be separated further (work not done).
To apply the proposed approach in single-measurement human PET studies, 12 subjects were imaged with [11C]CFN for 70 min and binned into 16 frames (4 × 0.5 min, 3 × 1 min, 2 × 2.5 min, 2 × 5 min, and 5 × 10 min). All scans were performed on a Siemens ECAT Exact HR+ scanner. Images were reconstructed after Fourier rebinning using 2D ordered subsets expectation maximization, 4 iterations, 16 subsets, and no postprocess smoothing, resulting in images with isotropic resolution of approximately 6 mm FWHM. Injected doses of [11C]CFN were 555 to 666 MBq (15 to 18 mCi) at a specific activity of > 2,000Ci/mmol. Carfentanil was administered as partial bolus followed by continuous infusion. Subject motion across frames was corrected using Neurostat (University of Michigan; Minoshima et al, 1994, 1993). Volumes of interest, including the occipital cortex reference region, were obtained using a standardized volume-of-interest template after reorientation and nonlinear warping to the stereotactic Talairach atlas (Talairach and Tournoux, 1988) using Neurostat routines. DVR parametric images were obtained using the proposed and existing methods. The PCA training set consisted of TACs from all the individual voxels in the brain volume and was used to obtain the principal components for each subject. The system matrix was then constructed using the first ‘q’ principal components (‘q’ selected from simulations). All individual voxel TACs in the image volume were then fitted to the system matrix to obtain smooth TACs with reduced noise. DVR values were estimated by substituting these smoothed TAC values in place of the noisy TAC values in the denominators of equation (5). DVR estimates were also obtained by using the existing bias-removal methods and compared with the proposed PCA approach.
To apply the proposed approach in intervention studies, four human studies were performed using [11C]CFN with a naloxone dose administered starting 40 min after CFN administration, to block tracer binding. Naloxone was administered as a 0.004 mg/kg bolus followed by 0.00004 mg per kg per min for the remainder of the study, a dosage designed to block approximately 50% of the specific binding sites. Four control subjects were also imaged where no intervention was given during the scan and thus no change in binding was expected between ‘early’ and ‘late’ DVR estimates. Data for these eight subjects were collected as in the single-measurement studies except that scans were acquired for 100 min and binned into 19 frames (4 × 0.5 min, 3 × 1 min, 2 × 2.5 min, 2 × 5 min, and 8 × 10 min). Except for the temporal range of data used, the procedure for estimating DVR using the proposed and existing methods in intervention studies is same as that for single-measurement studies as explained above.
The mean vector () and first three principal components (1, 2, and 3) obtained from the simulated CFN training data without blood volume component are shown in Figure 3. The first three components had eigenvalues of 1.24 × 103, 278.75, and 7.91, respectively, whereas the second last and last components had eigenvalues of 6.36 × 10−32 and < 10−50, respectively. This progression of eigenvalues (eigenvaluei–1 eigenvaluei) indicates that the first few components are the dominant vectors and capture most of the variance in the training set. The box plots in Figure 4 show the comparison of DVR estimates obtained using OLS with those estimated by increasing numbers of components for true DVR = 5 (PCAn G = [12 n]). In the proposed PCA approach, we fit the noisy TACs as a linear combination of the mean vector in addition to the chosen principal components (16 possible components). In all, 17 PCA-based models are possible; the model with the smallest degrees of freedom is PCA0 with the mean vector alone G = ), and the one with the largest degrees of freedom is PCA16 with mean vector and all the 16 principal components (G = [1 2 16]). Too few components were insufficient to model the data and yielded biased estimates with high precision, as seen for PCA0, which used only the scaling of the mean vector. As the number of components was increased, the bias was reduced at the expense of precision (PCA1, PCA2, and PCA3). However, as the number of components exceeded three, the bias returned, as the added components began to fit noise in the data. Using all the 16 principal components (PCA16), as expected, gave back the original noisy data, and the box plot in this case is identical to the box plot for OLS. As we include the mean vector in all the PCA models, PCA15 (the mean vector and 15 components) can perfectly describe the 16-frame data. Thus, PCA15 and PCA16 give identical results.
The PCA1 (G = [ 1]) provided the best bias-variance trade-off for true DVR = 5. This implied that a linear combination of only two curve shapes could sufficiently approximate the simulated test data with DVR = 5 (q = 1). To check the validity of this finding for all the simulated TACs in the test set, we analyzed the residuals between the true noise-less TACs and the PCA1-fitted TACs ( , 1≤ j ≤1,024). For the jth realization, the residual at the ith frame is given by (1≤ i ≤ p, where p = 16 is the number of frames in the study). These residuals were normalized by the true tissue curve values to obtain the percent-normalized residuals
For a model to be valid, the mean of these percent-normalized residuals must lie close to zero (Carson, 1986). Figure 5 shows that for TACs with receptor binding that were part of the training set (DVR = 5, 4, 3, and 2), the mean of percent-normalized residuals was close to zero (between + 3% and −3%). For reference region TAC (DVR = 1, k3 = 0) and for the atypical TAC not included in the training set (DVR = 8), the mean of the percent-normalized residuals was significantly different from zero, indicating that PCA1 was not a valid model for these curves. This makes intuitive sense as the reference region curve (shown in Figure 1, DVR = 1) has a more complicated shape compared with the curves with receptor binding because of rapid clearing of the tracer. The PCA3 (mean vector and three components) approach was found to be a valid model for the reference region. It should be pointed out that the actual reference region TAC used in the Logan plot calculations is not fitted using PCA but is the raw TAC. As the reference region is not derived from a single voxel, but is obtained from a volume of interest large enough to have minimal statistical noise, the raw reference region curve can be used for all DVR estimations without introduction of bias or loss of precision. Furthermore, the DVR estimates in regions of very low binding were found to be nearly unbiased even with nonzero residuals resulting for the PCA1 fit as seen in the last column of Table 2.
Invalidity of the PCA1 model for the atypical curve is also expected, as the success of the PCA-based approach depends on the diversity of the training set, and if a curve shape is significantly different from any in the training set, the PCA1 model may not be adequate. However, in human studies, the training set consists of all the TACs in the PET image volume. Thus, a region with unusual kinetics, if any, will contribute curves to the training set. If we include curves similar to the atypical curve in the simulated training set (1% of the total curves), the fit improves as seen in Figure 5 (DVR = 8+).
The standard deviations of the percent-normalized residuals for PCA1-based TACs calculated above were approximately 50% of that for original noisy tissue curves. Thus, from this result and Figure 5, we could conclude that PCA1 model was valid for the curves with receptor binding that were part of the training set and provided TACs with reduced noise.
More detailed results have been presented for the curve simulated with DVR = 5. The histograms in Figure 6 show the distribution of the DVR values estimated by each of the existing methods versus PCA1. A Gaussian kernel with unbounded support was used to smooth the histograms. These plots were useful for visualizing asymmetries and tails in the distribution of the estimated DVR values. The vertical line denotes the ‘true’ DVR value. The modes of the histograms of all the existing methods either lie below the true DVR value, show heavy tails, or both. The PCA approach has its mode at the true value but has a heavier tail than PCA1. The PCA approach was seen to be superior to all the existing methods. Figure 7 depicts the trade-off between precision and bias in the DVR estimates of the different approaches, again for DVR = 5. An ideal estimator has no bias and minimal variance, and would lie near the origin. The OLS method (superimposed with PCA16) is on the far right, denoting the bias present because of noise. The TLS method removes only a portion of the bias. The GLLS, LEGA, MA1, PCA2, and PCA3 methods reduce the bias but also have poorer precision. The PCA1 method shows the best performance and is closest to the origin than any other method. Figure 8 shows box plots for DVR estimates from all methods (DVR = 5). The median and range of DVR estimates using OLS are below the true value, as expected. The bias in the median is reduced only slightly using TLS with a decrease in precision. The GLLS, LEGA, MA1, and PCA2 methods reduce the bias almost entirely but have high variance. Principal component analysis 1 eliminates bias almost completely and yields the best precision estimates compared with any of the methods.
Table 2 summarizes the means and standard deviations of the estimated DVR values in units of percentage of true value for all the simulated test curves. Although the trends for DVR<5 (last four columns) are similar to those seen in Figure 8 (DVR = 5), bias for OLS is lower for curves with lower DVR values, as expected (Slifstein and Laruelle, 2000). The first two columns show the statistics for the atypical TAC with bias in the PCA1 estimate of DVR when the atypical curve is not part of the training set (DVR = 8, column 1), and reduction in bias when atypical curve shapes are included in the training set (DVR = 8a, column 2). Inclusion of curve shapes similar to the atypical curve in the training set (1% of total curves) improves the fit but does not remove the bias completely for PCA1. The representation of atypical curve shapes in the training set was required to be at least 5% for the bias to be removed completely for PCA1. Thus, very small regions with unusually shaped TACs, if identified before Logan analysis, could be fitted with PCA2, where bias is removed almost completely but at the expense of precision.
Figure 9 shows the box plot for DVR estimates from all methods for curves simulated with but not corrected for the blood volume component (DVR = 5). The box plots have similar bias–variance trends as in Figure 8, apart from the consistent bias seen in Figure 9 for all methods because of uncorrected blood volume component. This bias is minimal (<5%) compared with the noise-induced bias in OLS (>10%) for DVR = 5 and is even smaller for lower DVR values.
Analyses similar to those shown in Figures 4–9 were performed for intervention studies to determine the appropriate number of principal components for fitting the noisy data. Figure 10 shows percent absolute bias versus percent standard deviation in the DVR estimates for both pre-intervention (‘true’ DVRpre = 5) and post-intervention (‘true’ DVRpost= 3) data. As the bias in the estimated DVR values is higher for larger DVR values, in general the bias is higher in DVRpre than DVRpost. It can be seen that PCA1(pre) and PCA3(post) provide the best bias– variance trade-off (by virtue of their proximity to the origin) compared with all the other methods. In fact, DVR estimates using all PCA-based approaches tested (PCA1, PCA2, and PCA3) cluster close to the origin and reduce bias while maintaining good precision. All the existing bias-removal methods reduce bias at the expense of precision.
For the 12 human scans in the single-measurement case, eight regions of interest were placed on the estimated DVR images of each subject obtained from all the methods. The means of DVR values and means of the standard deviations within these regions from 12 subjects are reported in Table 3. The top 50 percentile pixels in these regions of interest were used to obtain the statistics. The performance of the methods in the human data was seen to match closely with the performance in simulation studies. The TLS estimates are higher but quite close to that of OLS for all regions. The estimates of GLLS, LEGA, PCA1, PCA2, and PCA3 are higher than OLS as seen in single-measurement simulation studies indicating reduction in bias. Mean values for MA1 are higher than the other methods because of the presence of outliers. The means of the standard deviations (across-subject means of the within-region standard deviations) also follow the same trend as seen in simulations with PCA1 having the lowest and MA1 having the highest variance compared with the other methods. It is important to note that in general the PCA method yields the smallest standard deviations, which is the strength of the approach.
Figure 11 shows DVR images obtained from the data obtained before and after 40 min of scanning in one of the control subjects where no intervention was given; hence, DVRpre and DVRpost images should be identical. The OLS, TLS, and GLLS methods show a decrease in global DVR, although no receptor blocking in present. The LEGA and MA1 methods give noisy DVR images, especially for the ‘late’ measure, and hence should not be used to create parametric images in intervention studies. The PCA1 images show little to no change between ‘pre’ and ‘post’ images and do not show any bad fits. The PCA2 and PCA3 images show higher variance but little to no change between early and late DVR images. To quantify these effects, scatter plots of pre and post binding potential estimates for voxels with receptor binding (BPND > 1) were made for all the subjects, and the mean values of correlation coefficient, slope, and intercept are reported in Table 4 (n = 4 for both groups). As errors exist in both abscissa and ordinate, the regression slopes and intercepts were estimated using TLS optimization.
In the case of the controls, the ideal estimator will have correlation coefficient = 1 and slope = 1, with intercept at the origin indicating no variance and no change in occupancy. The PCA1 method behaves closest to the ideal estimator for the control case. The other methods show lower correlations and slopes less than 1, indicating higher variance and an apparent change in occupancy because of higher bias in ‘late’ than in ‘early’ DVR estimates. Figure 12 shows DVR images before and after intervention in one of the subjects who received naloxone. Noise properties of DVR images are similar to those in the no intervention case. The OLS, TLS, and GLLS methods show the expected decrease but are noisier. Similar to the results in the control study, LEGA and MA1 images are much noisier and would have limited utility for intervention studies. The PCA1 images have the lowest variance and show the expected decrease in DVR values; PCA2 and PCA3 again have higher variance, but show the expected effects of naloxone on DVR estimates. The results for scatter plots in the naloxone intervention cases are summarized in Table 4 (right columns). The PCA-based methods have higher correlations than the other methods and the mean slope for PCA1 = 0.54 corresponds well with the expected ~50% occupancy of receptor sites at the chosen naloxone dose.
The performance of each method varied little from subject to subject, as seen by the small across-subject standard deviation values for the mean statistics as reported in Table 4. The performance of all the methods depends primarily on the noise level in the TACs (higher bias and variance with higher noise). As all the subjects were imaged with the same PET scanner, the noise level in all scans was similar and so was the performance of each method from subject to subject.
The negative bias in the DVR estimates from Logan plots is because of the noise in the PET TACs that propagates as correlated errors in dependent and independent variables of the Logan plot equation. Although several methods for reducing this bias have been proposed, these methods either remove only a portion of the bias, or reduce bias at the expense of precision. Parameter estimation by TLS removes only a portion of the bias. The GLLS, LEGA, and MA1 methods are quite effective in removing bias, but show poorer precision, especially for intervention studies where two DVR estimations are required. In the proposed approach, fitting the dynamic PET data to a low-dimension PCA-based linear model maintains the appropriate kinetic shape of the TACs while removing noise and hence the source of the correlated errors in the Logan plot equation. This provides DVR estimates with minimal bias and reduced variance. In single-measurement simulation studies, the linear combination of just two curve shapes (the mean vector and first principal component) was successful in modeling the simulated test data curves sampled from the training set and reduced bias and variance associated with their DVR estimates. If a test data curve was not part of the training set, PCA1 was insufficient to model it, leading to bias in DVR estimates as seen for the atypical curve (DVR = 8). In human studies, however, TACs from all the voxels in the image volume are used as training data to obtain the PCA-based linear model. A region of unusual kinetics, if present, will contribute to the total TACs in the training data. Thus, a region of unusual kinetics contributing 1% of the total TACs to the training data was simulated. By including curve shapes similar to the one with DVR = 8 in the training set (1% of total curves), the bias in DVR estimates was reduced but not completely removed. Small regions with atypical TACs, if identified before DVR estimation, could be fitted with PCA2, where bias will be removed at the expense of precision.
Adding another component to the PCA1 model (PCA2) also removed bias, but increased the variance of the DVR estimates. Although the increase in variance between PCA1 and PCA2 is expected because of an increase in degrees of freedom, the magnitude of this increase will vary from tracer to tracer and will depend on the particular curve shapes under consideration, and thus needs to be evaluated for each radiotracer application.
Another source of bias, unrelated to the noise-induced bias, exists for reference region approaches because of lack of accounting for the blood-borne component of the PET signal. This was observed to cause a small but consistent bias in all the methods, but was minor compared to the noise-induced bias in OLS.
It should be pointed out that in both the simulation studies and the human scans, we used a partial bolus followed by continuous infusion for administration of the radiotracer. This may reduce the kinetic complexity of the TACs compared to bolus injection studies, and hence limit the number of principal components required to describe the shape of the TACs. When applying the PCA approach to other PET ligands or radiotracer administration protocols, the optimal number of components for describing the TAC shapes needs to be reevaluated. However, in our experience, the mean vector plus at most two principal components seem sufficient for modeling standard single-measurement PET studies for most radiotracers, whereas the mean vector plus at most three components are sufficient for interventional PET studies that may have more complex kinetic behavior.
Obtaining the components from training data has a closed-form solution and hence is not computationally expensive. While the simulation studies are a necessary step in the validation of any newly proposed approach, the final evaluation must involve real data. Results presented in Tables 3 and and44 and Figures 11 and and1212 clearly show the utility of the PCA approach in human data. In particular, it should be emphasized that for intervention studies, the PCA approach is the only approach that produced voxel-by-voxel parametric images with acceptable noise levels and with little observable bias. Figure 11 shows that for the no-intervention case, only PCA1 yielded ‘early’ and ‘late’ DVR values that were of the same magnitude and of the same statistical quality. A similar advantage of PCA1 compared with the other methods for the intervention case is seen in Figure 12, with the expected decrease in DVR because of naloxone. Scatter plots of pre- versus post-intervention binding potential estimates show that PCA-based methods yield the highest correlation and have closer to the expected slopes compared with the other methods (Table 4). These results show both the sensitivity and specificity of the PCA-based approach. The PCA1 approach produced unchanged ‘early’ and ‘late’ measures when no intervention was given, and also showed the ability to detect intervention-induced changes when present. The PCA1 method, however, may not have the necessary degrees of freedom to model very large displacements. The PCA2 or PCA3 methods will capture these large displacements, but at the expense of higher variance.
The proposed PCA method is simple to implement and computationally fast. For relatively slowly equilibrating tracers like [11C]CFN, PCA1 produced DVR images that had less bias, lower variance, or both, compared with existing bias-removal methods. Advantages of the PCA approach were seen not only for standard single-measurement PET studies, but were especially impressive when applied to dual-measurement intervention studies where pre- and post-challenge DVR estimates are required.
This work was supported by the Department of Energy Grant DE-FG02-87ER60561 and National Institutes of Health (NINDS) Grant NS-15655.