Source-to-Target Automatic Rotating Estimation (STARE) - A publicly-available, blood-free quantification approach for PET tracers with irreversible kinetics: Theoretical framework and validation for [18F]FDG.

INTRODUCTION: Full quantification of positron emission tomography (PET) data requires an input function. This generally means arterial blood sampling, which is invasive, labor-intensive and burdensome. There is no current, standardized method to fully quantify PET radiotracers with irreversible kinetics in the absence of blood data. Here, we present Source-to-Target Automatic Rotating Estimation (STARE), a novel, data-driven approach to quantify the net influx rate (Ki) of irreversible PET radiotracers, that requires only individual-level PET data and no blood data. We validate STARE with human [18F]FDG PET scans and assess its performance using simulations.
METHODS: STARE builds upon a source-to-target tissue model, where the tracer time activity curves (TACs) in multiple "target" regions are expressed at once as a function of a "source" region, based on the two-tissue irreversible compartment model, and separates target region Ki from source Ki by fitting the source-to-target model across all target regions simultaneously. To ensure identifiability, data-driven, subject-specific anchoring is used in the STARE minimization, which takes advantage of the PET signal in a vasculature cluster in the field of view (FOV) that is automatically extracted and partial volume-corrected. To avoid the need for any a priori determination of a single source region, each of the considered regions acts in turn as the source, and a final Ki is estimated in each region by averaging the estimates obtained in each source rotation.
RESULTS: In a large dataset of human [18F]FDG scans (N = 69), STARE Ki estimates were correlated with corresponding arterial blood-based Ki estimates (r = 0.80), with an overall regression slope of 0.88, and were precisely estimated, as assessed by comparing STARE Ki estimates across several runs of the algorithm (coefficient of variation across runs=6.74 ± 2.48%). In simulations, STARE Ki estimates were largely robust to factors that influence the individualized anchoring used within its algorithm.
CONCLUSION: Through simulations and application to [18F]FDG PET data, feasibility is demonstrated for STARE blood-free, data-driven quantification of Ki. Future work will include applying STARE to PET data obtained with a portable PET camera and to other irreversible radiotracers.

Introduction

Positron emission tomography (PET) allows for in vivo quantification of brain metabolism of different molecules and neurotransmitter system components such as receptors, enzymes and ion channels. Full quantification of dynamically-acquired PET data provides estimates of the amount of radiotracer that is specifically bound to a target of interest in the brain via binding potentials (e.g., BPF and BPP) in the case of radiotracers with reversible kinetics (e.g., [11C]UCB-J and [11C]PBR28) (Innis, 2007). For radiotracers with irreversible kinetics (e.g., [18F]FDG and [18F]FDOPA), we can quantify the total radiotracer uptake and metabolism via estimation of Ki, the net influx rate of radiotracer into tissue from the vascular compartment (Innis, 2007). Essential to obtaining these quantitative estimates, whether through kinetic compartment modeling or graphical approaches, is knowing the input function. This is the concentration of radiotracer and its radiometabolites in the blood compartment throughout the scan, which allows us to model the radiotracer concentration in tissue as a function of the radiotracer concentration in the blood (Innis, 2007). The best validated and most widely used source of an input function is arterial blood. While arterial blood sampling via arterial catheterization has been safely applied in numerous research studies, it adds patient burden and cost, and can be labor intensive. Therefore, efforts to develop, validate, and disseminate less-invasive PET quantification techniques, can enhance use of fully quantitative PET in both research and clinical settings. Here, we present Source-to-Target Automatic Rotating Estimation (STARE), a novel approach that performs full PET quantification of Ki for PET radiotracers with irreversible kinetics, using only the individual-level PET brain data in a completely data-driven manner, without requiring collection of blood. We introduce the theory behind this approach, and report its initial implementation and validation in 69 previously acquired human [18F]FDG scans (D. Devanand, 2010; Roccia, 2019) and in [18F]FDG-based simulations.

[18F]FDG, a glucose analog, is the most ubiquitously used PET radiotracer and yields information on glucose metabolism, which in the brain is considered a marker for neural activity (Pacák et al., 1969). Due to the requirements for full quantification of [18F]FDG data to estimate Ki and the corresponding metabolic rate of glucose (CMRglu) (Sokoloff, 1977; M. Phelps, 1979), semi-quantitative metrics, such as the standardized uptake value (SUV), have been frequently employed instead (Zasadny and Wahl, 1993). However, without strict standardization of SUVs (Boellaard, 2009; Vriens et al., 2010; Westerterp, 2007), quantification of Ki and CMRglu can be preferable, especially where highly sensitive PET metrics are required to detect subtle biological differences. Therefore, there has been much work to develop methods that can estimate Ki and CMRglu without relying on an arterial input function (AIF). Common classes of less- or non-invasive quantification techniques include reference region approaches, image-derived input functions (IDIFs), population-based input functions (PBIFs), and simultaneous estimation (SIME) of the input function (Guo et al., 2007; Wong et al., 2002; Wong et al., 2001; Chen, 1998; Ogden et al., 2010; Bohorquez, 2020; Zanotti-Fregonara et al., 2020; S. Takikawa, 1993; Cunningham, 1991; Hume, 1992).

Reference region approaches quantify PET outcome measures with respect to the tracer time activity curve (TAC) in a region assumed to be devoid of specific binding to the target of interest for reversible tracers or devoid of an irreversible trapping component for irreversible tracers (Cunningham, 1991). However, for many tracers, e.g., [18F]FDG, which is taken up by all living tissues, there is no valid reference region available. Progress has been made in developing approaches for reversible tracers that correct for the bias introduced by reference region assumption violations (Gunn, 2011; Turkheimer, 2012; Salinas et al., 2015); however, these approaches still yield relative outcome measures (e.g., BPND). To our knowledge, such new methods that would allow for absolute quantification of Ki without sampling any blood do not yet exist for irreversible tracers.

