A Strategy for Removing the Bias in the Graphical Analysis Method

Abstract

The graphical analysis method, which transforms multiple time measurements of plasma and tissue uptake data into a linear plot, is a useful tool for rapidly obtaining information about the binding of radioligands used in PET studies. The strength of the method is that it does not require a particular model structure. However, a bias is introduced in the case of noisy data resulting in the underestimation of the distribution volume (DV), the slope obtained from the graphical method. To remove the bias, a modification of the method developed by Feng et al. (1993), the generalized linear least squares (GLLS) method, which provides unbiased estimates for compartment models was used. The one compartment GLLS method has a relatively simple form, which was used to estimate the DV directly and as a smoothing technique for more general classes of model structures. In the latter case, the GLLS method was applied to the data in two parts, that is, one set of parameters was determined for times 0 to T₁ and a second set from T₁ to the end time. The curve generated from these two sets of parameters then was used as input to the graphical method. This has been tested using simulations of data similar to that of the PET ligand [¹¹C]-d-threo-methylphenidate (MP, DV = 35 mL/mL) and ¹¹C raclopride (RAC, DV = 1.92 mL/mL) and compared with two examples from image data with the same tracers. The noise model was based on counting statistics through the half-life of the isotope and the scanning time. Five hundred data sets at each noise level were analyzed. Results (DV) for the graphical analysis (DV_g), the nonlinear least squares (NLS) method (DV_nls), the one-tissue compartment GLLS method (DV_f), and the two part GLLS followed by graphical analysis (DV_fg) were compared. DVFG was found to increase somewhat with increasing noise and in some data sets at high noise levels no estimate could be obtained. However, at intermediate levels it provided a good estimation of the true DV. This method was extended to use a reference tissue in place of the input function to generate the distribution volume ratio (DVR) to the reference region. A linearized form of the simplified reference tissue method of Lammertsma and Hume (1996) was used. The DVR generated directly from the model (DVR_fl) was compared with DVR_fg (determined from a “smoothed” uptake curve as for DV_fg) using the graphical method.

Keywords

Positron emission tomography Kinetic modeling Graphical analysis Distribution volume

The interpretation of positron emission tomography (PET) data in physiologic terms is dependent upon the estimation of model parameters. These model parameters which incorporate the effects of both the tissue uptake and the plasma delivery of the tracer are used to quantify PET data and thus allow intersubject comparisons. The nonlinear least squares (NLS) methods can provide statistically accurate estimates of model parameters. These methods, which are based on a particular model structure, generally require considerable computation time. The linearized version of the standard compartment models (Blomqvist, 1984; Evans, 1987) provide a more efficient method of parameter estimation. Graphical analysis represents a further simplification in which the set of linear equations is transformed into a linear plot (Logan et al., 1990). Although this is applicable to a multicompartment system, only two parameters are determined, the slope and the intercept, which are combinations of the model parameters. In the case of reversibly binding ligands, the slope is the distribution volume (DV). The strength of the graphical analysis is that it does not require a particular model structure since in many cases one model may not equally fit all data sets from the same structure and region of interest. When the model structure does not quite fit the data, a bias can be introduced into the model parameter of comparison, generally the distribution volume (DV). In the case of noisy data, the linearized equations can also introduce a bias, because the error term at any given time point also contains the error terms at the earlier time points (Feng et al., 1993, 1996). As a result, the graphical method on average will underestimate the DV and this underestimate is greater with larger DVs (Hsu et al., 1997; Slifstein and Laruelle, 1999; Abi-Dargham et al., 2000). However, the effect on any given data set depends upon the nature of the noise; for example, when one particular data set deviates in a seemingly nonrandom way from the “true” data (Figs. 2A and 3A), all analysis methods will exhibit an apparent bias. This can only be addressed by reducing the noise source in the data or averaging over data sets, etc.

Feng et al. (1996) have developed a solution to the bias problem in the case of specific compartmental models called the generalized linear least squares (GLLS) method. In the case of the one-tissue compartment model, there are two parameters to be estimated and the GLLS equations for parameter estimation take on a relatively simple form involving the inversion of a 2 × 2 matrix. This one-compartment GLLS method was adapted to use as a smoothing technique for more general classes of model structures. This is accomplished by applying the GLLS method to the data in two parts, that is, to determine one set of parameters for times 0 to T₁ and a second set from T₁ to the end time. The curve generated from these two sets of parameters can then be used as input to the graphical method to generate an unbiased estimate of the DV. This has been tested using simulations of data similar to that of the PET ligands [¹¹C]raclopride which binds to the D2 receptor (Volkow et al., 1993) and [¹¹C]-d-threo-methylphenidate which binds to the dopamine transporter (Volkow et al., 1995). This combination of the GLLS method and the graphical method provides the possibility of retaining the model independent type of analysis without the bias inherent in the linear methods while still maintaining a fairly simple method of analysis. Results from this method are compared with the DV calculated directly from the solution to the one-compartment GLLS method and to the NLS method as well as to graphical analysis without smoothing.

This method was also applied to the simplified reference tissue model of Lammertsma and Hume (1996) and Gunn et al. (1997). The equation was modified to allow a linear solution for k₂ as in Feng's method. Estimates of three parameters were generated in this case as opposed to two when the input function is measured. The same two part procedure was used to smooth the data as was done with the DV and the graphical method was applied to the smoothed data using the reference region and an average efflux constant (Logan et al., 1996). This was compared to the graphical method without smoothing and to the distribution volume ratio (DVR) determined directly from the linear solution of the simplified reference tissue model (Lammertsma and Hume, 1996).

The input function and reference tissue techniques were also applied to human PET data from a ¹¹C raclopride (RAC) and [¹¹C]-d-threo-methylphenidate (MP) study using regions from the basal ganglia where both tracers concentrate.

Theory

The one-tissue compartment model is described by

\begin{matrix} C p ⇄_{k_{2}}^{K_{1}} C_{1} \\ \frac{d C_{1}}{d t} = K_{1} C p (t) - k_{2} C_{1} \end{matrix}

(1a)

where Cp(t) is the plasma tracer concentration at time, t; C₁ is the tissue tracer concentration; K₁ and k₂ are the transport constants, plasma to tissue and tissue to plasma, respectively. The linear form of this model is the set of equations (for scan times t₁)

C_{1} (t_{i}) = K_{1} \int_{0}^{t_{i}} C p d t - k_{2} \int_{0}^{t_{i}} C_{1} (t) d t + ξ_{i}

(1b)

where the equation errors, ζ_i, are not statistically independent because each succeeding one depends upon the previous ones. This can result in biased parameter estimates (Feng et al., 1996). To overcome the bias problem in the solution of the one-compartment model, Feng et al. (1993) introduced the GLLS method that removes the bias. The GLLS form of Eq. 1b is

\begin{aligned} C_{1} (t_{i}) - {\hat{k}}_{2} e^{- {\hat{k}}_{2} t_{i}} \otimes C_{1} (t_{i}) = & K_{1} e^{- {\hat{k}}_{2} t_{i}} \otimes C p (t_{i}) \\ - k_{2} e^{- {\hat{k}}_{2} t_{i}} \otimes C_{1} (t_{i}) \end{aligned}

(2a)

where ⊗ denotes convolution and C₁ (t) is the measured tissue tracer concentration. This can be written in matrix form as

R = Z θ

where $θ^{T} = [K_{1}, k_{2}]$ and

\begin{aligned} R & = [\begin{matrix} C_{1} (t_{1}) - {\hat{k}}_{2} \int_{0}^{t_{1}} C_{1} (t) e^{- {\hat{k}}_{2} (t_{1} - t)} d t \\ ⋮ \\ C_{1} (t_{n}) - {\hat{k}}_{2} \int_{0}^{t_{n}} C_{1} (t) e^{- {\hat{k}}_{2} (t_{n} - t)} d t \end{matrix}] \\ Z & = [\begin{matrix} \int_{0}^{t_{1}} C p (t) - e^{- {\hat{k}}_{2} (t_{1} - t)} d t, \int_{0}^{t_{1}} C_{1} (t) e^{- {\hat{k}}_{2} (t_{1} - t)} d t \\ ⋮ \\ \int_{0}^{t_{n}} C p (t) - e^{- {\hat{k}}_{2} (t_{n} - t)} d t, \int_{0}^{t_{n}} C_{1} (t) e^{- {\hat{k}}_{2} (t_{n} - t)} d t \end{matrix}] \end{aligned}

(2b)

The solution is

