2.1.

Model Representation of the AIF

A physiologically based pharmacokinetic Markov chain model^9,10 is used to describe tracer atom movement with the systematic circulation (see Fig. 1). There are eight different states, each representing different parts of the body: (1) RV, (2) LV (arterial blood), (3) lung vasculature, (4) lung extravascular space, (5) body vasculature, (6) body extravascular space, (7) loss blood, and (8) venous blood. The random movement of a tracer atom between model states on each time step (heart beat) is governed by a transition matrix. The Markov chain transition matrix is

Fig. 1

(a) Whole-body circulation model for tracer atom dynamics within the body. The model has eight states detailed at the beginning of Sec. 2. States 1 to 6 are indicated with numbers in schematic, the other two states are loss output (lower-left, state 7) and venous blood (blue, state 8). The parameters in the model, $\mathbf{\theta }$ in text, determine the various transition probabilities associated with the black arrows. (b) Illustrative time-activity curve profiles for regions—associated states and numbers indicated—are displayed for ${\mathrm{H}}_2\mathrm{O}$ (${}^{15}\mathrm{O}$-labeled ${\mathrm{H}}_2\mathrm{O}$) and FDG (${}^{18}\mathrm{F}$-labeled FDG). The injection profile durations are 5 s (${\mathrm{H}}_2\mathrm{O}$) and 2 min (FDG). In the plots, the RV profile is normalized to have a maximum of unity with other profiles scaled relative to it.

It is also required that $r=r_v^L+r_v^B$ and $d=1-r_e^L+1-r_e^B-\ell$ be nonnegative. The input to the Markov model is represented by the convolution of a specified injection profile (typically a square wave of known duration) and a gamma density with shape and scale $(\alpha ,\beta )$—the choice of the gamma density, here, is based on empirical performance. In addition, there is a nonnegative scale ($A$) and a delay parameter ($\mathrm{\Delta }$) that adapts to the timing of the injection relative to the measured PET data. Using a combination of logarithmic and logistic transforms, the model is parameterized by a vector $\mathbf{\theta }=({\theta }_{\mathrm{\Delta }},{\theta }_A,{\theta }_{\alpha },{\theta }_{\beta },{\theta }_c,{\theta }_r,{\theta }_{or},{\theta }_d,{\theta }_{od},{\theta }_e,{\theta }_{e^{\prime }},{\theta }_s)$. The elements of $\mathbf{\theta }$ are mapped to the unknowns according to

2.2.

Estimation Method

Based on the above, the estimation of the AIF is replaced by a problem of estimating $\mathbf{\theta }$. Information from the available image data as well as a-priori past experience with the tracer is both used in a combined Bayesian-like estimation strategy. It is important to appreciate that while the circulation model references the whole body, the extraction procedure only requires data from the field of view of the scan—not the whole body. The imaging data are reduced by segmentation/clustering to obtain a statistically representative collection of time courses. The exact number of time courses is selected adaptively based on the characteristics of the image data. Details of the segmentation techniques are given in Refs. 11 and 12—a combination of split-and-merge segmentation and $k$-means clustering is used. Note the present AIF extraction, including the segmentation analysis, is applied separately to the ${\mathrm{H}}_2\mathrm{O}$ and FDG studies. A combined approach could potentially enhance the performance. Each of the segmentation-derived time courses, denoted as $z_j(t)$ with $j$ indicating the segment number and $t$ the time frame, is modeled as a positive linear combination of the eight profiles generated by the PKPB model as well as a cumulative AIF profile, i.e., a Patlak term, representing tissue retention. For a given $\mathbf{\theta }$ value, the model-derived prediction of the segment time course is denoted as ${\widehat{z}}_j(t|\mathbf{\theta })$. Note that ${\widehat{z}}_j(t|\mathbf{\theta })$ is the prediction of $z_j(t)$ based on the best positive linear combination of the $\mathbf{\theta }$-based PKPB model profiles. By fitting the PKPB model to directly sample arterially time courses in historical studies, i.e., prior experience, an a-priori characterization of the population behavior of the parameter vector $\mathbf{\theta }$ can be created. Using $p$ to signify population, we compute the population mean (${\mathbf{\mu }}_p$) and covariance matrix (${\mathrm{\Sigma }}_p$) of ${\mathbf{\theta }}^{*}$. Here, ${\mathbf{\theta }}^{*}$ is defined as the part of $\mathbf{\theta }$ with scale (${\theta }_A$) and delay (${\theta }_{\mathrm{\Delta }}$) elements removed. Our estimation of $\mathbf{\theta }$ is based on an objective function that combines the information from the imaging data at hand and the historical experience relevant to the tracer. The optimal parameter choice, $\widehat{\mathbf{\theta }}$, is chosen to minimize the objective function