Other proposed less-invasive methods, i.e., IDIFs, PBIFs and SIME, generally seek to recover a proxy for the AIF that is typically “anchored” or scaled to the individual in question, commonly by using one or more blood samples. Obviously, such an approach does not entirely eliminate the need for blood sampling. IDIFs rely on extracting the radioactivity within vasculature in the PET field of view (FOV) and PBIFs rely on blood data previously acquired with the same tracer in other subjects, but both still currently require individual blood-based anchoring for practical application. SIME of the input function achieves less-invasive quantification by fitting the proper tissue compartment model to multiple brain regions’ TACs simultaneously. This allows the free parameters of the model, that is the parameters requiring estimation, to be estimated simultaneously, under the usual assumption that the AIF is the input function common to all regions (Guo et al., 2007; Wong et al., 2002; Wong et al., 2001; Ogden et al., 2010; Bohorquez, 2020; Riabkov and Di Bella, 2002; Feng et al., 1997; Sari, 2018; Zanderigo, 2018; Maroy et al., 2020). In SIME of the input function, the model free parameters are both those describing the tracer kinetics in the tissue (e.g., for an irreversible tracer, the micro-parameters K1, k2, and k3 for each brain region) and those describing the AIF (e.g., the parameters of the model often used for the AIF, which is the sum of three decreasing exponentials). However, to ensure identifiability of all free parameters at the individual level, SIME of the input function also requires “anchoring” the solution using at least one blood sample acquired during scanning (Ogden et al., 2010).

Although theoretically, a single arterial sample for IDIFs, PBIFs, or SIME could be acquired using an arterial puncture, the procedure can generate a sudden reaction in the subject under scanning, potentially causing head motion and alterations in blood pressure and cerebral blood flow that may impact tracer delivery to and washout from the brain. Thus, it is preferable to avoid collecting even a single arterial sample. For [18F]FDG, where the radioactivity in venous and arterial blood closely approximate each other after 40 min post-injection (M.E. Phelps, 1979; Wakita, 2000), SIME, IDIFs, and PBIFs have been anchored with one or more venous plasma samples acquired late in the scan (Guo et al., 2007; Wong et al., 2001; Chen, 1998; Bartlett, 2019; S. Takikawa, 1993; Zhou, 2012; Chen, 2007; Naganawa, 2005; Zanotti-Fregonara, 2009). This approach still requires placement of a second intravenous catheter, in addition to the line used for radiotracer injection, and measurement of venous blood activity in a well counter, adding complexity to PET acquisition.

There has been some success in developing completely blood-free full PET quantification approaches. These solutions often require combinations of multiple techniques to achieve acceptable performance. One solution anchors a PBIF with IDIF information derived from whole-body PET scanning (Naganawa, 2020). Deep learning has also been leveraged to obtain blood-free quantification of [18F]FDG; however, to our knowledge this has yet to be validated with human scans (Wang et al., 2020). Further, machine learning applied to precompiled electronic health record (EHR) data has been combined with SIME of the input function to quantify [18F]FDG without the use of any blood samples (Roccia, 2019). However, these solutions are thus far situation specific, i.e., with whole body scanning, with large sets of biological variables in the form of EHR data, or with large datasets acquired from many subjects with the same radiotracer for training and validation of machine learning algorithms.

We now propose STARE (Source-to-Target Automatic Rotating Estimation), a new, blood-free, data-driven approach to quantification of PET tracers with irreversible kinetics that relies only on individual-level dynamic PET data. STARE utilizes a source-to-target tissue model, where the tracer radioactivity curve in a “target” region is expressed as a function of a “source” region, to eliminate the dependency of compartmental modeling on arterial blood. This source-to-target tissue model must be adapted to allow us to disentangle the parameters of the target region from those of the source region. We do this by considering multiple target regions at once, as a function of the common source, and fitting the source-to-target model across all target regions and the source simultaneously. This approach allows STARE to separate Ki in target regions from Ki in the source region. Differently from SIME of the input function, STARE does not use data from blood samples to “anchor” the solution to the given individual but instead, uses bootstrapped, PET image-derived measures of concentration in the vasculature present in the FOV. We validate STARE in a large set of human [18F]FDG scans in comparison to AIF-based estimation and using simulations.

Materials and methods

Theoretical framework

We first present the theoretical framework for STARE, then the details of our implementation in this set of human [18F]FDG scans in 2.2 Implementation.

STARE is based on a reformulation of the standard two-tissue irreversible compartment model (2TCirr). The 2TCirr model expresses the concentration of radiotracer in a target region of interest (C(t)) as a function of the concentration of radiotracer in the arterial plasma (C(t)) as follows (Phelps, 1979a; Sokoloff, 1977): where t is time, IRF is the impulse response function for the target region, K1, k2, k3 are the micro-parameters for the target region, and ⊗ denotes convolution. If C(t) were available, fitting the model in Eq. (1) to the TAC in a target region would result in estimates of the micro-parameters K1, k2, and k3, and thus, of K as K = K1k3∕(k2 + k3).

Without acquisition of arterial plasma samples throughout the scan, an estimate of C(t) is not available. In the case of a PET tracer with irreversible kinetics, such as [18F]FDG, Eq. (1) typically holds for TACs from any brain region. Therefore, it is possible to reformulate Eq. (1) so that C(t) in a target region is expressed as a function of the TAC in another region (C(t)), denoted here as the “source” region and thus, obviating the need to know C(t). To do so, we first express both C(t) and C(t) with the 2TCirr model as in Eq. (1):

By applying Laplace transformation, substitution, and subsequent transformation back into the time-domain to Eq. (2), C(t) can be substituted out so that the target region TAC, C(t), is expressed as a function of: (1) its own micro-parameters, (2) the micro-parameters describing the source region, and (3) the source region TAC itself, C(t), as follows (see Appendix for full derivation):

Which, we denote the source-to-target tissue model, where t is time and θ are the free parameters (the 2TCirr micro-parameters K1, k2, and k3 for the target and the source region). As shown in Eq. (3), θ comprises the macro-parameters L, M, υ, and ε, which are combinations of the 2TCirr free micro-parameters (see Appendix for full derivation).