\hat{θ} = (Z^{T} Z)^{- 1} Z^{T} R .

The parameter ${\hat{K}}_{2}$ is an initial estimate which can be obtained from the linear analysis of the noisy data, Eq. 1b or from the intercept of the graphical analysis (Logan et al., 1990). The solution estimates K₁ and a new k₂ which can then be used in the next iteration as ${\hat{K}}_{2}$ . This procedure provides unbiased estimates for K₁ and k₂ within one or two iterations (Feng et al., 1993).

Because the one-compartment model is somewhat restrictive, the above procedure was adapted to be applicable to a general PET data set and does not require a particular model structure. This is to be accomplished by applying the GLLS method for one-compartment to the data in two parts-that is, determine one set of parameters for times 0 to T₁ and a second set of parameters for T₁ to the end. The equation analogous to Eq. 2a for times ti > T₁ is

\begin{array}{l} C_{1} (t_{i}) = \exp (- {\hat{k}}_{2} (t_{i} - T_{1})) C_{1} (T_{1}) - {\hat{k}}_{2} \exp (- {\hat{k}}_{2} t_{i}) \int_{T_{1}}^{t_{i}} \exp ({\hat{k}}_{2} t) C_{1} (t) d t \\ = K_{1} \exp (- {\hat{k}}_{2} t_{i}) \int_{T_{1}}^{t_{i}} \exp ({\hat{k}}_{2} t) C p (t) d t - k_{2} \exp (- {\hat{k}}_{2} t_{i}) \int_{T_{1}}^{t_{i}} \exp ({\hat{k}}_{2} t) C_{1} (t) d t \end{array}

(2c)

Using the two sets of ĥ to generate a “smoothed” time-activity data set, the simple GLLS model of Eq. 2 will be tested to see if it can approximate data from more complex models. Applying the graphical method for reversible tracers to the smoothed data set should generate an unbiased estimate of the distribution volume (DV), if this smoothing technique is valid.

This procedure can be extended to the graphical analysis with a reference region. Using for the reference region

\frac{d C_{R E F}}{d t} = K_{1}^{R E F} C p - k_{2}^{R E F} C_{R E F}

(3)

and substituting for Cp(t) from Eq. 3 gives

\begin{aligned} C_{1} (t) & = K_{1} e^{- k_{2} t} \otimes C p (t) \\ = \frac{K_{1}}{K_{1}^{R E F}} C_{R E F} (t) - \frac{K_{1}}{K_{1}^{R E F}} (k_{2} - k_{2}^{R E F}) \int_{0}^{t} e^{- k_{2} (t - τ)} C_{R E F} (τ) d τ \end{aligned}

(4)

which is the form of the equation used by Lammertsma (1996) in which the k₂ here is in place of k₂a = k₂/(1 + BP) in that article. Expressing Eq. 4 in the form of Eq. 2 gives

\begin{aligned} C_{1} (t) = & {\hat{k}}_{2} e^{- {\hat{k}}_{2} t} \otimes C_{1} (t) = - k_{2} e^{- {\hat{k}}_{2} t} C_{1} (t) \\ + \frac{K_{1}}{K_{1}^{R E F}} C_{R E F} (t) - \frac{K_{1}}{K_{1}^{R E F}} (k_{2} - k_{2}^{R E F}) e^{- k_{2} t} \otimes C_{R E F} (t) \end{aligned}

(5)

which allows for an iterative linear solution for k₂ given an initial guess ( ${\hat{K}}_{2}$ ). There are now two other constants to be determined in addition to k₂, K₁/K₁^REF, and (K₁/K₁^REF)(k₂ -k₂^REF). In this case the solution designated as θ^T_REF, is given by θ^T_REF = [k₂, K₁/K^REF₁, (K₁/ K^REF₁)(k₂ - k^REF₂)]. The matrix Z of Eq. 2b now contains three columns corresponding to the three terms on the righthand side of Eq. 5. The constants can be used to generate a fit to the data at times t through the righthand side of Eq. 4. This can be applied in two time intervals as described previously, generating a fit to the data (C₁(t)) to which the graphical method can be applied. Alternatively, if the simple model is applicable-that is, the same parameters are used throughout the whole scanning period-the DVR (the ratio of the DV of the region with specific binding to that of the reference region) is given directly by the constants so that

D V R = \frac{K_{1} / k_{2}}{K_{1}^{R E F} / k_{2}^{R E F}}

(6)

as in the simplified reference tissue model (Lammertsma and Hume, 1996). Both of these techniques are explored.

MATERIALS AND METHODS

The above procedures have been applied to data simulated from the two-tissue compartment model given by

\begin{array}{l} C p ⇄_{k_{2}}^{K_{1}} C_{1} ⇄_{k_{4}}^{K_{3}} C_{2} \\ \frac{d C_{1}}{d t} = K_{1} C p (t) - (k_{2} + k_{3}) C_{1} + k_{4} C_{2} \\ \frac{d C_{2}}{d t} = k_{3} C_{1} - k_{4} C_{2} \end{array}

(7)

using measured plasma input functions, Cp(t). The radioactivity corresponding to a simulated region of interest (ROI) is given by C₁ + C₂. Two examples were used for these simulations which are similar to data from studies with RAC and MP and the plasma input functions were taken from actual PET studies. The model constants used for the RAC simulation data were K₁ = 0.15 mL/min/mL, k₂ = 0.36 min ¹, k₃ = 0.18 min⁻¹, and k₄ = 0.05 min⁻¹ for a DV = 1.92 mL/mL. For the reference region simulations, a one-tissue compartment model was used with K₁ = 0.08 mL/min/mL and k₂ = 0.192 min⁻¹. For the MP simulated data, K₁ = 0.6 mL/min/mL, k₂ = 0.06 min⁻¹, k₃ = 0.5 min⁻¹, and k₄ = 0.2 min⁻¹ for a DV = 35.0 mL/mL For the reference tissue, the model parameters were K₁ = 0.62 mL/min/mL and k₂ = 0.062 min⁻¹.

For comparison with the simulated data, human PET data from one RAC and one MP study were used. The PET scanner was a Siemens HR+, 63 slices, 4.5 × 4.5 × 2.4 mm full width half maximum (FWHM) in three-dimensional mode (see Wang et al. (2000) for experimental details using human subjects on the HR+). An arterial input function was measured for both tracers and plasma samples were analyzed for unchanged tracer (see Ding et al. (1997) and Volkow et al. (1993) for plasma sampling protocols). For the RAC study, 7.69 mCi was injected and 7.3 mCi was injected for the MP study. The scanning times were the same as used in the simulations. To increase the signal on each plane, contiguous planes were summed giving images with 4.8 mm FWHM. For both RAC and MP, one plane at the level of the basal ganglia was analyzed combining right and left. This gave a region containing 193 pixels for RAC and 237 pixels for MP. Data were analyzed both as the average of the pixels, which was taken to be the “true” value, and as the individual pixels.

There are a number of sources of noise in the PET image. The noise model used here is random, related to counting statistics. Noise is increased in the later time frames that have lower count rates because of radioactive decay. Taking into account scan times and the half-life of the isotope, the following formula was used to introduce random noise

d e v (t) = (x x) S c \sqrt{\frac{e^{- λ t} R O I (t)}{Δ t}}

(8)

where ROI(t) is the original (noise free) simulated radioactivity (C₁ + C₂) at time t, e^−λt ROI(t) refers to radioactivity not corrected for decay, e^λt refers to the decay correction where λ is the half-life of the isotope which is either 110 or 20.4 minutes for the simulation of ¹⁸F or ¹¹C. xx is a pseudo random number from a gaussian distribution with zero mean and variance of one, Sc is a scale factor that determines the level of noise, and Δt is the scan length. The simulated ROI with noise (ROI_N) is given by

R O I_{N} (t) = R O I (t) + e^{λ t} d e v (t)

where ROI(t) corresponds to either RAC or MP and e^λt dev(t) is the noise contribution at time t. The effect of introducing the decay factor is to take into account the difference in noise properties of the PET isotopes ¹¹C and ¹⁸F with their different half-lives. Although the examples used here are both from ¹¹C studies, this provides a mechanism for determining the effects of temporal differences in noise levels because the later scans of a ¹¹C study will have significantly more noise relative to the earlier scans than an ¹⁸F study given the same scanning schedule. The amount of noise introduced in the simulations depends upon both the scale factor, Sc, and the half-life used in the decay factor. For the examples presented here, Sc ranged from 0.25 to 8. The average noise in the different sets of simulations is compared by