This can be viewed as a Bayesian technique.⁶ $\lambda$ is a positive scale factor, which is adjusted to ensure that the estimate remains consistent with a-priori experience. Specifically, the value of $\lambda$ is adjusted to ensure that the Mahalanobis distance of the nonscale and delay parts of the estimate (${\widehat{\mathbf{\theta }}}^{*}$) from the a-priori mean (${\mathbf{\mu }}_p$), i.e., ${({\widehat{\mathbf{\theta }}}_{*}-{\mathbf{\mu }}_p)}^{\prime }{\mathrm{\Sigma }}_p^{-1}({\widehat{\mathbf{\theta }}}_{*}-{\mathbf{\mu }}_p)$, does not exceed the 90th percentile of corresponding distance values observed in historical data. The weights [$w_j(t)$] used reflect the reliability of the PET data—this is based on the approximate Poisson nature of the measurement.¹³ We also adjust the weighting to place more emphasis on the early part of time course (within the first 4 min for ${\mathrm{H}}_2\mathrm{O}$ and 15 min for FDG) where vascular signals are more dominant. The entire approach is implemented with substantial reliance on open-source software, specifically AMIDE¹⁴ and R.¹⁵

2.3.

Scaling

The above estimation process does not reliably identify the absolute scale of the AIF. A separate process is used to scale. Since arterial and venous blood profiles ultimately equilibrate,⁸ it is possible to scale our AIF curve using either arterial or venous sampled data. As long as the timing of these measurements is known relative to the injection time of the tracer, there is no need to restrict to later sampling time points. In particular, the sampling times do not have to be taken at a time when the arterial and venous activities can be considered to have equilibrated. This flexibility allows the blood sampling to occur when the activity levels are high and the relative error in the blood measurements is lower. In the context of blood flow-metabolism mismatch studies, a scaling based on data recovered from measurements made in LV-ROIs or possibly even the descending aorta is certainly feasible, including at early- or late-time points where spillover from adjacent structure does not pose a challenge (early) or arterial–venous equilibrium minimizes the impact of spillover between arterial and venous structures. In the LABCa studies reported below, AIFs were scaled using late-time blood time-course data recovered from the LV-ROI, and this is to ensure better agreement with the LV-ROI-extracted AIFs.

3.1.

Clinical Study and Data

A dataset from a study of PET imaging in LABCa patients is used. The study is fully described in a number of reports, see, for example, Dunnwald et al.² and Mankoff et al.¹⁶ Briefly, this was an institutional review board-approved prospective study at the University of Washington Medical Center. Patients with histologically confirmed breast carcinoma and scheduled for neoadjuvant chemotherapy prior to surgical resection were eligible. Selected patients were imaged at baseline before chemotherapy and also at mid- and posttherapy. Patients underwent definitive breast and axillary surgery posttherapy where any remaining tumor was resected. Our analysis focuses on the baseline data from the set of 53 patients. As previously reported, the patient data include standard clinical diagnostics, including the number of positive axillary nodes and the estrogen receptor (ER) and progesterone receptor (PR) status.¹⁶ Outcome data for treatment for each patient were recorded as the presence or absence of residual tumor at posttherapy surgery [no pathologic complete response (pCR) versus pCR], disease-free survival, and overall survival (OS). Survival and disease progression were assessed¹⁶ from the time of diagnosis. The date of the last follow-up is May 4, 2009, and patients who had not died or progressed at that point were considered censored in the analysis.