Because the reformulation in Eq. (3) by itself, with a single target and the source region, only allows for quantification of PET outcome measures “relative” to the selected source region, we model multiple target regions at once as a function of the same common source region and simultaneously estimate the parameters for all targets and the source at once. This, together with STARE anchoring, allows for absolute estimation of the target and source micro-parameters, and of Ki for each. This estimation of absolute measures is one of the ways STARE differs from reference region approaches, which also make TAC-based substitutions to eliminate C(t), but only yield outcome measures “relative” to the selected reference region (e.g., BPND).

This strategy is analogous to SIME of the input function. In that context, the parameters describing the unknown AIF are estimated together with the free parameters describing the tracer kinetics in the tissue by fitting all regions’ TACs at once, under the assumption that the AIF parameters are in common to all regions. Similarly, here we model multiple target regions at once, under the assumption that the source region parameters that they are expressed as a function of, are in common to all regions. The parameters for the target regions and the source region are then estimated at once with simultaneous estimation as follows.

Once each target region is expressed as a function of the source region (according to Eq. (3)), the weighted sum of squared residuals is used across N target regions, where the residuals are the distances between each measured target TAC (C(t)) (with T = 1,…,N indicating the different target regions) and the corresponding modeled target TAC (f (t, θ)) at each t time point (m = 1, … , n), as follows:

Where, ω in Eq. (4) indicates a set of known weights for the different PET frame durations (as is standard in PET imaging).

Similar to SIME of the input function, however, minimization of Eq. (4) will not yield unique estimates of the free parameters for both the target regions and source region. This is due to the fact that there exist multiple combinations of such free parameters that yield equivalently good TAC fits. For SIME of the input function, this identifiability problem is solved by “anchoring” the solution to a blood sample acquired from the subject during scanning (Wong et al., 2002; Ogden et al., 2010; Bartlett, 2020).

Analogously, in STARE, we ensure identifiability by effectively “anchoring” the estimation process, not to data from blood samples, but to PET-derived measures of activity in the vasculature in the FOV (as are common with IDIF methods), as described in detail in the Implementation section. To do this, an additional penalty term is added to the weighted sum of squared residuals in Eq. (4):

This penalty term enforces identifiability by constraining the solution, and K, the K in region T, to a subject-specific neighborhood around estimates derived from the signal in the vasculature within the PET FOV. As is described in 2.2.1 Implementation: STARE Anchoring, this is estimated in a data-driven manner based on the individual’s PET data. The PET signal is automatically extracted from vasculature within the PET FOV, as in IDIF techniques. However, the variability in kinetics across voxels within this vasculature region is then bootstrapped to generate a wide range of possible micro-parameter estimates for each brain region for a given participant, from which the K values are derived. This range of micro-parameter estimates is also used to bound all free parameters to be estimated during the minimization of the STARE cost function (Eq. (5)). The parameter λ is introduced in Eq. (5) for generalizability to other datasets to balance the contribution of the model goodness of fit and anchoring terms, because the goodness of fit term is dependent on TAC radioactivity units, whereas the anchoring term is independent of TAC units. . In this initial [18F]FDG dataset, λ was arbitrarily set to 1; however, its value should be optimized for other scanners and/or radiotracers.

Unlike reference region approaches, the only required theoretical assumption of STARE for the source region TAC is that it follows the 2TCirr model, as shown in Eqs. (2) and (3), where the source and target region TACs are both assumed to follow a 2TCirr model. Although, any region whose TAC follows the 2TCirr model could be used as the source region, to ensure that final estimates do not rely on an arbitrary choice of source, we elected to allow each region to act in turn as the source region. Final Ki estimates thus result from averaging the Ki estimates obtained for each source “rotation”. The theoretical framework and implementation of STARE anchoring is described in Fig. 1.

Fig. 1.

Graphical representation showing the theoretical framework and implementation of STARE.

Implementation

STARE was applied to 69 [18F]FDG human brain scans (D. Devanand, 2010; Roccia, 2019; Bartlett, 2019) and to simulated data considering TACs from six regions of interest (ROIs): cerebellum, cingulate cortex, hippocampus, parietal cortex, medial prefrontal cortex, and parahippocampal gyrus. STARE implementation is shown in Fig. 1. Matlab 2016b (The MathWorks, Natick, MA) was used for implementation and all subsequent processing.

STARE anchoring

STARE anchors the estimation process to a unique solution for each individual via the penalty term in Eq. (5) and via data-driven upper and lower bounds that are generated and imposed on the free parameters in the model (each region’s micro-parameters). Here we describe how this anchoring can be fully automated.

To summarize, we first use a two-step k-means clustering approach to automatically extract a vasculature cluster by parsing the dynamic PET data into characteristic regions, such as background, brain tissue with irreversible uptake, and vasculature. This is followed by partial volume correction (PVC) of the final vasculature cluster (Fig. 1). To anchor the STARE solution in the correct subject-specific “neighborhood” of the free parameter space, we capitalize on the variability of signal within the extracted, and corrected for partial volume, vasculature cluster by bootstrapping the voxel TACs within the vasculature cluster and obtaining distributions of possible micro-parameter estimates for each region. More specifically:

K-means Step 1.

To automatically select the optimal number of clusters to be extracted, based solely on an individual’s PET data, rather than a priori assuming that a pre-set number of clusters will optimally partition all scans, k-means clustering runs multiple times, each with a different number of extracted clusters. For [18F]FDG, we used a generous range from 6 to 40 clusters. An optimal vasculature cluster is then automatically selected from each k-means runs by: (1) eliminating all clusters where the maximum value of the average TAC within the cluster corresponds to the end-point of the curve because that indicates an irreversible kinetic, which most likely represents tissue; (2) eliminating all clusters whose average TAC shows negative values because that most likely represents background voxels; and (3) selecting from the remaining clusters, the one whose average TAC shows the highest peak value at the earliest time of peak because that most likely represents PET signal arising from blood vasculature.

K-means Step 2.