< f_{N}^{R O I} >_{N, t} = a v g \{\frac{a b s (R O I_{N} (t) - R O I (t))}{R O I (t)}\}

(9)

where ROI(t) is the original (“true”) data, ROI_N(t) is the data with noise, and avg is the average over all time points (t) and over all data sets (N). A measure of how well the estimated ROI data approximates the original data is given by

< f_{N}^{E S T} >_{N, t} = a v g \{\frac{a b s (R O I_{E S T}^{N} (t) - R O I (t))}{R O I (t)}\}

(10)

where ROI^N_EST is the “smoothed” data set.

The scanning sequences used were the same as for the PET experiments. For RAC this was a total of 20 frames, 10 frames of 60 seconds and 10 frames of 300 seconds for 20 frames. For MP there were 22 frames, 10 × 60 seconds, 4 × 300 seconds, and 8 × 450 seconds. Data were analyzed graphically taking the average of slopes between 25, 30, to 75 and 25, 30, to 79 minutes for MP and 23, 28, to 53 and 23, 28, to 58 minutes for RAC. Region of interest estimates were determined using the one-tissue compartment fit for the entire time course and by the two part analysis. For the two part analysis, T₁ was taken to be 10 minutes. However, to assure a smooth transition between the two regions, the fit to the first part was taken to be somewhat longer, 15 to 18 minutes, and the second part was begun at 8 minutes. The choice of 8 minutes insured a sufficient number of points in each region and also times for graphical analysis were all included in the second region. The DVs were redetermined using the smoothed data. Whether the estimated ROI data used for the analysis was from the two part or one part fit (the original Feng formulation) depended upon which gave the best fit to the (noisy) data set in a least squares sense. The DV was also determined by optimizing the model parameters for the two-compartment model for RAC and for the one-compartment model for MP. The starting parameters were taken to be the “true” values ± a random amount of 25% of the parameter value. The DV based on the optimized parameters is given by (K₁/k₂)(1 + k₃/k₄) for RAC and K₁/k₂ for MP. The one-compartment model was sufficient to recover the true distribution volume for the simulated MP data.

The DVR was also calculated using the reference region as described in Logan et al. (1996) and using the simplified reference model of Lammertsma et al. (1996) adapted to the linear iterative form of Eq. 5, which requires only an initial estimate of k₂ designated ${\hat{K}}_{2}$ ; the 2 other parameters and subsequent values of k₂ are determined directly from the linear solution of Eq. 5. Similar to Gunn et al. (1997), the initial estimate, ${\hat{K}}_{2}$ , was taken to be -1/intercept from the graphical analysis, or if this was negative as occasionally occurs with very noisy data, an average ${\hat{K}}_{2}$ based on the original data was used. The DVR was determined graphically from the data with simulated noise using the average values of k₂^REF, 0.165 min⁻¹ (RAC), and 0.05 min⁻¹ (MP) (Logan et al., 1996). Data were “smoothed” and the estimated ROI data then were used to determine the DVR graphically. In the case of MP, which can be described by a one-compartment model, k₂^REF determined from simplified reference model using the noise free data for both regions was also used in the DVR calculation. As above, the ROI estimate was based on either the one or two part fit, whichever gave the best fit to the noisy data set {ROI_N (t)}. The DVR was also calculated from the parameters of the one part fit (Eq. 6).

Solution of the differential equations and optimization of model parameters were accomplished using routines in Numerical Recipes in C (Press et al., 1988). Numerical integrations of the convolutions in Eqs. 2, 4, and 5 were performed using the routine “qtrap” also from Numerical Recipes in C (Press et al., 1988). A linear interpolation between time points was used for the ROI data. Equal weights were used in all of the linear analyzes.

RESULTS

