Abstract
Many quantitative imaging protocols that make use of a metabolitecorrected arterial input function require the use of a mathematic model to describe the rate of metabolism of the radioligand. Commonly, parametric models are fit to metabolism data and then the fitted model is used to correct the plasma input function. ^{11}CWAY 100635 is a rapidly metabolized radioligand used extensively in mapping the 5hydroxytryptamine receptor 1A system. Methods: To evaluate the adequacy of fit of 4 metabolite models, we examined data from 92 subjects who received an injection of ^{11}CWAY 100635, were imaged with PET, and underwent measurement of total plasma concentration and metabolites. The performance of these models was assessed according to residual plots, as well as fit and information criteria. Results: The study showed that the choice of model has a substantial effect on the resulting estimates of outcome measures. Conclusion: Among the models considered, the Hill model provides the best fit across all criteria.
The radioligand ^{11}CWAY100635 is used extensively in imaging the 5hydroxytryptamine receptor 1A system (1–5). In such studies, the method for estimation of binding potential (maximum number of binding sites divided by dissociation constant) that involves the fewest assumptions requires determination of the metabolitecorrected arterial input function. Typically, total radioactivity counts in the plasma are corrected for the rate of metabolism of the injected “parent” (unmetabolized radiotracer). In the case of the PET ligand ^{11}CWAY100635, metabolism is rapid, with nearly 80% of the parent compound being metabolized within 5–15 min.
For accurate quantification of the binding potential, we have shown that a minimum of 110 min of emission scanning are required (6). Metabolites cannot be determined reliably at late times because of low levels of parent compound; low counting rates decrease confidence in count estimates; and ^{11}C decays relatively quickly: For these reasons, the metabolite curve is typically fit by a parametric model and extrapolated to late times according to the fitted model.
Models that have been used for this purpose include 1exponential (6), 2exponential with 1 time constant constrained (7), 3exponential (8,9), multiexponential (10,11), and a Hill function (12,13).
We examined 4 candidate models for the rate of metabolism of the parent compound in plasma over time. Model fit was assessed using residual sum of squares, Akaike information criterion (AIC) (14), and residual plots.
MATERIALS AND METHODS
Data Description
We considered data from 92 subjects recruited in an ongoing neuroreceptor imaging study. Imaging and measurement of metabolites were performed as described previously (6). Briefly, after acetonitrile extraction of the plasma, samples were separated into parent and 1 polar metabolite peak by highpressure liquid chromatography on 5 samples. Extraction was measured at every time point and did not vary with time. The rate of metabolism is presented as a ratio of the area under the curve of the unmetabolized compound peak divided by the total unmetabolized and metabolized peak. In other words, the fraction of unmetabolized compound iswhere v_{i} is the radioactivity count, corrected for background and decay, under peak i. The metabolites are eluted from the column first and collected in tubes (fractions) numbered 1 through 3. The parent is eluted later, with no overlap between the parent and metabolite, and collected in fractions 4 and 5. The fractions are then counted in the well counter, and the percentage of unmetabolized compound is calculated as above.
Modeling Metabolite Data
We considered 2 types of models. One approach is to model the proportion of unmetabolized parent ligand. An alternative is to model the concentration of the parent compound directly with a physiologic model of metabolism. At time t, we used C_{p}(t) and C_{met}(t) to denote the concentrations of parent compound and metabolite, respectively. Total concentration is given by . Thus, “proportion models” are fit to describe the fraction ; the kinetic modeling alternative is to fit the measurements directly.
Proportion Models.
Observed fractions taken at corresponding time points are modeled as , where represents the noise. Proportion models are constrained in 3 ways. Specifically, because the compound is injected at time zero, it is assumed that the compound is not metabolized and hence (the first constraint). After that, the compound is continuously metabolized (the second constraint), and its metabolism is irreversible (the third constraint).
The first proportion model considered is the 1exponential model:Eq. 1The restrictions are satisfied provided that and .
Another common model is the 2exponential model:Eq. 2We arbitrarily let represent the larger of the 2 time constants (i.e., ). Restrictions on this model are satisfied if .
The third model for metabolite fraction considered is the Hill model (15):Eq. 3for which we require , , and .
A Kinetic Model.
Knowledge of the kinetics of the metabolism process may improve our ability to model the data. To our knowledge, no published paper using WAY100635 has applied a kinetic model to fit the metabolite function, but it seems important to consider at least 1 kinetic model for comparison's sake. We considered a basic 2compartment kinetic model that allows for 1 metabolite as diagrammed in Figure 1. Though alternative kinetic models could be fit, model complexity is limited by the small number of data points per subject.
The total concentration of radioactivity in the plasma is measured over time through plasma analysis. The concentration of metabolite is measured through metabolite analysis, but cannot be observed directly. In our study, the metabolite data consist of the fraction of unmetabolized compound at specific time points (), and the plasma data consists of time points and corresponding total concentrations of radioactive substances (i.e., ).
The system of differential equations corresponding to Figure 1 in terms of rate parameters isThe solution for is given byEq. 4where ⊗ denotes convolution.
The total concentration in Equation 4 was estimated from a plasma data model. Our plasma model comprised 2 pieces: After the injection, the concentration of the compound in plasma increased linearly for a short time to its peak value; after peaking, the concentration fell off according to a sumofexponentials model (6). The linear rise in the plasma function until the peak was a fairly crude representation of the true plasma concentration curve, but because this interval was small relative to total scanning time (about 1%), it makes little difference in the resulting model estimates. Once the total plasma data are fit to this model, the estimated plasma concentration function is plugged in for estimating .
Weighting.
The data used for the fitting are computed as ratios of radioactive counts in each fraction. In general, the fraction of unmetabolized (parent) compound is computed for a given time aswhere I represents an index set. For our data, indices range from 1 to 5 and I = {4, 5}.
The γcounter (1480 Wizard 3M Automatic γCounter; Wallac) provides both the values and an estimate of their SD. Using these quantities, estimates of the SD of the fractions can be computed using the socalled delta method. If each denotes the SD of the corresponding measurement and the 's are taken to be independent, then the variance of Y is approximated according toEq. 5This is used to determine the weights used in model fitting.
ModelFitting Algorithm.
Standard iterative nonlinear regression methods are used for the model fitting:Eq. 6where β is the vector of the parameters in the models and is computed as according to Equation 5.
RESULTS
Visual Analysis of Model Fits
Figure 2 gives example fits for each model for a specific dataset. For the kinetic model, to equalize the modelfitting procedure (to allow a more direct comparison), we computed the corresponding fraction bywhere comes from Equation 4 and then is fit using Equation 6. Though the fitting is done using the fractional data, results are displayed using metabolite concentration data. The 2exponential model clearly undershoots the later time points, because the 2exponential model does not contain a constant term and is thus constrained to go to zero rapidly, even when the data do not justify this (as in this example). The data were also fit with a 2exponential model with constant (results not shown), but the fitting procedure was rather unstable (many parameters with few data points) and for almost all subjects, one of the exponential terms had a zero coefficient, providing no improvement over the 1exponential model.
Comparisons Based on Fit Diagnostics
We also compared weighted residual sum of squares and AIC for each model across all datasets. AIC (14) attempts a compromise between model fit and parsimony by imposing a penalty for additional model parameters. For normally distributed errors, the criterion is AIC = nlog(weighted residual sum of squares/n) + 2M, where n is the number of data points and M is the number of model parameters.
Figure 3 shows box plots for AIC values corresponding to each candidate model; corresponding box plots for log(weighted residual sum of squares) are similar and not shown. On the basis of these measures, the Hill model fits the data best among the models considered.
Residual Analysis
If the model fits well, at least in large datasets, residuals for each dataset should have an approximate normal distribution with a mean of zero. We had no more than 6 observations in each dataset and thus could not expect residuals to have a distribution close to normality. However, for an appropriate choice of model, we would expect residuals to be centered on zero and to be nearly symmetric.
In a typical regression analysis, many observations are pooled to fit a model and all residuals are plotted together in a single plot. In our present study, data from each subject were fitted separately. In such a situation, we search for patterns among independent residuals by pooling across all subjects separately for each time point.
All residuals from all 92 subjects are displayed simultaneously in Figure 4. On the basis of both absolute size of residuals and distribution of residuals about zero, the Hill model would be preferred.
Effect on Kinetic Outcome Measures
A crucial question about metabolite modeling is whether the choice of model makes a difference for subsequent regionofinterest or voxellevel analysis. Therefore, we compared the resulting biologic parameters binding potential (BP_{1}) and volume of distribution (V_{T}) that resulted from different metabolite models.
For each candidate model, the resulting estimated metabolitecorrected plasma concentration curve was used as the input function for subsequent kinetic modeling of PET data. This procedure was applied to all 87 complete datasets to estimate total equilibrium V_{T} (using a 1tissue kinetic model) for cerebellum, and BP_{1} (using a 2tissue constrained model) (6) for 4 regions: amygdala, hippocampus, cingulate, and parietal lobe.
The effect of the metabolite modeling is summarized in Table 1. To test for significance of the differences, paired t tests and Wilcoxon signedrank tests were applied. Results are given in Table 2. All P values are small, indicating that the choice of metabolite model had a highly significant impact on resulting outcome measures derived from the imaging data.
DISCUSSION
Having an arterial input function for modeling of dynamic PET data can be a great help, provided that the input function is measured and modeled appropriately. Key to the validity of the input function is correctly accounting for the metabolism process of the ligand. Thus, it is of vital importance to ensure that the model used to describe metabolism is adequate.
We considered only 1 kinetic model, though others are possible. Huang et al. (16) presented the kinetic model building procedure and model formulas for the general case, including the case of multiple metabolites. Andree et al. (17) showed that WAY100635 has more than 1 metabolite in humans. However, they also pointed out that only 1 metabolite (^{11}CdesmethylWAY100635) can be used as a radioligand for central 5hydroxytryptamine receptor 1A; other metabolites are more polar and not likely to pass the blood–brain barrier. Osman et al. (18) observed another significant polar metabolite in some experiments—^{11}Ccyclohexanecarboxylic acid—along with some other metabolites that were more polar still. Thus, with the exception of ^{11}CdesmethylWAY100635, no metabolites are expected to bind strongly to 5hydroxytryptamine receptor 1A. Even if the radioactive acid metabolite binds nonspecifically, it is in low levels and cleared rapidly (in about 5 min). On the basis of these considerations, as well as the obvious limitation on the complexity of models that can be fit with 6 data points, we concluded that the kinetic model considered in the data would be an appropriate candidate for these data.
CONCLUSION
Among those considered, the Hill function was the model of choice to describe the metabolism of ^{11}CWAY100635, and the impact of model choice was confirmed by the analysis on BP_{1} and V_{T} data, demonstrating the need to carefully select an appropriate metabolite model in analyzing PET data.
Acknowledgments
This work was supported in part by MH62185 and NARSAD.
Footnotes

COPYRIGHT © 2007 by the Society of Nuclear Medicine, Inc.
References
 Received for publication November 14, 2006.
 Accepted for publication March 16, 2007.