K-means clustering then runs again only on those voxels within the vasculature cluster selected during Step 1, with the assumption that this cluster represents a gross estimate of the vasculature within the FOV, which might be corrupted by some spill-in from nearby tissue (especially late in the scan) and spill-out of vasculature signal (especially early in the scan). Because it might be that this initial gross estimate is comprised of signal from arteries, veins, sinuses, and tissue, in our implementation, we elected to extract 4 clusters from the Step 1 vasculature cluster in Step 2. Among the extracted clusters, similar to Step 1, the voxels belonging to clusters whose average TAC has the highest peak value are selected as the final vasculature cluster. Prior studies have shown that with [18F]FDG, as well as with a multitude of other radiotracers and drugs, the signal arising from the arterial vasculature early in the scan is higher than the signal arising from the venous vasculature, as well as any other tissues (Bartlett, 2019; Bartlett, 2020; Wakita et al., 2000; Chiou, 1989). Therefore, the k-means cluster with the highest early-scan peak is likely to be the closest approximation of the true activity in arterial blood.

Partial volume correction.

PVC via Single Target Correction (STC) is then applied to the voxels within the final vasculature cluster from Step 2 (Sari, 2017). STC was previously optimized and validated, where PVC is performed on a voxel-wise basis for a single region (i.e., the final vasculature cluster), accounting for the voxel-wise spill-in and spill-out of radioactivity from the vasculature cluster (Sari, 2017). Given the reconstruction parameters of the [18F]FDG dataset considered here, a point-spread function of 5.9 mm full width at half maximum (FWHM) was used for PVC via STC (D. Devanand, 2010; Sari, 2017). This parameter should be selected and optimized based on the scanner resolution and reconstruction parameters (e.g., post-reconstruction smoothing) for the dataset at hand.

Data-driven extraction of Eq. (5) penalty term and parameter space bounds.

To derive the individualized penalty term in Eq. (5), specifically K for each target region T, and the lower and upper bounds for the model micro-parameters [LB, UB] for each micro-parameter for each target region T, we leveraged the variability of the TACs within the final, partial volume-corrected vasculature cluster. Our approach extends beyond standard IDIF approaches that often simply use the average TAC in a vasculature cluster, because these standard approaches typically require blood-based scaling (Zanotti-Fregonara et al., 2011).

First, we simulate many instances of vasculature signal curves (C(t), with b = 1, … , B and B = 1000 in our implementation), by bootstrapping curves that fall between one standard deviation below and one above the average TAC across voxels in the partial volume-corrected vasculature cluster, according to the following: where t is time, μ(t) is the average TAC in the vasculature cluster, σ(t) is the frame-wise standard deviation of voxels in the vasculature cluster, and X is a uniformly distributed random number in the interval [0,1] at each time t for each bootstrapping iteration b. Each bootstrapped vasculature signal curve is then fit with a 3-decreasing exponential model (F(t)), commonly adopted in PET to describe the post-peak blood input function (holding at μ(t) from time zero to the time of peak, to avoid any early-scan non-physiological uptake patterns) (see Fig. 1 for an example of the generated bootstrapped curves).

Each region’s TAC is then fit with the 2TCirr model with the F(t) curves serving as proxies for the input function, yielding B sets of estimated 2TCirr micro-parameters per region. The probability density estimates of these micro-parameter sets are then obtained via nonparametric estimation of the probability density function (Matlab function “ksdensity”). See Fig. 1 for an example of the hippocampus k2 density estimate. The bounds for each micro-parameter, [LB, UB], are then automatically set using the FWHM of the probability density estimate (g) of each k parameter, such that in each region:

See Fig. 1 for an example of [LB, UB] for hippocampal k2 values.

The B sets of estimated 2TCirr micro-parameters are then also used to compute B corresponding K values, from whose distribution, probability density function is also derived in the same manner. The K value used to penalize Eq. (5), is then derived for each region as .

STARE cost function minimization

For each source rotation, Eq. (5) is minimized in parallel via simulated annealing (Pincus, 1970). Simulated annealing was implemented using the Matlab function “simulannealbnd” with default options for initial temperature (100), reannealing interval (100), function tolerance (1e-6), and maximum iterations (Inf). The simulated annealing initial guesses were set randomly in the range [LB, UB] for each parameter, and the search space was confined to [LB, UB] for each parameter.

Vascular correction

The source-to-target tissue model in Eq. (3) assumes that the TACs of target and source regions are corrected for vascular contribution. However, in the absence of measurements of the radiotracer activity in whole blood or plasma, such vascular correction is not easily achieved. We, therefore, implemented an optional vascular correction procedure within STARE. The partial volume corrected average TAC in the vasculature cluster (μ(t)) can be used to perform vascular correction of the TACs in all regions according to the following: where t is time, C(t) is the measured tissue radioactivity from the PET camera, C(t) is the tissue radioactivity corrected for vascular contribution, and V is a user-modifiable vascular volume fraction. We investigated the effect of neglecting vascular correction (V = 0.00), and including STARE’s implementation of vascular correction (with V = 0.05 (Leenders, 1990)), according to Eq. (8), on STARE performance in quantifying Ki.

AIF-based quantification of the available [18F]FDG dataset

An available set of 69 previously acquired and published [18F]FDG scans (D. Devanand, 2010; Roccia, 2019; Bartlett, 2019) was considered, which included participants with mild cognitive impairment, mild Alzheimer’s disease and healthy controls. Per the data sharing agreement, data could be made available by request to Drs. J John Mann/Davangere P Devanand. As previously described, written informed consent was obtained from all participants (D. Devanand, 2010; Roccia, 2019; Bartlett, 2019) and the study was approved by the Institutional Review Board of the New York State Psychiatric Institute and Columbia University. Acquisition details are previously described (D. Devanand, 2010). Arterial plasma was sampled throughout the scan via arterial catheterization as previously described (Devanand, 2010a). To generate a “gold-standard” AIF, the measured values of tracer total radioactivity in arterial plasma were interpolated from time 0 to the time of the plasma peak. After the peak, a sum of three decreasing exponentials was used to fit the radioactivity data via non-linear least squares. The AIF then was used as the input function to the 2TCirr compartmental model to fit each considered TAC, estimate the model micro-parameters, and calculate corresponding Ki estimates. These blood-based 2TCirr Ki estimates were considered to be the “gold-standard” comparison for STARE-based estimates of Ki, although we acknowledge that AIF data, and thus AIF-based Ki estimates, may be prone to noise and error (Graham, 1997). For comparison with STARE, AIF-based 2TCirr was run either while neglecting vascular correction (V = 0.00) or with standard vascular correction, using the plasma AIF (because the radioactivity curve in whole blood was not measured in this previously acquired dataset) with V = 0.05 (Leenders, 1990), according to Eq. (8).