3.2.

PET Imaging and Data Processing

Details of PET radiotracer production and dynamic tumor imaging protocols are given in the earlier reports.^2,16 For ${\mathrm{H}}_2\mathrm{O}$, 725 to 1902 MBq of [${}^{15}\mathrm{O}$]-labeled ${\mathrm{H}}_2\mathrm{O}$ in a 1- to 4-mL volume was injected as a bolus; with FDG, 218 to 396 MBq of [${}^{18}\mathrm{F}$]-labeled FDG in 7- to 10-mL volume was injected over 2 min with a constant infusion pump. The ${\mathrm{H}}_2\mathrm{O}$ and FDG imaging studies were conducted in the same session on a 35-plane GE-Advance scanner, starting with the ${\mathrm{H}}_2\mathrm{O}$ study. The axial extent of the scanner field of view is roughly 20 cm. Dynamic ${\mathrm{H}}_2\mathrm{O}$ imaging was carried out over 8.75 min according to the following acquisition sequence (number of frames and their durations given): 1 (1 min) preinjection, 15 (2 s), 15 (5 s), 12 (10 s), 8 (15 s), and 6 (20 s). The FDG injection and dynamic imaging followed ${\mathrm{H}}_2\mathrm{O}$. The FDG acquisitions were conducted over 61 min according to the sequence: 1 (1 min) preinjection, 4 (20 s), 4 (40 s), 4 (1 min), 4 (3 min), and 8 (5 min). All data were acquired in two-dimensional mode using a filtered-backprojection reconstruction algorithm. Circular ROIs of $\sim 1.5\text{-}\mathrm{cm}$ diameter were drawn over tumor and contra-lateral normal breast, a similarly sized ROI was placed over a portion of the LV. These ROIs provided tissue and blood time-course (LV-ROI) data for analysis. The SAIF extraction method, reviewed above, was applied to the data to produce a second AIF. AIFs were scaled using late-time LV blood time-course data. A voxel-level analysis of the dynamic data was also applied¹⁷ to obtain a mapping of the local blood volume for ${\mathrm{H}}_2\mathrm{O}$ and FDG. Tissue time-course data were analyzed using the standard 1- and 2-compartmental models of Kety for ${\mathrm{H}}_2\mathrm{O}$¹⁸ and Huang–Sokoloff for FDG.¹⁹ All numerical computations, including segmentation, SAIF extraction, and compartmental model fitting by nonlinear least squares were carried out using the open-source R Statistical Computing Software.¹⁵ The AIFs produced by the LV-ROI and SAIF were used as inputs for the compartmental analysis of tissue time-course data. Estimates of blood flow ($K^H$, $\mathrm{mL}/\min /\mathrm{g}$) and the tissue FDG-based glucose metabolic rate (${\mathrm{MR}}_{\mathrm{glc}}^F$, $\mathrm{mg}/\min /100\mathrm{g}$ of tissue) were then constructed. In addition, the blood flow-metabolism mismatch ratio $\frac{{\mathrm{MR}}_{\mathrm{glc}}^F}{K^H}$, measuring the glucose utilization rate per unit blood flow, was also created.

3.3.

Statistical Analysis

Pair-wise plots and correlations between LV-ROI- and SAIF-derived kinetic variables, including the blood flow-metabolism mismatch, were evaluated. The impact of the alternative AIFs on the prediction of patient outcomes, pathologic response, and OS from Dunnwald et al.² was also considered. For survival, a multivariate Cox regression analysis, incorporating standard clinical diagnostics and either LV-ROI- or SAIF-derived blood flow-metabolism mismatch, were compared. Formal comparisons were based on the performance of the prognostic model, likelihood, and concordance values were considered. The statistical Bootstrap methodology²⁰ was used to evaluate differences between model performance criteria, appropriately adjusted for bias. Standard errors for differences between LV-ROI and SAIF prognostic models are also reported. A similar approach was used for pathologic response (CR). A multivariate logistic regression incorporating standard clinical diagnostics and either LV-ROI- or SAIF-derived blood flow-metabolism mismatch was compared. The comparisons were based on the likelihood performance (deviance) of the logistic model. All analyses were conduced in the R statistical package.¹⁵