Results from the analyzes of the MP simulated data using the ¹¹C and ¹⁸F factors in noise generation are given in Tables 1A and 1B, respectively. Plots are also shown of the mean and standard deviation of DV_G, DV_F, and DV_FG for ¹¹C in Fig. 1A. The average values over the 500 data sets of f_N^ROI>_t, <f_N^EST>_t, the graphically determined DV, DV_G, are reported as well as DV_FG, which refers to the DV determined from the smoothed ROI data (ROI_EST), whether by two part smoothing operation or from the one-compartment fit, whichever gave the best fit to ROI_N. (At the higher noise levels this was generally the two part fit because of its greater flexibility.) Also given are DV_F, the DV determined using the one-compartment solution of Eq. 2a (Feng's method), and DV_NLS, determined from the model fit using the differential equations directly. The corresponding coefficients of variation (CV) reported as the standard deviation/mean are given in parentheses. As expected, DV_G decreases with increasing noise. From Table 1A, (¹¹C) DV_G shows a 18% decrease in the mean as f_N^ROI>_N,t increases from 0.008 to 0.11. At the highest noise level tested (f_N^ROI>_N,t =0.26), the decrease in DV_G was 40%. Similar results were observed by (Abi-Dargham et al., 2000) in simulations based on [¹¹C]NNC 112. The DV generated using the 2 part fit (DV_FG) shows an increase in average value with increasing noise, 10% over the true value at the highest noise level. Both DV_F and DV_NLS based on the one-compartment model gave the same results with a much smaller increase in both the DV and coefficient of variation. At intermediate noise levels (f_N^ROI> = 0.10), DV_FG, DV_F, and DV_NLS give similar results. For all methods, the CV increases with noise but the 2 part fit has a somewhat larger CV because of the increased flexibility of the fit. Also, at the highest noise levels, the GLLS method fails in some cases in the determination of the solution vector, θC . For Sc = 8, 35 out of 500 DV determinations were either out of range (<0 or >100) or failed to generate a solution θC . For the ¹⁸F decay factor (Table 1B), the average value for DV_FG changes very little with increasing noise and the CV is smaller than that for ¹¹C noise, basically the same as for DV_F and DV_NLS. The decrease in DV_G with increasing noise is also less, a decrease of 11% with increase in f_N^ROI>_N,t from 0.023 to 0.14. When f_N^ROI>_N,t is on the order of 0.14, the CVs for the ¹¹C simulations are 50% greater than those of the ¹⁸F simulations (0.23 and 0.09, respectively for DV_FG). The ROI estimates, {ROI^N_EST(t)}, generated by the smoothing technique provide a better approximation to the original data than do the noisy data sets, {ROI_N(t)}, as measured by the ratio f_N^EST>^Nt/(f_N^ROI>^N_t which is on the order of 0.4 to 0.5 for all noise levels in both simulations.

FIG. 1.

(A) Plots of data from Table 1A (MP), DV versus f_N^ROI>: DV_G (⋄), DVF (▲), DV_FG (♦). (B) Plots of data from Table 2A (RAC), DV versus f_N^ROI>: DV_G (⋄), DV_FG (♦). MP, [¹¹C]-d-threo-methylphenidate; DV, distribution volume; RAC, ¹¹C raclopride.

TABLE 1.

Distribution volumes for MP simulated data with “true” DV of 35 using ¹¹C and ¹⁸F noise from Eq. 8 with Sc as shown

A. ¹¹C
Sc	<f_N^ROI>	<f_N^EST>	DV_G	DV_FG	DV_F	DV_NLS
0.25	0.0082	0.005	34.9 (0.012)	35.0 (0.01)	34.8 (0.007)	34.8 (0.007)
0.50	0.016	0.007	34.8 (0.024)	35.0 (0.02)	34.8 (0.03)	34.8 (0.015)
1.00	0.032	0.0135	34.3 (0.046)	35.1 (0.04)	34.8 (0.044)	34.8 (0.030)
1.50	0.056	0.019	33.5 (0.072)	35.2 (0.056)	34.9 (0.059)	34.9 (0.043)
2.0	0.064	0.026	32.5 (0.093)	35.3 (0.081)	34.9 (0.074)	34.9 (0.059)
2.5	0.080	0.0308	31.3 (0.109)	35.5 (0.097)	35.1 (0.09)	35.0 (0.072)
3.0	0.096	0.037	30.1 (0.12)	35.7 (0.12)	35.2 (0.11)	35.0 (0.091)
3.5	0.11	0.043	28.8 (0.14)	36.0 (0.15)	35.3 (0.12)	35.2 (0.11)
4.0	0.13	0.051	27.6 (0.15)	36.6 (0.194)	35.3 (0.12)	35.3 (0.12)
4.5	0.14	0.057	26.4 (0.16)	37.1 (0.23)	35.4 (0.16)	35.4 (0.14)
5.0	0.16	0.063	25.4 (0.16)	37.5 (0.26)	35.6 (0.16)	35.5 (0.16)
6.0	0.19	0.075	23.7 (0.17)	38.6 (0.34)*	35.9 (0.20)	35.9 (0.20)
8.0	0.26	0.096	21.3 (0.19)	38.5 (0.35)*	37.7 (0.40)	37.9 (0.40)

B. ¹⁸F

Sc	<f_N^ROI>	<f_N^EST>	DV_G	DV_FG	DV_F	DV_NLS

1.00	0.023	0.009	34.83 (0.02)	35.05 (0.014)	34.78 (0.013)	34.79 (0.013)
1.50	0.034	0.013	34.66 (0.03)	35.07 (0.022)	34.77 (0.020)	34.80 (0.021)
2.0	0.046	0.017	34.42 (0.04)	35.1 (0.029)	34.80 (0.026)	34.81 (0.027)
2.5	0.057	0.021	34.14 (0.049)	35.1 (0.036)	34.82 (0.033)	34.85 (0.034)
3.0	0.068	0.025	33.80 (0.059)	35.1 (0.043)	34.83 (0.040)	34.83 (0.040)
3.5	0.080	0.029	33.4 (0.067)	35.2 (0.051)	34.86 (0.046)	34.85 (0.047)
4.0	0.091	0.033	32.9 (0.075)	35.24 (0.058)	34.89 (0.054)	34.90 (0.053)
4.5	0.10	0.037	32.44 (0.083)	35.5 (0.066)	34.95 (0.061)	34.95 (0.061)
5.0	0.11	0.041	31.9 (0.091)	35.4 (0.074)	35.0 (0.068)	34.97 (0.069)
6.0	0.14	0.050	30.77 (0.106)	35.5 (0.09)	35.1 (0.084)	35.07 (0.085)

Mean for DV_FG was determined by eliminating those for which the method failed or had values <0 and >100, that were clearly out of range (35 out of 500 data sets for SC = 8 and 8 out of 500 for SC = 6). Plots of DV_G, DV_F, and DV_FG with their standard deviations are also shown in Fig. 1A. Note the larger standard deviations associated with the ¹¹C simulations in A.

Figure 2 (♦) illustrates simulated data (MP) generated with noise based on ¹¹C. For these simulations f_N^ROI>_t = 0.10 and 0.095 for Fig. 2A and 2B, respectively. DV_G was 43 and 30 mL/min/mL for Fig. 2A and 2B, respectively, whereas the smoothed data gave 44 and 35.9. The smoothed data in Fig. 2A differed from the “true” data significantly at later time points although still providing a good fit to the simulated data. In Fig. 2B, the smoothed data provides a good fit to the “true” data. Distribution volumes from the 1-compartment model (DV_NLS) were 38.8 and 36.7 mL/mL.

FIG. 2.

Simulated time-activity data for example [¹¹C]-d-threo-methylphenidate (MP). Original data (data without added noise) indicated by the solid line has distribution volume (DV) = 35 mL/min/mL. Data with noise added and ¹¹C decay factor is indicated by ♦. The estimated region of interest (ROI) data is indicated by diam;. In both cases, the smoothed data appears to describe the noisy data but in A it differs from the “true” data. (A) DV_G = 43.3, DV_FG = 44.0, DV_F = 38.9, and DV_NLS = 38.8 mL/mL and f_N^ROI> = 0.10. (B) DV_G = 30.00, DV_FG = 35.9, DV_F = 36.7, and DV_NLS = 36.7 mL/mL and f_N^ROI> = 0.09.

Summaries of the data for RAC are given in Tables 2A and 2B. At the lowest noise level, the smoothing operation actually introduces more noise because f_N^EST>_N,t is greater than f_N^ROI>_N,t (0.012 and .035, respectively). The RAC data indicate a 26% decrease from the true value in DV_G for f_N^ROI>_N,t = 0.21, the highest noise level for ¹¹C (Table 2A). For comparable noise levels, as measured by f_N^ROI>_N,t, DV_G is 13% less for ¹¹C and 6.4% less for ¹⁸F (Table 2B) (f_N^ROI>_N,t = 0.12). This indicates the greater sensitivity to noise in the later scans (because of ¹¹C vs. ¹⁸F) as was also observed with MP. DV_FG values from both sets of noise simulations show an increase with an increase in noise. The DV_F in both cases underestimates the DV by approximately 10%. The CVs for all measures of DV increase with increasing noise. As with MP, the method fails more frequently as the noise increases. The DV_NLS increases somewhat with increasing noise. To test the effect of varying the starting values of the model parameters, DVs were compared using both the true and the randomly displaced as initial values in the optimization procedure in the direct solution of the differential equations. For Sc = 2, there was a small improvement in using the “true” values, DV_NLS = 1.93(0.08) mL/mL; but for Sc = 3.5, DV_NLS = 1.99(0.22) mL/mL, the same as when the initial values were displaced from the true values.

TABLE 2.

Distribution volumes for RAC simulated data with “true” DV of 1.92 using ¹¹C and ¹⁸F noise from Eq. 8 with Sc as shown

A. ¹¹C

Sc	<f_N^ROI>	<f_N^EST>	DV_G	DV_FG	DV_F	DV_NLS
0.25	0.012	0.035	1.88 (0.014)	1.93 (0.013)	1.71 (0.007)	1.92 (0.025)
0.40	0.019	0.036	1.88 (0.022)	1.93 (0.021)	1.71 (0.011)	1.92 (0.027)
0.60	0.028	0.034	1.87 (0.032)	1.92 (0.031)	1.71 (0.028)	1.93 (0.031)
0.80	0.037	0.036	1.86 (0.043)	1.92 (0.042)	1.71 (0.038)	1.92 (0.038)
1.00	0.046	0.038	1.84 (0.052)	1.93 (0.053)	1.72 (0.038)	1.92 (0.044)
1.20	0.056	0.040	1.82 (0.061)	1.93 (0.065)	1.72 (0.033)	1.93 (0.06)
1.60	0.074	0.045	1.77 (0.079)	1.95 (0.092)	1.72 (0.043)	1.94 (0.09)
2.0	0.093	0.051	1.71 (0.092)	1.97 (0.12)	1.72 (0.059)	1.95 (0.105)
2.5	0.12	0.059	1.64 (0.106)	2.00 (0.16)	1.74 (0.070)	1.97 (0.13)
3.0	0.14	0.067	1.56 (0.18)	2.05 (0.20)	1.73 (0.085)	1.96 (0.15)
3.5	0.16	0.075	1.50 (0.25)	2.07 (0.26)*	1.74 (0.10)	1.99 (0.21)
4.5	0.21	0.093	1.40 (0.14)	2.06 (0.31)*	1.75 (0.13)*	2.05 (0.28)

B. ¹⁸F

Sc	<f_N^ROI>	<f_N^EST>	DV_G	DV_FG	DV_F	DV_NLS

0.50	0.023	0.025	1.88 (0.015)	1.90 (0.013)	1.71 (0.008)	1.91 (0.025)
0.75	0.030	0.027	1.88 (0.021)	1.90 (0.020)	1.71 (0.012)	1.92 (0.033)
1.00	0.040	0.037	1.87 (0.028)	1.93 (0.026)	1.71 (0.015)	1.92 (0.033)
1.25	0.050	0.036	1.86 (0.035)	1.92 (0.032)	1.71 (0.019)	1.92 (0.035)
1.40	0.056	0.037	1.86 (0.039)	1.92 (0.036)	1.71 (0.0211)	1.92 (0.040)
1.60	0.064	0.039	1.85 (0.043)	1.92 (0.042)	1.71 (0.021)	1.92 (0.045)
1.80	0.072	0.043	1.84 (0.049)	1.94 (0.048)	1.71 (0.025)	1.93 (0.051)
2.0	0.080	0.044	1.83 (0.054)	1.94 (0.053)	1.71 (0.028)	1.94 (0.06)
2.5	0.10	0.049	1.80 (0.065)	1.95 (0.068)	1.71 (0.036)	1.94 (0.08)
3.0	0.12	0.054	1.76 (0.076)	1.96 (0.083)	1.72 (0.047)	1.95 (0.091)
3.5	0.13	0.058	1.73 (0.08)	1.95 (0.096)	1.72 (0.05)	1.95 (0.095)

(f_N^ROI) and (f_N^EST), defined in Eqs. 9 and 10, are measures of the difference between the original data and the data with noise and between the original data and the “smoothed” data, respectively. DV_G is the graphical distribution volume, DV_FG is the DV calculated using the smoothing technique described before applying the graphical analysis, DV_F is the DV determined using Feng's method for the 1-compartment model (Eqs. 2a and 2b), and DV_NLS is determined from the one-compartment model directly by solving the differential equation and optimizing the model parameters. The coefficients of variation (CV) (standard deviation/mean) are given in parentheses. Plots of DV_G and DV_FG with their standard deviations are also shown in Fig. 1. Note the larger standard deviations associated with the ¹¹C simulations. *Mean for DV_FG was determined by eliminating those data sets for which the method failed or had values <0 or >12, that were clearly out of range (3 out of 500 data sets for Sc = 3.5 and 15 out of 500 for Sc = 4.5).

Figure 3 illustrates examples from RAC (♦) generated with noise based on ¹¹C (f_N^ROI>t = 0.10). In Fig. 3A, the smoothed data underestimates the “true” data and in Fig. 3B, the smoothed data is close to the “true” data. DV_G is 1.58 (Fig. 3A) and 1.70 (Fig. 3B) mL/mL. The smoothed data (diam;) has f_N^EST>t = 0.051 (Fig. 3A) and 0.034 (Fig. 3B). DV_FG = 1.74 (Fig. 3A) and 1.90 (Fig. 3B) mL/min/mL and DV_NLS from the 2-compartment model was 1.73 (Fig. 3A) and 1.85 (Fig. 3B) mL/mL. DV_F = 1.61 (Fig. 3A) and 1.71 (Fig. 3B) mL/mL.

FIG. 3.

Simulated time-activity data for example ₁₁C raclopride (RAC). Original data (data without added noise) indicated by the solid line has DV = 1.92 mL/mL. Data with noise added and ₁₁C decay factor is indicated by ♦. The estimated region of interest (ROI) data is indicated by ⋄. In both cases the smoothed data appears to describe the noisy data but in A it differs from the “true” data. DV_G = 1.58, DV_FG = 1.74, DV_F = 1.61, and DV_NLS = 1.73 mL/mL and f_N^ROI> = 0.103. (B) DV_G = 30.00, DV_FG = 1.90, DV_F = 1.71, and DV_NLS = 1.85 mL/mL and f_N^ROI> = 0.099.

Tables 3 and 4 present summaries of simulations of DVR calculations using reference tissue data. Average values over all 500 data sets of DVR_G with the average k₂ (K₂^-REF) (Logan et al., 1996), DVR_FL calculated from the parameters of the simplified reference tissue model (Lammertsma and Hume, 1996), DVR_FG calculated graphically using the smoothed data, and k₂^-REF and DVR_FG calculated graphically without k₂^-REF. For MP, DVR_FG calculated using the k₂^REF determined from the reference region and the “true” data was also reported (last column).

TABLE 3.

Distribution volume ratios for MP simulated data with “true” DVR of 3.5 using ¹¹C noise from Eq. 8 with Sc as shown

Sc	<f_N^EST>	DVR_Gk̄₂^REF	DVR_FL	DVR_FGk̄₂^REF	DVR_FG (no k₂)	DVR_FGk₂^REF
0.25	0.006	3.67 (0.014)	3.53 (0.014)	3.68 (0.012)	2.87 (0.005)	3.53 (0.015)
0.50	0.008	3.65 (0.028)	3.53 (0.030)	3.69 (0.024)	2.88 (0.011)	3.53 (0.030)
1.00	0.015	3.58 (0.056)	3.54 (0.061)	3.68 (0.050)	2.87 (0.061)	3.53 (0.051)
1.50	0.021	3.48 (0.082)	3.58 (0.098)	3.69 (0.074)	2.85 (0.032)	3.52 (0.07)
2.00	0.028	3.35 (0.104)	3.63 (0.14)	3.71 (0.10)	2.88 (0.043)	3.54 (0.09)
2.50	0.035	3.21 (0.122)	3.71 (0.20)	3.74 (0.13)	2.89 (0.054)	3.56 (0.11)
3.00	0.043	3.05 (0.14)	3.86 (0.31)	3.77 (0.16)	2.89 (0.065)	3.58 (0.14)
3.50	0.050	2.90 (0.14)	3.82 (0.29)	3.75 (0.17)	2.89 (0.062)	3.61 (0.16)
4.00	0.063	2.75 (0.16)	3.94 (0.37)*	3.79 (0.20)	2.89 (0.08)	3.64 (0.19)
4.50	0.071	2.62 (0.17)	3.88 (0.34)*	3.92 (0.27)	2.90 (0.09)	3.68 (0.22)
5.00	0.078	2.51 (0.18)	3.81 (0.34)*	3.94 (0.30)	2.90 (0.10)	3.69 (0.24)
6.00	0.094	2.30 (0.19)	3.77 (0.37)*	4.00 (0.38)	2.92 (0.12)	3.76 (0.32)
8.00	0.125	2.11 (0.20)	4.82 (1.1)*	3.93* (0.71)	2.92 (0.17)	3.78* (0.57)

(f_N^EST), defined in Eq. 10, is a measure of the difference between the original data without noise and the “smoothed” data. DVR_G is the graphical DVR using a reference ROI and an average $\overset{ˉ}{K}$ ₂^REF = 0.05 min⁻¹ (“true” k₂ = 0.063 min⁻¹). DV_FL is the DVR determined using Feng's method for the 1-compartment model adapted to use the reference region. DVR_FG is the DVR calculated using the smoothing technique before applying the graphical analysis with the reference region, and $\overset{ˉ}{K}$ ₂^REF and DVR_FG (no k₂) is determined using the reference region graphical analysis on the smoothed data (without k₂). DVR_FG k₂^REF is the DVR determined graphically using the smoothed data and k₂ determined from the original data and the cerebellum (k₂^REF = 0.061 min⁻¹). The coefficients of variation are given in parentheses.

Indicates that values out of range (<0 or >10) were eliminated.

TABLE 4.

Distribution volume ratios for RAC simulated data with “true” DVR of 4.6 using ¹¹C noise from Eq. 8 with Sc as shown

Sc	<f_N^EST>	DVR_Gk̄₂^REF	DVR_FL	DVR_FGk̄₂^REF	DVR_FG (no k₂)
0.25	0.026	4.55 (0.014)	5.02 (0.020)	4.61 (0.020)	4.40 (0.010)
0.40	0.027	4.54 (0.023)	5.03 (0.032)	4.64 (0.021)	4.42 (0.016)
0.60	0.029	4.52 (0.034)	5.05 (0.049)	4.66 (0.039)	4.44 (0.024)
0.80	0.031	4.48 (0.044)	5.08 (0.067)	4.69 (0.050)	4.42 (0.032)
1.00	0.033	4.43 (0.054)	5.11 (0.088)	4.71 (0.060)	4.42 (0.040)
1.20	0.036	4.38 (0.063)	5.15 (0.112)	4.73 (0.070)	4.43 (0.049)
1.60	0.042	4.26 (0.080)	5.28 (0.200)	4.79 (0.092)	4.46 (0.067)
2.00	0.050	4.11 (0.094)	5.42 (0.70)	4.83 (0.108)	4.48 (0.079)
2.50	0.061	3.93 (0.109)	5.47 (1.0)	4.83 (0.16)	4.49 (0.10)
3.00	0.071	3.75 (0.12)	6.63 (1.2)	4.92 (0.20)	4.59 (0.14)
3.5	0.078	3.59 (0.13)	5.7 (1.9)	5.09 (0.26)	4.63 (0.17)
4.5	0.10	3.22 (0.14)	5.8 (0.48)*	5.1 (0.28)*	4.76 (0.22)

(f_N^EST> is defined in Eq. 10. DVR_G is the graphical DVR using a reference ROI and an average $\overset{ˉ}{K}$ ₂^REF = 0.165 min⁻¹ (“true” k₂ = 0.19 min⁻¹). DVR_FL is the DV determined using Feng's method for the 1-compartment model adapted to use the reference region. DVR_FG is the DVR calculated using the smoothing technique before applying the graphical analysis with the reference region and $\overset{ˉ}{K}$ ₂^REF and DVR_FG (no k₂) is determined using the reference region graphical analysis on the smoothed data (without k₂). The coefficients of variation are given in parentheses.

Indicates that values out of range (<0 and >10) were eliminated (40 out of 500).

For MP (Table 3), the DVR_G is somewhat larger (5%) than the true value at the lowest noise level. This is because k₂^-REF for the reference region differs from the actual k₂ used in the simulations. For the same reason DV_FG (k₂^-REF) also somewhat overestimates the DVR. DVR_FL provides a good estimate at low noise levels, but the CVs increase more rapidly with noise than the other methods. At the highest noise levels, the means and CVs reported in Table 3 (*) were obtained by eliminating values <0 and >10 (approximately 7% of the values at the highest noise level). For MP, which conforms to a one-tissue compartment model, the linearized simplified reference tissue method does provide the correct value of k₂^-REF (k₂^-REF = 0.062 ± 0.017 for Sc = 3, and 0.064 ± 0.027 for Sc = 6). Using 0.062 for k₂^REF gives the results reported in the last column of Table 3. Although DVR_FG (no k₂) underestimates the DVR, the CVs are less and there is very little variation in mean from the lowest to highest noise level.

Besides $\overset{ˉ}{K}$ ₂^REF, the other parameter estimates from the linearized simplified reference tissue model for MP also agree well with the “true” values. For Sc = 2, these are k₂>_N = 0.0159 min⁻¹ and K₁/K₁^REF − 0.92 compared with “true” values of 0.016 min⁻¹ and 0.98, respectively.

For RAC at the lowest noise level (Sc = 0.25, Table 4) all DVRs except DVR_FL are within 5% of the true value. The means of DVR_G, DVR_FG ( $\overset{ˉ}{K}$ ₂^REF), and DVR_FG (no k2) provide good estimates of the true DVR at the low to moderate noise levels. At higher noise levels DVR_G underestimates the DVR. The most stable value over the range of noise levels is DVR_FG (no k₂), which also has the smallest CVs. DVR_FL increases significantly (as does the CV) as the noise increases. This is because this model is not appropriate for RAC. The individual parameter estimates from the data set for Sc = 0.25 are k₂^REF>_N − 0.092 min⁻¹, k₂>_N = 0.024 min⁻¹, and K₁/K₁^REF 1 3, which are not in agreement with the “true” values of 0.192, 0.078, and 1.87 min⁻¹, respectively. Using the original data without noise, the 1-compartment reference tissue method gives similar results, k₂^REF = 0.092 min⁻¹, k₂ = 0.024 min⁻¹, and K₁/K₁^REF = 1.31 so that DVR_FL = 5.05, which is approximately 10% too high. The values for DVR_FG( $\overset{ˉ}{K}$ ₂^REF) and DVR_FG(no k₂) for RAC do not differ significantly indicating that the term containing $\overset{ˉ}{K}$ ₂^REF is not necessary. DVR_FG also fails in some instances at the higher noise levels.

Figures 4 and 5 illustrate data generated using the reference tissue method for MP and RAC. The solid line indicates the original data. The input data were taken to be the original ROI and reference region data to see how well the methods work without the effects of added noise. In Fig. 4, the fit generated using the linearized reference tissue method (⋄) reproduces the original data set very closely (DVR = 3.6 mL/mL). In Fig. 5, the ROI data generated using the 2 part fit described previously is indicated by ⋄ and the fit generated using the reference tissue method is indicated by ♦. The true DVR was 4.6 and the DVR_FG generated graphically from the 2 part fit was 4.62. The DVR calculated from the model parameters using the reference tissue method was 5.05 mL/mL, however, the DVR calculated graphically using the ROI estimated by the reference tissue method was 4.82 mL/mL. For DVR_FG, the average value of k₂ (0.165 min⁻¹) was used, although for this data neglecting k₂ introduces only a small error.

FIG. 4.

Data generated using the reference tissue method for [¹¹C]-d-threo-methylphenidate (MP) (⋄) from Eq. 4. The original data is indicated by the solid line. The input data were taken to be the original region of interest data and reference region data without noise. The reference tissue method reproduces the original data.

FIG. 5.

Data generated using the reference tissue method (without added noise) for ¹¹C raclopride (RAC) (♦). The distribution volume ration (DVR) based on the calculated parameters was 5.05 from Eq. 6. The solid line is the original data that has DVR = 4.6. The region of interest (ROI) data indicated by ⋄ represent the estimated ROI fit using the two part fit with the reference region input. The DVR generated physically using this data was 4.66 using the average k₂ and 4.5 without k₂, both of which are very close to the true value and better estimates than the original reference tissue method.

Table 5 presents the analysis of image data from a RAC study on the HR+. One hundred ninety-three pixels from one plane including right and left basal ganglia were used. In Table 5A, results from the ROI analysis of this data are given. Both 1- and 2-tissue compartment models were used to analyze the data from the basal ganglia (BG) and the cerebellum (CB), which was used as the reference region. For both regions, the two-tissue compartment gave a better fit to the data and agreed more closely with the graphical analysis. This was only a small difference in the BG, but the DV for the CB was 15% greater using the 2-compartment model. The binding potential (BP), defined as k₃/k₄, is 10% less than its estimation as DVR-1, using either the graphical or 2-compartment model for BG and CB DVs. Based on the 1-compartment model for CB, DVR-1 is 25% greater. The Feng method gave the same result as the NLS method for a one-compartment model. The DVR estimates were in close agreement for all methods. Results from the analysis of time-activity curves of the individual pixels are given in Table 5B. The average values of all methods were the same as for the ROI analysis except for DV_G, which was 3% less. The binding potential, k₃/k₄, was in agreement with ROI analysis although with a CV of 0.28. The DVR data was also in good agreement with the ROI data. Although it is not possible to know the “true” data, using the average ROI data as a measure of the “true” data, f_N^ROI>_N,t was found to be 0.11. A plot of data from one pixel is shown in Fig. 6, with the image data indicated by ♦ and smoothed data is indicated by ⋄ with the solid line. The DVs from the various methods are in good agreement (within 6%).

FIG. 6.

Data from a raclopride study representing one pixel at the level of the basal ganglia (♦). The smoothed data is indicated by ⋄. DV_G = 1.8, DV_FG = 1.91, DV_F = 1.83, and DV_NLS = 1.91 mL/mL.

TABLE 5A.

ROI analysis (193 pixels from basal ganglia)

Plasma input function
Method	DV (BG)	BP (k3/k4)	DV (CB)	DV (BG/CB)
NLS (2 parm)*	1.69		0.367	4.60
NLS (4 parm)†	1.78	2.78	0.435	4.10
DV_G	1.77		0.432	4.10
DV_FG	1.78
DV_F	1.69

Reference region input

Method				DVR
DVR_G (avg k₂)				4.04
DVR_G (no k₂)				3.95
DVR_FL				3.99
DVR_FG (avg k₂)				4.03
DVR_FG (no k₂)				3.95

TABLE 5B.

Pixel analysis for BG

Plasma input function
Method	DV (BG)	BP (k3/k4)	DV (CB)	DV (BG/CB)
NLS (2 parm)	1.69 (0.14)		0.367	4.60
NLS (4 parm)	1.79 (0.14)	2.83 (0.28)	0.435	4.11
DV_G	1.72 (0.14)		0.432	4.10
DV_FG	1.78 (0.14)
DV_F	1.69 (0.13)

Reference region input

Method				DVR
DVR_G (avg k₂)				3.94 (0.14)
DVR_G (no k₂)				3.91 (0.13)
DVR_FL				4.03 (0.13)
DVR_FG (avg k₂)				4.07 (0.13)
DVR_FG (no k₂)				3.98 (0.13)

A. Region of interest (ROI) analysis of the average over pixels in a region of interest from the basal ganglia (right and left combined) for a study with ¹¹C raclopride. B. Analysis of individual pixels within the region. ^ROI> = 0.11 for the difference between ROI(t) taken from the individual pixels and the average over pixels in the basal ganglia (BG). BP, binding potential; DV, distribution volume; DVR, distribution volume ratio, NLS, nonlinear least squares.

K₁ = 0.122, k2 = 0.2557, k3 = 0.217, k4 = 0.0784 min⁻¹.

†

K₁ = 0.083, k2 = 0.0496, DV = 1.69.

Table 6 presents the analysis of image data from a MP study on the HR+. Two hundred thirty-seven pixels from one plane including right and left basal ganglia were used. In this case, the image appears more noisy than the RAC study with f_N^ROI>_N,t = 0.18, also using the average time-activity curve as a measure of the “true” data. An example from this data set is shown in Fig. 7. In this example there was a small difference between DVs from the one- and two-compartment models for both BG and CB. The two-compartment model gave a slightly better fit to the data for the BG (the ratio of X² for two-compartment to the one-compartment model was 0.64). This ratio was 0.10 for CB. DV_G was closer to the 2-compartment value. Parameter values from the linearized reference tissue method were k₂^REF = 0.045, k₂ = 0.02 min⁻¹, and K₁/K₁^REF = 1.19, which is in excellent agreement with values from 1-compartment model (NLS), 0.046, 0.019 min⁻¹, and 1.16, respectively. In the pixel analysis, DV_F agrees well with NLS (2 parm) (23.5 and 22.8 mL/mL, respectively), as expected because they both represent a 1-compartment model. Similarly, DV_FG and NLS (4 parm) agree (25.1 and 25.6 mL/mL). The binding potential (k₃/k₄) calculated from the pixel analysis using the 2-compartment model has a large CV but is in reasonable agreement with DVR-1 based on DV_FG or the DVR from the model solutions. The binding potential calculated from the ROI analysis is somewhat larger. DV_G in the pixel analysis underestimates the DV by 20%. Similarly, DVR_G in the pixel analysis underestimates the DVR, whereas in the ROI analysis, DVR_G is in good agreement with the other methods. In the pixel analysis, DVR_FG using k₂ = 0.045 min⁻¹ is close to the ROI values. Without k₂, the DVR_FG underestimates the DVR by 20%, assuming the “true” value to be 1.8. DVR_FL is much more sensitive to noise with a large CV. This perhaps could be improved by using k₂^REF determined from the ROI analysis.

FIG. 7.

Data from a [¹¹C]-d-threo-methylphenidate (MP) study representing one pixel at the level of the basal ganglia (♦). The smoothed data is indicated by ⋄. DV_G = 20.4, DV_FG = 26.7, DV_F = 27.5, and DV_NLS = 26 mL/mL.

TABLE 6A.

ROI analysis (237 pixels from basal ganglia)

Plasma input function
Method	DV (BG)	BP (k₃/k₄)	DV (CB)	DVR (BG/CB)
NLS (2 parm)*	22.5		8.16	2.76
NLS (4 parm)†	24.4	2.28	8.85	2.76
DV_G	24.5		8.59	2.85
DV_FG	24.0			2.79
DV_F	22.3			2.60

Reference region input

Method				DVR
DVR_G (k₂ = 0.045)				2.74
DVR_G (no k₂)				2.37
DVR_FL				2.71
DVR_FG (k = 0.045)				2.71
DVR_FG (no k₂)				2.34

TABLE 6B.

Pixel analysis for BG

Plasma input function
Method	DV (BG)	BP (k3/k4)	DV (CB)	DVR (BG/CB)
NLS (2 parm)	23.5 (0.19)		8.16	2.88
NLS (4 parm)	25.6 (0.34)	1.93 (.99)	8.85	2.89
DV_G	18.7 (0.20)		8.59	2.18
DV_FG	25.1 (0.25)			2.92
DV_F	22.8 (0.18)			2.65

Reference region input

Method				DVR
DVR_G (k₂ = 0.045)				2.22 (.21)
DVR_G (no k₂)				2.21 (.16)
DVR_FL				2.65 (.67)
DVR_FG (k₂ = 0.045)				2.77 (.21)
DVR_FG (no k₂)				2.26 (.15)

A. Region of interest (ROI) analysis of the average over pixels in a ROI from the basal ganglia (right and left combined) for a study with ¹¹C d threo methylphenidate. B. Analysis of individual pixels within the region. <f_N^ROI> = 0.18 for the difference between ROI(t) taken from the individual pixels and the average over pixels in the basal ganglia (BG).

DV, distribution volume; BP, binding potential; CB, cerebellum; DVR, distribution volume ratio.

K₁ = 0.439, k₂ = 0.0186 min_-1.

†

K₁ = 0.484, k₂ = 0.0654, k₃ = 0.205, k₄ = 0.090 min⁻¹.

DISCUSSION

Graphical analysis of PET data has the advantage of being a rapid, easily applied method. In addition, it is model-independent and will yield the total distribution volume of a ROI whether it conforms to a 1-, 2-, or multi-tissue compartment model. However, the primary drawback is that it is subject to bias in the DV estimates when the time-activity data contain noise. This is less of a problem with small DVs and for ROI data; however, for the generation of DV images and in the case of noisy ROI data with higher DV values it would be useful to have a method that preserves the model-independent approach but estimates the total DV without a bias. To this end we have tested an adaptation of Feng's bias-free linear method for parameter estimation of the one-compartment model using simulation data based on the PET ligands, raclopride (DV = 1.92 mL/mL) and d-threo-methylphenidate (DV = 35 mL/mL), with varying amounts of noise. We also have performed some preliminary tests on image data using the same tracers for comparison. The kinetic parameters used in the simulations for these two ligands are very different illustrating two general examples. For MP, the binding constants, k₃ = 0.5 and k₄ = 0.2 min⁻¹, are rapid and the large value of the DV is due to the slower tissue to plasma efflux, k₂. In the other example (RAC), k₂ is large (0.36 min⁻¹) but the ligand receptor dissociation constant k₄ in this example is considerably smaller (0.05 min⁻¹). A result of these different kinds of kinetics is that a 2-compartment model is required to describe RAC, but an accurate measure of the total DV can be obtained from a 1-compartment model for MP. Why 1 set of parameters requires a 2-compartment model and the other does not has to do with the impulse response function of the 2-compartment model given by A₁ exp(-α₁t) + A₂ exp(-α₂t) (Carson et al., 1998), where A_1,2 and α_1,2 are combinations of K₁, k₂, k₃, and k₄,

\begin{aligned} α_{1, 2} & = (γ \mp \sqrt{γ^{2} - 4 k_{2} k_{4}}) / 2 γ = k_{2} + k_{3} + k_{4} \\ A_{1} & = \frac{K_{1}}{α_{1} - α_{2}} (k_{3} + k_{4} - α_{1}) \\ A_{1} & \frac{K_{1}}{α_{1} - α_{2}} (α_{2} - k_{3} - k_{4}) \end{aligned}

If one of the exponential terms dominate, a one-compartment model will adequately describe the data. Following (Carson et al., 1998) whether a two-compartment fit is required can be determined by considering the fraction of the area of the response function due to the second term for time T, that is

\frac{A_{2} (I - \exp (- α_{2} T)) / α_{2}}{A_{1} (1 - \exp (- α_{1} T)) / α_{1} + A_{2} (1 - \exp (- α_{2} T)) / α_{2}}

(6)

For RAC (simulated data), the ratio of Eq. 6 is 0.38 at 5 minutes, 0.168 at 20 minutes, 0.13 at 30 minutes, and 0.100 at 60 minutes. For MP (simulated data), the corresponding values are 0.017, 0.005, 0.0036, and 0.002. Therefore, MP appears as a one-compartment model, whereas RAC represents a two-compartment model.

For the image data (RAC), the ratios are 0.25, 0.10, and 0.085 for 5, 20 and 30 minutes and 0.078, 0.029, and 0.0208 for the same times for MP image data. The RAC image data still requires a 2-compartment model, but MP image data is not as clear as the simulation data and small differences between DVs are seen for the 1- and 2-compartment models. This introduces some ambiguity into what should be taken as the “true” value for the DV in this example. Also, attempts to estimate the BP directly are likely to be subject to large errors when the integrated response function ratios are small.

Because the RAC data requires a two-compartment model, the authors must verify that the method used here can in fact describe the two-compartment model. This is illustrated in Fig. 8. The solid line is the two-compartment solution to Eq. 7 and ⋄ represents the two part fit using Feng's one-compartment model. The two part fit differs slightly from the true model at early time points but is identical with the model at later times, which determine the DV. An “exact” fit to the model (indicated by ♦ in Fig. 8) can be obtained by including a third term in the smoothing operation of the first time interval so that Eq. 2a has an additional term on the righthand side consisting of Q exp(- ${\hat{K}}_{2}$ t), where Q is a constant (in addition to K₁ and k₂) to be determined from the linear solution of the matrix equations. The solution now requires inversion of a 3 × 3 matrix. In practice with noisy data the inclusion of an extra term provided results no different from the procedure using only two terms, so the data reported here is based on the two term solution as described. Thus, the “smoothing” method proposed here adequately represents the two-compartment model.

FIG. 8.

Data generated using the two part fit (without added noise) for 11C raclopride (RAC). ** indicates the fit using 3 parameters and ** indicates the fit using 2 parameters. The solid line is the original data. Distribution volumes were the same for both (1.91).

In many instances, data from RAC studies can be described by a one-compartment model, in which case the Feng method for a one-compartment model would be adequate. There are many cases in which little difference exists between the fits for a one- or two-compartment model (Koeppe et al., 1991; Carson et al., 1997). There are also instances in which the “nonspecific” reference regions are better described by a 2-compartment model. This has been observed for some studies with the radioligands ¹¹C raclopride (for example, Logan, 2000) and ¹⁸F spiperone (Logan et al., 1987). Also, Abi-Dargham et al. (2000) observed that a two-compartment model gave a somewhat better fit to cerebellar data for the D₁ ligand, [¹¹C]NNC 112. In cases such as these this method would be applicable without the necessity of determining which model to apply.

Although values close to the original DV or DVR were recovered in the mean, there are large variances at the higher noise levels. Certainly some data sets have DVs that differ considerably from the mean (for example, Figs. 2A and 3A). Some difficulties are encountered with all methods in the presence of significant noise. This smoothing appears to do a good job in the intermediate range. Furthermore, the failure rate increases at high noise levels. When the data differ from the “true” signal, the noise even when random can introduce an apparent bias because of the finite number of measurements. By combining several data sets (as in an ROI) a value closer to the true value is obtained. In the case of the graphical analysis, the combining must be performed before the analysis to eliminate the noise and the bias. Also, the bias can be reduced by combining time frames.

The advantage of the linear methods, the graphical methods or the GLLS method of Feng, over that of the nonlinear least squares method is that they provide a considerable savings in computing time, which is especially important in constructing images because of the large numbers of pixels and multiple slices. The time required to generate solutions in the nonlinear approach is dependent upon how close the initial guess is to the best solution and how noisy the data is. The method proposed here is not as rapid as the graphical approach, but because it converges in one or two iterations it provides a considerable savings in time over NLS methods. The reference region method of Lammertsma and Hume (1996) and Gunn et al. (1997), implemented in the linear form of Feng, is also a rapid method that works well when data conform to a one-compartment model, as in the example MP. However, for the two-compartment example, RAC, it overestimates the DVR somewhat at low noise and the overestimate increases at higher noise levels with the CV becoming very large. Similar problems have been observed with the simplified reference tissue model when either the reference region or the receptor containing region cannot be described by a one-tissue compartment model, although use of the more general reference tissue model (Lammertsma et al., 1996) may not have these problems. In any case, the two part smoothing using the same simplified method followed by graphical analysis gives the original DVR. There is some susceptibility to noise at high noise levels resulting in larger CVs.

This method should be able to provide unbiased estimates of the DV (or DVR) for ROI data or images with intermediate noise levels. At high noise levels particularly with ¹¹C tracers, the DV or DVR will be overestimated with a large CV. If a 1-tissue compartment model is appropriate, then the Feng analysis with a measured input function appears to work well even at higher noise levels. The DVR calculations using the smoothing without k₂ provides estimates with the lowest CVs but can underestimate the DVR for ligands with slow efflux such as MP. This can be improved by extending the analysis times for the graphical analysis. For RAC, however, the DVR estimate without k₂ is good even at the highest noise levels. The smoothing method does work well for the two examples of image data presented here. More extensive testing on images with different tracers and different DVs will be performed to establish limits of applicability. The difficulties of parameter estimation in the presence of significant noise are illustrated in Figs. 2A and 3A in which the smoothed data appears to provide a good description of the noisy data but does not provide accurate estimations of the DV.

The noise level in PET data is most frequently a problem for short-lived isotopes (such as ¹¹C) at longer time periods because of radioactive decay. In this study, we investigated the feasibility of a combination of the GLLS and the graphical methods to provide a bias-free estimate of the DV. However, because these methods have significant CVs at high noise levels, it also would be desirable to reduce the noise through improved counting statistics for the generation of parametric images.

References

Abi-Dargham

Martinez

Mawlawi

Simpson

Hwang

Slifstein

Anjilvel

Pidcock

Guo

Lombardo

Mann

Van Heertum

Foged

Halldin

Laruelle

(2000) Measurement of striatal and extrastriatal dopamine D1 receptor binding potential with [¹¹C]NNC 112 in humans: validation and reproducibility. J Cereb Blood Flow Metab 20:225–243

Blomqvist

(1984) On the construction of functional maps in positron emission tomography. J Cereb Blood Flow Metab 4:629–632

Carson

Breier

de Bartolomeis

Saunders

Schmall

Der

Pickar

Eckelman

(1997) Quantification of amphetamine-induced changes in [11C]raclopride binding with continuous infusion. J Cereb Blood Flow Metab 17:437–447

Carson

Kiesewetter

Jagoda

Der

Herscovitch

Eckelmen

(1998) Muscarine cholinergic receptor measurements with [18F] FP-TZTP: control and competition studies. J Cereb Blood Flow Metab 18:1130–1142

Evans

(1987) A double integral form of the three-compartment model, four rate-constant model for faster generation of parameter maps [abstract]. J Cereb Blood Flow Metab 17 (suppl 1):S453

Ding

Y-S

Fowler

Volkow

Dewey

Wang

G-J

Logan

Gatley

S J

Pappas

(1997) Chiral drugs: comparison of the pharmacokinetics of [C-11]d-threo and l-threo-methylphenidate in the human and baboon brain. Psychopharmacology 131:71–78

Feng

Wang

Huang

S-C

(1993) A study on statistically reliable and computationally efficient algorithms for the measurement of local cerebral blood flow with positron emission tomography. IEEE Trans Med Imag 12:182–188

Feng

Huang

S-C

Wang

(1996) An unbiased parametric imaging algorithm for nonuniformly sampled biomedical system parameter estimation. IEEE Trans Med Imag 15:512–518

Gunn

Lammertsma

Hume

Cunningham

(1997) Parametric imaging of ligand-receptor binding in PET using a simplified reference region model. Neuroimage 6:279–287

10.

Hsu

Alpert

Christian

Bonab

Morris

Fishman

(1997) Noise properties of a graphical assay of receptor binding. J Nucl Med 38:204P

11.

Koeppe

Holthoff

Frey

Kilbourn

Kuhl

(1991) Compartmental analysis of [11C]flumazenil kinetics for the estimation of ligand transport rate and receptor distribution using positron emission tomography. J Cereb Blood Flow Metab 11:735–744

12.

Lammertsma

Bench

Hume

Osman

Gunn

Brooks

Frackowiak

(1996) Comparison of methods for analysis of clinical [¹¹C]raclopride studies. J Cereb Blood Flow Metab 16:42–52

13.

Lammertsma

Hume

(1996) Simplified reference tissue model for PET receptor studies. Neuroimage 4:153–158

14.

Logan

Wolf

Shiue

Fowler

(1987) Kinetic modeling of receptor-ligand binding applied to positron emission tomographic studies with neuroleptic tracers. J Neurochem 48:73–83

15.

Logan

Fowler

Volkow

Wolf

Dewey

Schlyer

Macgregor

Hitzmann

Bendriem

Gatley

Christman

(1990) Graphical analysis of reversible radioligand binding from time-activity measurements applied to [N-11C-methyl-(-)-cocaine PET studies in human subjects. J Cereb Blood Flow Metab 10:740–747

16.

Logan

Fowler

Volkow

Wang

G-J

Ding

Alexoff

(1996) Distribution volume ratios without blood sampling from graphical analysis of PET data. J Cereb Blood Flow Metab 16:834–840

17.

Logan

(2000) Graphical analysis of PET data applied to reversible and irreversible tracers. Nucl Med Biol 27:661–670

18.

Press

Flannery

Teukolsky

Vettering

(1988) Numerical recipes in C: the art of scientific computing. Cambridge, England: Cambridge University Press

19.

Slifstein

Laruelle

(1999) Statistical bias in the Logan graphical method. J Nucl Med 40:35P

20.

Volkow

Fowler

Wang

G-J

Dewey

Schlyer

MacGregor

Logan

Alexoff

Shea

Hitzemann

Angrist

Wolf

(1993) Reproducibility of repeated measures of carbon-11-raclopride binding in the human brain. J Nucl Med 34:609–613

21.

Volkow

Ding

Fowler

Wang

G-J

Logan

Gatley

Schlyer

Pappas

(1995) A new PET ligand for the dopamine transporter: studies in the human brain. J Nucl Med 36:2162–2168

22.

Wang

G-J

Volkow

Fowler

Franceshi

Logam

Pappas

Wong

Netusil

(2000) PET studies of the effects of aerobic exercise (running) on human striatal dopamine release. J Nucl Med 41:1352–1356