Assessing STARE performance relative to AIF-based quantification

STARE accuracy

Agreement between STARE Ki and AIF-based 2TCirr Ki estimates was assessed using linear regression and Pearson correlation. Regressions and correlations were computed for: (1) all regions and scans together, (2) scan-by-scan, across all regions, and (3) region-by-region, across all scans. Signed percent differences were also computed as: (STARE Ki − 2TCirr Ki)/2TCirr Ki∗100. To assess the possibility of regional dependence in STARE performance, a linear mixed effects model was fit with outcome = the natural logarithm of Ki and fixed effects = region and quantification method (STARE vs. AIF-based 2TCirr). Participant was modeled as a random effect.

STARE precision

Because STARE contains non-deterministic algorithms (i.e., k-means clustering and simulated annealing), the stability in estimating Ki was tested. STARE was run 10 times per scan for a random subset of the N = 69 [18F]FDG scans. Stability of STARE across runs and regions was assessed with the coefficient of variation of Ki estimates (COV = standard deviation / mean).

STARE performance across diagnostic groups

Given that the [18F]FDG dataset includes a transdiagnostic sample, we assessed whether STARE’s performance varies by disease state (participants with mild cognitive impairment (MCI), mild Alzheimer’s disease (AD) and healthy controls (HC)). A linear mixed model was fit as in 2.4.1 STARE accuracy, but with diagnostic group added as a fixed factor (in addition to region and quantification method). The two-way interaction of diagnosis by method to test if Ki varied by quantification method on a diagnosis-specific basis was examined. Statistics were performed in R version 4.0.3 (Bates et al., 2014; Team, 2013).

Simulations

A set of simulations was designed to examine the sensitivity of STARE to variations in its anchoring procedure used to determine K and the upper and lower bounds [LB, UB] of the model free parameters. One set of simulations, (A), investigated the effect of variations in μ(t) (the average TAC across the voxels in the final vasculature cluster) by manipulating its area under the curve (AUC), the curvature of its tail, and its overall shape (in this last case, while holding the AUC constant). Another set of simulations, (B), investigated the effect of variations in σ(t) (the frame-wise standard deviation of voxel activities within the final vasculature cluster). Across all simulations, the same representative [18F]FDG scan was used as a starting point. To quantify how much each simulation altered μ(t) or σ(t) relative to the original curves, a difference score was computed as the absolute summed difference across time-points between the original μ(t) or σ(t) and the simulated μ(t) or σ(t) curves. This was compared to the percent difference in K between the original STARE run and the STARE run following the given manipulation.

Effects of μvasc(t) on STARE performance

Alter μvasc(t) AUC.

The AUC of μ(t) was manipulated while maintaining the shape of the curve, by solving for ad d in the following equation: where t is time and F is the scaling factor that alters the AUC (with F in the range [0.1, 3]). Ad d was estimated at each F instance, and then added to the representative scan’s μ(t).

Alter μvasc(t) tail shape.

The shape of the tail of μ(t) was manipulated by altering the exponential decay of the tail. To do this, the extracted μ(t) was fitted to the usual 3 decreasing exponential model, and the exponential term with the smallest decay constant, which models the late-scan kinetics of μ(t), was modulated using a multiplicative scaling factor to increase or decrease the rate of decay. The same F values as in Eq. (9) were used. The corresponding AUC of the simulated curve was allowed to change with the varying F values.

Alter μvasc(t) overall shape.

The shape of the μ(t) curve was then manipulated while holding the AUC constant, by using the same μ(t) fit as in 2.5.1.2 above. Sets of three decay constants, one per each decreasing exponential, were randomly generated in each simulation and combined with the original corresponding initial magnitude values from 2.5.1.2 above. This new simulated curve was then divided by a factor, diυ, that was estimated in each iteration using a similar procedure to Eq. (9), in order to hold the final simulated curve’s AUC constant at the representative scan’s original value.

Effects of σvasc(t) on STARE performance

The standard deviation of the tracer radioactivity in the voxels within the final vasculature cluster, σ(t), was then manipulated via additive scaling using the same approach as in (Eq. (9)). In this case, however, σ(t) substituted for μ(t) in Eq. (9) and the estimated ad d at each F was added to the representative scan’s σ(t).

Results

STARE accuracy

Blood-free STARE Ki estimates were highly correlated with AIFbased Ki estimates (regression slope (b)=0.88, y-intercept=0.004, Pearson’s r = 0.80, Fig. 2(A), Table 1)). Although the regression slope was less than one, the intercept was greater than zero and STARE Ki estimates were on average modestly overestimated relative to AIF-based estimates (signed percent difference: 5.07% ± 18.14%; Table 1). STARE Ki estimates were 0.00091 ± 0.0041 greater than AIF-based Ki estimates (Fig. 2(B)).

Fig. 2.

Blood-free STARE Ki vs. Ki obtained via arterial blood-based two-tissue irreversible (2TCirr) compartmental modeling in 69 human [18F]FDG scans. In both (A) & (B): Each color corresponds to a single scan, where Ki is quantified for six regions (cerebellum, cingulate cortex, hippocampus, parahippocampal gyrus, parietal cortex, and prefrontal cortex). (A) Scatter plot, with linear regression slope (b), y-intercept, and Pearson’s correlation coefficient (r) reported across all scans and regions. (B) Bland-Altman plot with the mean Ki difference (i.e., overall bias) shown in blue and the 95% confidence interval for the mean difference estimate shown in light gray. Limits of agreement (95% confidence intervals of mean difference) are shown in orange.

Table 1

STARE Ki Performance Compared to 2TCirr Ki.

	Regression Slope (b)	y-intercept	Pearson’s r	Percent Difference (STARE Ki – 2TCirr Ki ) / 2TCirr Ki * 100 (mean ± standard deviation)