4.1.

Illustrative AIF Reconstructions

Sample SAIF reconstructions for an ${\mathrm{H}}_2\mathrm{O}$ and FDG study are shown in Fig. 2. Included in this are some of the image-derived time courses produced by the automated data-segmentation procedure during the extraction process. The selected time courses shown originate in regions where arterial (LV), prearterial (RV and lung), and postarterial venous signals dominate. The ability of the SAIF methodology to harvest information from these diverse sources is a key feature of the approach. Comparisons between the LV-ROI and SAIF are shown in Fig. 2. The recovered venous time-course information, a by product of the SAIF technique, is also shown. For ${\mathrm{H}}_2\mathrm{O}$, the SAIF is not as peaked as the LV-ROI data. This may be due to early spillover from the RV, which receives the initial bolus signal at the heart. With FDG, the tracer injection was over a 2-min period. We expect less impact from the RV on the LV data. The FDG SAIF has a higher peak than the LV-ROI-based AIF. The SAIF is scaled by the tail of the LV measurement. With FDG, a case could be made to scale by the earlier LV data. Such an approach could bring the SAIF into better alignment with the initial phase of the LV-ROI AIF while also reducing its tail height.

Fig. 2

Sample SAIF reconstructions for ${\mathrm{H}}_2\mathrm{O}$ and FDG studies in a breast cancer patient: (a) estimated SAIF (red) and venous (blue) times-courses, sampled LV data (points); (b) local blood volume maps, produced by residue analysis,¹⁷ are overlaid on the PET attenuation scan; and (c) segment time-course data (points) with vascular (purple) and extravascular (maroon) components indicated.

4.2.

Kinetic Analysis

Time-course data and corresponding kinetic model analysis obtained using the LV-ROI and SAIF AIFs are shown in Fig. 3. The data are from the same study as shown in Fig. 2. In the case of FDG, the LV-ROI and SAIF AIFs deviate by a maximum of 10%; with ${\mathrm{H}}_2\mathrm{O}$, the maximum deviation is 25%. Figure 3 shows differences between the fits of kinetic models based on the LV-ROI and SAIF inputs. This is most apparent for the ${\mathrm{H}}_2\mathrm{O}$ data, perhaps to be expected given the variability of the time-course data. Figure 3 also presents comparisons between the ${\mathrm{H}}_2\mathrm{O}$-derived blood flow, the FDG-derived glucose metabolic rate, and the blood flow-metabolism mismatch. High correlations are observed for all three parameters. Not surprisingly, the strongest relation ($r=0.99$) is found for glucose metabolic rate. The greater divergence is observed for the ${\mathrm{H}}_2\mathrm{O}$-derived blood flow ($r=0.92$). Again, in light of variation in the time-course data for ${\mathrm{H}}_2\mathrm{O}$, the increased divergence is not surprising. The observed correlation for the blood flow-metabolism mismatch is the weakest ($r=0.89$), which is still highly satisfactory. Further comparisons between kinetic parameters, including computed blood and tissue volumes of distributions, also show strong association between values produced using LV-ROI and SAIF as input functions.

Fig. 3

Comparison of LV-ROI and SAIF AIFs in terms of kinetics analysis for tissue data. Top two rows show ROI data versus time (in minutes) for normal (middle) and tumor (right) regions in one patient. AIFs and model fits of ROI data with LV-ROI (solid-red) and SAIF AIFs (dashed-blue) are shown. Kinetic parameter comparisons for the glucose metabolic rate as determined by FDG, MR(FDG) in $\mathrm{mg}/\min$ per 100 units, blood flow, flow (${\mathrm{H}}_2\mathrm{O}$) in $\mathrm{mL}/\min /\mathrm{g}$ units, and the associated blood flow-metabolism mismatch ratio in all 106 ROIs examined are shown in the bottom row (SAIF values on the $x$-axis)—solid/open dots represent tumor/normal ROIs. Linear correlation coefficients are indicated in boxes.