Across all regions (Fig. 2 (A))	0.88	0.00	0.80	5.07 ± 18.14
Cerebellum	0.71	0.01	0.66	4.83 ± 18.3%
Cingulate cortex	0.80	0.01	0.68	4.83 ± 18.0%
Hippocampus	0.83	0.00	0.71	5.97 ± 18.7%
Parietal cortex	0.82	0.01	0.72	4.63 ± 17.83%
Medial prefrontal cortex	0.85	0.00	0.72	4.82 ± 18.18%
Parahippocampal gyrus	0.72	0.01	0.66	5.36 ± 18.44%

At a scan-by-scan level, across the six regions considered, the agreement between STARE and AIF-based Ki estimates was assessed and slopes ranged from 0.59 to 1.58, y-intercepts ranged from −0.0036 to 0.0046, and Pearson’s r ranged from 0.95 to 1.00. Although the r values ranged from 0.66 to 0.72 and the slopes ranged from 0.71 to 0.85 for the comparison of individual regions’ STARE Ki estimates vs. AIF-based Ki estimates region-by-region, across all scans (Table 1), there was no significant statistical evidence for a regional dependence in STARE’s estimation of Ki relative to AIF-based Ki estimates (p = 0.999).

In comparison to AIF-based Ki estimates, Ki,vasc values (i.e., the peak of the probability density function of the bootstrapped Ki estimates used for anchoring Eq. (5)) were slightly more biased than final STARE-estimated Ki values, based on regression slopes (see Supplemental Materials). Notably, Ki values obtained by extracting an unscaled, non-bootstrapped, average curve from the clustering algorithm (i.e., the average TAC across all voxels within the extracted vasculature cluster) and then using it as the input function to the 2TCirr model, showed overestimation with respect to AIF-based Ki estimates (slope = 1.12, see Supplemental Materials), and a percent difference with respect to AIF-based Ki estimates that is higher than that of STARE-estimated Ki (percent difference for STARE Ki = 5.07 ± 18.14%; percent difference for Ki from vasculature cluster curve = 15.81 ± 21.37%; Supplemental Materials).

STARE precision

Although STARE includes non-deterministic algorithms, STARE stably estimated Ki across runs (COV = 6.74 ± 2.48%).

STARE performance across diagnostic groups

We tested whether STARE’s estimation of Ki was consistent across diagnostic groups (AD, MCI, and HC). The interaction between diagnosis (AD vs MCI vs HC) and quantification method (STARE vs AIF-based) was non-significant (p = 0.17), indicating that the difference in Ki values estimated with the two methods did not significantly vary across diagnoses.

STARE vascular correction

While the level of correlation between STARE Ki and corresponding AIF-based estimates is not affected by whether the vascular correction strategy is considered (r = 0.78) or not within STARE (r = 0.79), their agreement slightly varies. More specifically, when comparing STARE Ki estimates obtained without the vascular correction strategy to AIF-based Ki estimates obtained with VB = 5%, the slope of the regression line is b = 0.83, (Fig. 3(A)). When the proposed vascular correction strategy is used within STARE and compared again to AIF-based Ki with VB = 5%, the slope of the regression line is b = 0.87 (Fig. 3(B)), suggesting that considering the vascular correction strategy is favorable, at least in this [18F]FDG dataset.

Fig. 3.

Effect of vascular correction on STARE estimates of Ki. In both (A) & (B): Each color corresponds to a single scan, where Ki is quantified for six regions (cerebellum, cingulate cortex, hippocampus, parahippocampal gyrus, parietal cortex, and prefrontal cortex). Linear regression slopes (b), y-intercepts, and Pearson’s correlation coefficients (r) are shown. (A) STARE Ki estimates with VB = 0% (no vascular correction, y-axis) relative to AIF-based 2TCirr Ki estimates with VB = 5% (x-axis) and (B) STARE Ki estimates with VB = 5% (within STARE vascular correction, y-axis) relative to AIF-based 2TCirr Ki estimates with VB = 5% (x-axis).

STARE simulations

Simulation showed that altering the overall amplitude of signal within the vasculature cluster used to generate the STARE anchors (i.e. altering the AUC of μ(t) via additive scaling) had the most substantial impact on STARE Ki estimation (Fig. 4(A)1). To compare with human data, in the 69 [18F]FDG scans, the mean percent difference in AUC between STARE μ(t) and the AIF was −6.10 ± 15.69%. According to our simulations, when the simulated difference in the AUC of μ(t) was 10%, it yielded a −9.97% change in Ki estimates.

Fig. 4.

Simulation results for STARE. The data from a single, representative scan were used to build all simulations in (A) and (B). Simulations A.1 through A.3 alter scaling (area under the curve (AUC) and shape characteristics of the signal in the vasculature cluster (i.e. the mean of voxel radioactivities in the vasculature cluster, μvasc(t))), whereas Simulation B alters the standard deviation of the radioactivity of voxels within the vasculature cluster (σvasc(t)). These properties are used in the generation of the STARE anchors. The left columns of (A) and (B) use the same metric on the y-axes to assess the difference between the original and simulated μvasc(t) or σvasc(t) curves, which sums the difference between the curves across all time-points. The x-axes of the left columns of (A) and (B) are the percent difference in STARE-estimated Ki between that simulation iteration and the original result. If necessary, a center column is shown as a zoomed-in inset of the left column. The right columns of (A) and (B) show the simulated curves (with the original curves shown in black)). Within each simulation (across rows), the colors (on a yellow-to-blue scale) are consistent from left to right.

STARE was relatively robust to alterations in the kinetics of signal arising from the vasculature cluster (Figs. 4(A)2 and 4(A)3) and the level of voxel-wise variance in the vasculature cluster used to generate the STARE anchors (Fig. 4(B)).

More specifically across all simulations, Ki only changed by greater than 50% with respect to the original, non-simulated STARE K values under the following simulated conditions: (1) μ(t) AUC was increased by at least 200% (corresponding to a y-axis difference score of at least 5.40 in Fig. 4(A)1) or (2) μ(t) AUC was reduced by at least 40% (corresponding to a y-axis difference score of at least 2.16 in Fig. 4(A)1). From Fig. 4(A)1, we can observe that: (1) altering μ(t) AUC yielded the expected inverse effects on STARE Ki estimates, such that, for example, increasing the AUC yielded negative biases in STARE Ki; and (2) reducing μ(t) AUC had a greater impact on STARE Ki than increasing μ(t) AUC, presumably due to instabilities in the 2TCirr modeling in the free parameter bound generation because the input function tail approached 0 microCi/cc in activity and at times became negative.

In the second simulation in Fig. 4(A)2, which manipulated the kinetics of the amplitude of signal arising from the vasculature cluster (allowing the AUC to change with change in shape), the simulated Ki estimates changed by less than 50%, even when the rate of μ(t) tail decay was increased by 300% or decreased by 90% (Fig. 4(A)2). In Fig. 4(A).3, where the AUC of μ(t) was held constant, while the shape was varied, a much smaller impact on STARE Ki was observed, with all percent differences less than 15% (Fig. 4(A)3). Similarly, altering the level of voxel-wise variance within the vasculature cluster had little impact on STARE Ki (Fig. 4(B)).

Discussion