4.3.

Prognostic Comparisons

Table 1 reports results of multivariate Cox regression models for OS using established clinical diagnostics and the PET-derived blood flow-metabolism mismatch values produced using the LV-ROI and SAIF AIFs. The established clinical variables are the PR and ER receptor status (negativity) of the tumor (PR and ER) and the number of nodes involved (${\text{Nodes}}^{+}$). The differences between the Cox model regression coefficients obtained with the different AIFs are $<20\%$ of the associated estimated standard errors—statistically insignificant. The overall fit of the Cox model in terms of log-likelihood and concordance is also shown in Table 1. The Bootstrap estimates the difference in the Cox model fit are $0.15(\pm 1.90)$ and $0.007(\pm 0.025)$ for the log-likelihood and concordance, respectively. These differences are in the order of 1% and nonsignificant in formal statistical terms—$p$-values are shown in the table. Corresponding results for the analysis of pathologic response using the multivariate logistic regression model are given in Table 2. Again, the differences between the logistic regression model coefficients obtained with the different AIFs are $<40\%$ of the associated estimated standard errors—statistically insignificant. The difference in the deviance (negative logarithm of likelihood) of the models is $-3.63(\pm 4.18)$, showing 10% improvement with the SAIF variable, but as with survival analysis, the difference in deviances is statistically insignificant, $p$-value is 0.34.

Table 1

Multivariate Cox regression analysis of time to death. Based on 53 subjects, 10 of whom were reported dead at the last follow-up. Regression model: log(Risk of death)=PR+Nodes++ER+Mismatch.

Table 2

Multivariate logistic regression analysis of response. Based on 53 subjects, 11 had pathologic responses to chemotherapy. Regression model: log(Risk of death)=PR+Nodes++ER+Mismatch.

A direct comparison between the prognostic risk assessments corresponding to the survival and response models in Tables 1 and 2 is presented in Fig. 4. There is strong correlation between the risks evaluated using the LV-ROI and SAIF AIFs. As a further confirmation of the prognostic equivalence of the two AIFs, we compare Kaplan–Meier survival characteristics of patients stratified into high- and low-risk groups according to whether the model risk, calculated using coefficients in Table 1, is above or below the median risk for the sample. Note this stratification means that half the patients are designated as “high risk” and remaining patients are “low risk.” Their stratification produced using the risk calculation based on either LV-ROI or SAIF inputs are very similar (see Fig. 5). This is consistent with the quantitative analysis of likelihood and concordance reported in Table 1. A comparison between the risk distributions, calculated using coefficients in Table 2, among responders and nonresponders is presented in Fig. 5. The LV-ROI- and SAIF-based risks show similar separation between responders and nonresponders. Again, this is consistent with the quantitative analysis of model deviance reported in Table 2.

Fig. 4

Prognostic comparison. Model-defined risk computed from coefficients in Tables 1 (left) and Table 2 (right) resulting from LV-ROI and SAIF (automatic image extraction) generated blood flow-metabolism mismatch values. Correlation coefficients among risk values are shown in boxes.

Fig. 5

Prognostic comparison: (a) Kaplan–Meier survival plots versus time (in months) for patients stratified into high- (solid) and low- (dash) risk groups according to whether their estimated risk is above or below the median risk for the sample—risk values as shown in the left of Fig. 4. The stratification with LV-ROI mismatch is in red, and SAIF mismatch is in blue. (b) Comparison of risk values (red/blue for LV-ROI/SAIF) from logistic regression analysis of pathologic response. Patients who are found to respond to treatment are coded as “Yes” and nonresponders are coded “No.” The risk values are based on the formula in Table 2, as shown in the right of Fig. 4.