We present the theory and initial validation with human [18F]FDG scans, for a new, publicly-available PET quantification approach – STARE – that quantifies the net influx rate (Ki) of radiotracers with irreversible kinetics (https://github.com/elizabeth-bartlett/STARE). This method, in theory, allows for noninvasive quantification of PET data without requiring blood sampling during the course of a PET scan and operates in an automatic, data-driven, individual-subject way. In a large dataset of human [18F]FDG scans (D. Devanand, 2010), STARE Ki estimates showed modest overestimation with respect to, and strong correlation with, “gold-standard” AIF-based Ki estimates, and were precisely estimated. In simulations, STARE Ki was largely robust to deviations in the subject-specific STARE anchoring.

The goal of STARE is to provide blood-free full quantification of data acquired with PET tracers with irreversible kinetics. While validated here using [18F]FDG, STARE’s theoretical framework is based on a manipulation of the general 2TCirr model. Therefore, theoretically, STARE can be applied to any tracer whose kinetics can be fitted with an irreversible compartment model. For tracers with non-brain-penetrating radiometabolites, we expect the source-to-target modeling portion of STARE to perform well, because that portion completely eliminates the dependence on blood measurements and would solely be relying on activity of the parent tracer in tissue. However, the subject-specific anchoring portion of STARE that relies on leveraging the activity in the vasculature in the FOV, would be corrupted by activity from the radiometabolites, and thus, would require optimization and further validation, or development of a different approach to anchoring, depending on the tracer at hand. Key considerations when validating STARE for another tracer include PET signal heterogeneity and noise levels across the brain, according to the PET camera used to acquire the data. This validation, in comparison to AIF-based quantification, can be achieved with a modest sample size.

This highlights one of the advantages of STARE, which is that the approach does not depend on machine learning, for which large sample sizes are often required for adequate model training. In fields like dynamic PET imaging, where data sharing initiatives are still in their infancy, and PET acquisition with arterial blood sampling is costly and complicated, large datasets meeting the appropriate requirements for machine learning are rare, especially for novel radiotracers. One such machine learning method, noninvasive SIME (nSIME), which was previously validated with the same [18F]FDG dataset considered here, trained a model with 83 different predictors extracted from precompiled EHR data to estimate [18F]FDG Ki in conjunction with simultaneous estimation (Roccia, 2019). nSIME performed comparably with blood-free STARE (r = 0.80 STARE, r = 0.83 nSIME (Roccia, 2019); Bland-Altman plots appear qualitatively comparable across methods). STARE appears to exhibit less bias, despite the fact that STARE only requires an individual participant’s dynamic PET data and does not require other participants’ data for model training.

We note that there is inter-subject variation in STARE performance (range of regression slopes of STARE Ki relative to AIF 2TCirr Ki is 0.59 to 1.58 across subjects). Some inter-subject variability is to be expected, given that STARE is estimating a blood-based outcome measure without any blood data; however, this metric (i.e., regression slope at the individual level) is rarely reported in the literature for comparison. It would be expected that STARE would be more sensitive to PET data quality issues (e.g., substantial patient motion or poor signal-to-noise) than blood-based modeling, because STARE solely relies on the PET data, without the additional information provided by the blood data, as in the case of traditional AIF 2TCirr. Further work is required to determine: (1) how the present level of inter-subject variability impacts PET studies (e.g., test-retest and pre/post pharmacological intervention comparisons), (2) whether the current inter-subject variability is stable across PET scanners, tracers, and populations, and (3) whether the inter-subject variability can be reduced with further optimizations.

STARE also includes an option for vascular correction, where, at the user’s discretion, the source and target TACs can be corrected for any desired level of blood volume fraction. While the level of correlation between STARE Ki and corresponding AIF-based estimates is not affected by whether the vascular correction strategy is considered or not within STARE, their agreement slightly varies, and our results suggest that considering the vascular correction strategy within STARE is favorable, at least in this [18F]FDG dataset. The suitability of STARE’s vascular correction strategy should be examined for each tracer and in additional datasets, scanners, and populations, and future work could explore optimization of the proposed strategy or implementation of alternative strategies. Because using a fixed blood volume fraction (here 5%) may compromise quantification in diseases with potential blood volume variations (e.g., in tumors (Aronen, 2000), metabolic disorders such as hypothyroidism (Pak et al., 2020), or stroke (Derdeyn, 2002)), development of a STARE version that includes an estimable VB for each participant and/or brain region could be beneficial.

STARE performance was comparable across brain regions for both accuracy and precision. This finding highlights another key feature of STARE – the source rotation. Unlike reference region methods, where a priori knowledge is required to determine one appropriate region assumed to be devoid of tracer specifically bound to the target, STARE’s source rotation step does not require any a priori determination of a source region. Our findings show that STARE estimates Ki with equivalent accuracy when we compute final Ki estimates as averages across all rotations. We also provide preliminary evidence that STARE is robust to disease-specific uptake patterns, where we found that STARE performance did not significantly vary across HC, AD, and MCI groups, providing initial validation of the assumption that STARE can operate on any TACs/brain regions with irreversible kinetics.

STARE was applied here with six ROIs that we have considered in previous work with the same dataset (Roccia, 2019); however, any region can potentially be quantified with STARE, with the caveat that, as the number of ROIs increases, the dimensionality of the free parameter space also increases (three free parameters per ROI). Running STARE with different numbers of ROIs, and with different ROIs from the original six, produced comparable STARE Ki performance relative to blood-based 2TCirr estimation. The results of these analyses compared to estimation of Ki with AIF-based 2TCirr are as follows: (1) random subsets of the original six ROIs used in STARE, including using five, four, and three of the original ROIs (regression slope = 0.68 to 0.82 and r = 0.62 to 0.73 across ROIs), (2) using three original with three new ROIs (regression slope = 0.83 to 0.95 and r = 0.68 to 0.75 across ROIs), and (3) using six new ROIs (regression slope = 0.86 to 1.07 and r = 0.72 to 0.88 across ROIs).

In our assessment, our proposed strategy of having the source region rotate among the ROI set is favorable with [18F]FDG in brain tissue as applied, for example, to psychiatric or neurodegenerative disorders, where each region of the brain shows relatively similar kinetics and no source region seems to outperform the other. For other tracers or in other conditions (e.g., tumor uptake), one region may outperform another as a source, or there may be a priori reasons for selecting a specific region as the source, in which case, the rotation step can be avoided. Additionally, a weighting scheme could be optimized that weights Ki estimates from particular source rotations more heavily, if there is a specific rationale for it, rather than using the unweighted average approach proposed here.

This is the first presentation of STARE; therefore, there are limitations to the conclusions we can draw, as well as many potential future directions of this work. In this [18F]FDG dataset, λ, which balances the relative contributions of goodness of fit and anchoring terms in Eq. (5), was set arbitrarily to 1. The magnitude of the anchoring term relative to the goodness of fit term, i.e., the relative influence of the anchoring term, will change depending on TAC radioactivity units and thus, the choice of λ should be optimized and set accordingly. In future studies, we propose to assess the influence of λ and optimize it during application of STARE to new [18F]FDG scans acquired with two different types of scanners (CerePET from Brain Biosciences, Inc, Rockville, MD, USA and Siemens Biograph mCT from Siemens Healthineers, Knoxville, TN, USA) as part of grant R01EB026481 (PI: Francesca Zanderigo, PhD). Future work with different tracers will also help optimize the choice of λ. Also in the current implementation, defining ROIs requires a T1-weighted magnetic resonance imaging (MRI) scan. Future work can examine the use of clustering techniques to identify ROIs based only on the PET data. Further, in this first implementation we did not partial volume correct the source and target regions. Future work should explore the influence of partial volume effects on STARE performance, especially in the context of neurodegenerative disorders. Simulated annealing was selected here as the optimizer. However, simulated annealing is time consuming, taking ~20 min total for optimization across the six parallelized source rotations. A more computationally efficient implementation may allow STARE to be extended to voxel-level estimation and involves developing adaptations that are robust in the face of the higher noise level present at the voxel-level. This extension of STARE’s application is part of our planned future work.

In this initial validation, all scans were obtained on the same PET scanner and thus, validation will be required to assess STARE’s generalizability to other scanners and tracers. In simulations, STARE anchoring was robust to deviations in the signal within the vasculature cluster, but can be hindered by large inaccuracies in the mean signal within the vasculature cluster. For this reason, the two-stage k-means extraction process is critical, because it allows for automatic cluster number selection based on the kinetics and noise present within that particular scan. However, in the PVC step, the user must supply the FWHM of the point spread function for the specific scanner and reconstruction approach. It is essential that this value is selected appropriately, as grossly altering the AUC of the mean signal in the vasculature cluster, could begin to effect STARE performance. As detailed in the STARE manual, the user additionally has the opportunity for quality control and manual intervention of the vasculature cluster generation portion of STARE, where the remainder of STARE can proceed from a different vasculature cluster that is selected.

Although not used for this purpose in STARE, the novel vasculature clustering technique could be used to generate an unscaled IDIF. This approach will require optimization (particularly for the PVC step), and validation with external datasets. In this initial dataset, the Ki estimates obtained by extracting an unscaled, non-bootstrapped average curve from the vasculature cluster, have poorer correspondence to AIF-based Ki estimates than STARE-estimated Ki (Supplemental Materials). Furthermore, Ki,vasc estimates (i.e., the peak of the probability density function of the bootstrapped Ki estimates used for anchoring in Eq. (5)) are slightly more biased than final STARE-estimated Ki values. We do not expect Ki,vasc or an unscaled, non-bootstrapped curve from the vasculature cluster to perform as well in external datasets, where direct effects of scanner FOV, spatial resolution, smoothing during reconstruction, motion, etc., have yet to be investigated. However, the simulations conducted here suggest that the STARE algorithm as a whole will be robust to deviations in the accuracy of the anchoring procedure, which we will evaluate in ongoing studies.

Conclusions

STARE - Source-to-Target Automatic Rotating Estimation – is a novel approach for automated, blood-free quantification of the net influx rate (Ki) of PET radiotracers with irreversible kinetics, that relies solely on the individual’s dynamic PET data. Letting each of the brain regions for which quantification is desired to act, in turn, as a common “source” brain region for all other “target” regions, allows STARE to gain strength by exploiting the information in multiple regions at once with the goal of accurate and precise estimation of Ki across all regions. STARE is “anchored” for subject-specific identifiability using information derived from a novel vasculature cluster extraction and bootstrapping procedure. We validated STARE with a set of human [18F]FDG scans, finding stable estimation of Ki values that were correlated with “gold-standard” AIF-based estimates. STARE performance was also similar across diagnostic groups. In simulations, Ki estimates were largely robust to characteristics of the vasculature cluster of voxels used for STARE anchoring. With more validation, STARE can accelerate use of quantitative PET imaging in the clinic by simplifying acquisition, reducing cost, and affording individualized, non-invasive quantification.