Novel voxelwise residual analysis of [ 11 C]raclopride PET data improves detection of low-amplitude dopamine release

Abstract

Existing methods for voxelwise transient dopamine (DA) release detection rely on explicit kinetic modeling of the [¹¹C]raclopride PET time activity curve, which at the voxel level is typically confounded by noise, leading to poor performance for detection of low-amplitude DA release-induced signals. Here we present a novel data-driven, task-informed method—referred to as Residual Space Detection (RSD)—that transforms PET time activity curves to a residual space where DA release-induced perturbations can be isolated and processed. Using simulations, we demonstrate that this method significantly increases detection performance compared to existing kinetic model-based methods for low-magnitude DA release (simulated +100% peak increase in basal DA concentration). In addition, results from nine healthy controls injected with a single bolus of [¹¹C]raclopride performing a finger tapping motor task are shown as proof-of-concept. The ability to detect relatively low magnitudes of dopamine release in the human brain using a single bolus injection, while achieving higher statistical power than previous methods, may additionally enable more complex analyses of neurotransmitter systems. Moreover, RSD is readily generalizable to multiple tasks performed during a single PET scan, further extending the capabilities of task-based single-bolus protocols.

Introduction

Dopamine release in PET

Neurotransmitter release can be detected in vivo using dynamic (4 D) positron emission tomography (PET) imaging. To date, the most common mechanism for detection relies on competition for binding between endogenous neurotransmitter and a PET ligand; ¹ during release, a localized time-varying increase in endogenous neurotransmitter concentration results in a decrease in ligand binding in that region, ultimately causing a decrease in the PET time activity curve (TAC).

The majority of research thus far has focused on the neurotransmitter dopamine (DA), due to the important role DA plays in healthy brain functioning and to the favorable binding properties of the relatively commonly available PET ligand [¹¹C]raclopride (RAC), a D2/D3 receptor antagonist. RAC DA-detection frameworks mostly confine their search to the D2-rich regions of the striatum where the signal-to-noise ratio (SNR) is adequate for reliable detection. DA is involved in motor function, executive function, and reward,^{2 –4} making it a useful research target for neurodegenerative diseases, neuropsychiatric disorders, and neuropsychological tasks in the healthy and diseased brain.

In its earliest form, this competition mechanism was exploited in a double scan protocol (double RAC), where one scan occurred with the subject in the baseline condition and a second scan occurred with a sustained stimulus evoking DA release. ⁵ The change in binding potential estimate from the baseline to stimulus scan is a proxy for the evaluation of the magnitude of DA release. However, such a protocol has several disadvantages: scan-to-scan variability in the binding potential confounds stimulus-evoked changes; increased commitment and radiation dose is required for study participants; and the release magnitude is assumed to be constant over time for analysis purposes, which may not be true. More advanced double RAC analysis methods leveraged the residual behavior between the baseline and activated scans by use of singular value decomposition (SVD) to accurately determine DA release at varied task timings, ⁶ however this method is practically limited since N + 1 PET scans are required to analyze N tasks. ⁷

Transient release modeling

More recently, several competition models have been suggested to quantify the transient kinetics of neurotransmitter release within a single scan: the subject is scanned in baseline condition during a first part of the scan that includes tracer uptake, then the stimulus is presented later in the scan to observe and quantify any transient responses.^{8 –10} Such frameworks use compartmental models that include time-varying kinetic parameters to model the dynamic competition between the tracer and endogenous neurotransmitter.

Stimulus-evoked endogenous DA changes often occur over second or sub-second intervals. ¹¹ However, the clearance rate of RAC is on the order of minutes, meaning high frequency endogenous DA changes are mapped to low frequency RAC changes. This is similar—conceptually—to the hemodynamic response in functional magnetic resonance imaging (fMRI): neuronal activation occurs on the order of milliseconds, but the measured hemodynamic response via the Blood Oxygen Level Dependent (BOLD) signal evolves on the order of seconds.^12,13

Transient release detection

An initial voxelwise model comparison framework was formalized and optimized to detect and localize transient DA release.^14,15 The framework employed a voxelwise F-test between the multilinear reference tissue method (MRTM) ¹⁶ —incapable of modeling competition—and a method that models the competition, linear parametric neurotransmitter PET (lp-ntPET), ¹⁰

$$\begin{array}{cc}\mathbf{MRTM}:C_T\left(t\right)=R_1{C}_R\left(t\right)+k_2{\int }_0^t{C}_R\left(u\right)\text{d}u& \hspace{1em}\hspace{1em}\hspace{1em}-k_{2a}{\int }_0^t{C}_T\left(u\right)\text{d}u\end{array}$$

(1)

$$\begin{array}{c}\mathbf{lp}-\mathbf{ntPET}:\\ C_T\left(t\right)=R_1{C}_R\left(t\right)+k_2{\int }_0^t{C}_R\left(u\right)\text{d}u\\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}-k_{2a}{\int }_0^t{C}_T\left(u\right)\text{d}u-\gamma {\int }_0^t{C}_T\left(u\right)h\left(u\right)\text{d}u\end{array}$$

(2)

where $C_T\left(t\right)$ is the target region tissue concentration; $C_R\left(t\right)$ is the reference region tissue concentration; $h\left(t\right)$ is the pre-specified shape of the PET response to DA release characterized by the onset timing of release $t_D$ , peak timing, $t_P$ , and sharpness of release, $\alpha$ ,

$$h\left(t\right)=\theta \left(t-t_D\right){\left(\frac{t-t_D}{t_P-t_D}\right)}^{\alpha }\text{exp}\left[-\alpha \left(\frac{t-t_P}{t_P-t_D}\right)\right]$$

(3)

where $\theta$ is the Heaviside function. The coefficients of each term in (1) and (2) are determined through linear regression,^14,15 then computation and thresholding of the voxelwise F-statistic provides binary prediction maps of DA release. HighlY constrained backPRojection (HYPR) ¹⁷ denoising was applied to the data, which has been shown to improve detection sensitivity by indirectly decreasing temporal noise levels. ¹⁵ Applying this methodology to cigarette smoking, spatiotemporal differences in the release profiles of men and women have been detected. ¹⁸ A separate study was able to map out consistent spatiotemporal patterns of release in response to cannabis smoking. ¹⁹

More recently several additional approaches have been developed in an effort to increase the accuracy and reliability of detection: direct reconstruction with lp-ntPET, ²⁰ a machine learning denoising approach using a feed-forward perceptron neural network, ²¹ a machine learning detection method that trains personalized neural networks on simulated lp-ntPET data, ²² a Monte Carlo extension to the F-test approach, ²³ and a Bayesian framework for lp-ntPET parameter estimation. ²⁴ At the core, though, all of these methods rely on the accuracy of a kinetic model—namely lp-ntPET—and a reliable fit of this model to the data.

Disadvantages of kinetic model-based detection

The lp-ntPET model incorporates a family of basis functions that describe physiologically-plausible release timecourses. Thus, the number of degrees of freedom used to parametrize the F-distribution used for significance testing needs to be modified based on the number of basis functions chosen; failing to do so will alter the expected false positive rate (FPR). ²⁵ Additionally, the F-test itself assumes that the data are normally distributed, but image noise in PET is best described by a gamma distribution. ²⁶ Furthermore, detection of release is greatly improved by significant post-processing, ¹⁵ which presumably further violates the normality assumption of the F-test.

While the model comparison approach has produced high detection sensitivity in simulations of relatively high levels of DA release expected in pharmacological challenges, such as cigarette smoking,^14,15 model comparison in simulations of lower DA release levels—which typically are present in motor, cognitive, or neuropsychological tasks—^11,27,28 proved less successful. ²³ The major goals of this work are thus three-fold:

Develop a DA release detection method that does not rely on direct kinetic modeling.

Improve detection sensitivity for lower levels of DA release.

Demonstrate feasibility of the method for a single bolus injection of RAC.

Kinetic models like lp-ntPET were initially developed to accurately quantify release and not specifically to be embedded within a detection framework. ¹⁰ Thus, we hypothesize that decoupling kinetic modeling from detection will increase the reliability and localization of DA release, after which lp-ntPET or a similar model may be used to accurately quantify the release profiles at these locations. If our hypothesis is verified—especially for lower levels of DA release—the ability of PET to detect transient fluctuations in DA levels in response to motor, cognitive, or neuropsychological tasks will enhance studies of the dopaminergic system.

Isolating DA-induced perturbations via residuals analysis

As an initial effort to decouple detection from kinetic modeling, we propose to shift detection to a residual space where DA-induced perturbations to the TAC are isolated. The percentage difference between an observed voxel TAC and a data-driven hypothetical baseline TAC define the voxelwise residuals. Any change to the observed TAC during the task will then appear as structured deviations in the residual space. This mapping to the residual space has conceptual similarities to previous work, ²⁹ and acts to remove the confound of baseline pharmacokinetics and provide a constant signal during baseline. A model comparison via F-test is thus no longer required and simpler detection methods—like those used in fMRI—may be applied. We hypothesize that the isolation of DA-induced perturbations will increase detection sensitivity compared to the model comparison approach by overcoming some of its shortcomings. In principle, this further enables the detection of multiple-task, condition-specific DA release using residual analysis.

Due to typically high levels of noise in voxel-level TACs it is necessary to apply denoising to the data before applying subsequent analyses on voxel residuals. To best characterize the voxel residuals we rely on a combination of novel denoising methods, namely HYPR4D kernelized reconstruction,^30,31 along with denoising using IHYPR4D, ³² both of which extend the HYPR framework to the temporal domain (Figure 1). HYPR4D kernelized reconstruction has been shown to outperform the standard and current state-of-the-art clinical reconstructions in terms of 4 D noise reduction while preserving spatiotemporal patterns within data. IHYPR4D has been shown to provide improvements over HYPR3D in joint noise and bias reduction, and inherently bypasses limitations imposed by HYPR3D; for specific details see literature^{30 –32}.

Figure 1.

Simulated noisy voxel TACs with different reconstruction and denoising approaches. Gaussian noise is added to simulated single-task release ground truth voxel TACs, with noise levels matching real noise measured on the GE SIGNA PET/MR for each reconstruction protocol. Results in this paper use PSF-HYPR4D-K-TOFOSEM+IHYPR4D (right panel).

This paper proceeds as follows. We first outline how to calculate the residuals from acquired data in the context of a DA release experiment. We then propose a detection method in residual space inspired by fMRI ³³ : using a simpler General Linear Model (GLM) rather than kinetic modeling in image space. Collectively, we refer to this process as Residual Space Detection (RSD). Next, we test and optimize RSD on simulated single-task DA release experiments under the conditions of homogeneous and heterogeneous RAC binding; previous simulation studies assume homogenous RAC baseline binding across the striatum, which may not be valid in diseased populations that are of interest for DA studies, such as Parkinson’s disease (see Supplemental Fig. 2). Finally, the optimized RSD method is applied to single-task human scans as a proof-of-concept. Simulation and human results from RSD are compared with results from the traditional lp-ntPET model comparison method.

Theory

Transforming time activity data to a residual space

The first step of this work is to define a hypothetical baseline TAC for each voxel that may be observed in the absence of DA release. DA release in response to non-pharmacological interventions is likely to be localized, ³⁴ and of relatively low amplitude, meaning that a significant fraction of striatal voxel TACs are unlikely to contain DA-induced perturbations. Therefore, on first approximation it is reasonable to assume that the mean striatal TAC appears mostly devoid of DA-induced perturbations and can approximate baseline conditions.

Under the assumption that the mean striatal TAC and measured voxel TACs behave very similarly prior to a DA release-inducing task, a first attempt to define the voxel-level baseline TAC is to use the mean striatal TAC and fit it to a given voxel TAC for timepoints prior to task initiation. To attain a full scan-length voxel-level baseline TAC, the scaling factor from this pre-task fit can be applied to the entire mean striatal TAC, hence extrapolating baseline behavior to later time points. Comparing this extrapolated fit with the measured voxel TAC yields structured residuals in voxels where DA release is present (Figure 3).

Two immediate issues arise when attempting this protocol: (i) TACs of voxels near the boundaries of the striatum have poor fits to the mean striatum TAC due to partial volume effects, and (ii) in the presence of heterogeneous tracer binding within the striatum, defining a single mean striatal TAC as a regressor will result in biased fits and extrapolations for voxels that have higher or lower non-displaceable binding potentials ( $B{P}_{ND})$ than average. Further, some contamination from TACs corresponding to DA release voxels can still be present in spite of the above mentioned assumptions.

The first issue is remedied by extracting the mean TAC from a reference region devoid of specific binding, in the case of RAC using the cerebellum, and adding this as a regressor to the pre-task fit. The second issue may not arise in general as the density of D2 receptors in the healthy striatum is relatively uniform ³⁵ ; this is however not always true, particularly for conditions affecting the striatum such as Parkinson's disease (PD). To resolve this, we propose applying k-means clustering to pre-task voxel TACs normalized by their peak value, with k = 3 to produce three spatial regions representing relatively low, mid, and high $B{P}_{ND}$ within the striatum (various k choices are explored further in the Supplemental Figure 4). Mean TACs for each k-region can then be used as regressors in a non-negative least squares (NNLS) fitting framework for each voxel TAC (Figure 2), assuming a sufficiently large number of voxels are found in each region. To address DA release voxels contaminating these baseline regressors, an iterative approach is introduced below.

Figure 2.

K-means results for noisy simulated data with an anterior-posterior/superior-inferior Bmax gradient. (a) Concentration of a simulated noise-free PET frame (frame 10, t = 25 min) displaying the effect of a simulated B_max gradient in the striatum. (b) Spatial k-regions found by running k-means on pre-task normalized TACs from a noisy simulation. Normalized mean TACs of each k-region are displayed on the plot, demonstrating spatially-varied binding characteristics. These TACs, along with the mean cerebellum TAC, are used to define the voxelwise set of baseline TACs. Voxel and TAC colors are consistent.

Before the commencement of a DA release-inducing task, a striatal voxel TAC, $C_i\left(t\right)$ , is thus approximated as

$$\begin{array}{cc}{C}_i\left(t\right)\approx {\widehat{C}}_i\left(t\right)=\sum _{k=1}^3{w}_{s,k}{\overline{C}}_{s,k}\left(t\right)+w_R{\overline{C}}_R\left(t\right);& t<t_D, w_{s,k},w_R>0\end{array}$$

(4)

i.e., ${\widehat{C}}_i\left(t\right)$ is the hypothetical baseline TAC, determined by performing NNLS fitting to the ith measured voxel TAC prior to the task start time $t_D$ . The NNLS regressors are $\{{\overline{C}}_{s,k}\left(t\right):k=1\text{,}2,3\}$ , the set of mean TACs over the three k-regions, and ${\overline{C}}_R\left(t\right)$ , the mean reference region TAC, here used to account for spill-in of non-target-related activity from surrounding tissue. The corresponding regressor weights are $\{w_{s,k}:k=1\text{,}2,3\}$ and $w_R$ . These weights, obtained from applying NNLS over $t<t_D$ , are applied to their regressors defined over the entire scan duration to extrapolate ${\widehat{C}}_i\left(t\right) \forall t$ . The residuals are then defined as

$$R_i\left(t\right)=100\%\times [1-C_i\left(t\right)/{\widehat{C}}_i\left(t\right)]; \forall t$$

(5)

Note, by this construction, a negative DA-induced perturbation in the voxel TAC post-intervention, $C_i(t>t_D)$ manifests as positively structured residuals (Figure 3).

Figure 3.

Summary of the Residual Space Detection (RSD) and lp-ntPET methods. RSD: Voxelwise residuals are extracted by taking the percent difference between voxel TACs and their predicted baseline TAC (equations (4) and (5)). A predicted residual response is defined (equation (6)) and GLM is computed, assigning $\beta$ weights to each voxel. A single iterative update is performed in which $\beta$ values from the first iteration define new baseline TAC regressors and residual predictors. GLM is recomputed and $\beta$ maps are thresholded to achieve a desired FPR. lp-ntPET: Each voxel TAC is subject to MRTM and lp-ntPET fitting by use of the mean TAC from the cerebellum (equations (1) and (2)). From the set of pre-specified basis functions $h\left(t\right)$ , the function providing the best fit to the voxel TAC is selected, and the F-statistic computed and thresholded. For both methods, in-plane edge voxels are removed from binary maps due to partial volume effects, and a cluster size threshold is applied.

Residual space detection (RSD)

Transforming the PET TACs to residual space allows for analysis and subsequent detection of voxel TACs that deviate from the extrapolated mean baseline behavior as a consequence of task-evoked response. The basis of RSD is to compare residuals to a predictor, $P\left(t\right)$ that approximates an expected residual response to the task (Figure 3). We initially define $P\left(t\right)$ based on an assumed shape of DA release in response to a block task design; however, to avoid potential bias in results introduced by an incorrect initial assumption, we iteratively re-define $P\left(t\right)$ using both the data and the results from the first iteration.

Inspired by fMRI detection approaches—where task events/blocks are convolved with a hemodynamic response function to yield predictors—we first define the predictor as the convolution between an event vector, $E\left(t\right)$ , and a gamma variate response function, $H\left(t\right)$

$$\begin{array}{c}P\left(t\right)=\lambda E\left(t\right)*H\left(t\right)\\ =\lambda \left(\theta \left(t-t_D\right)-\theta \left(t-t_F\right)\right)*\left(\frac{t}{L}\text{exp}\left(-\frac{t}{L}\right)\right);& \underset{\lambda }{\max }\left(P\right)=1\end{array}$$

(6)

where $\theta$ is the Heaviside function, $t_D$ is the delay time indicating the start of the task relative to the scan start time, $t_F$ is the task end time ( $t_D<t_F$ ), $\lambda$ is a normalization constant, and $L$ is a parameter governing the broadness of the gamma variate response function. Empirically, we found that a gamma variate function with $L=10\hspace{0.17em}$ min reasonably matched the shape of residuals derived from noiseless simulations.

A simple means of comparing residuals to predictors is to employ GLM. ³⁶ GLM in the case of RSD is written as

$$R_i\left(t\right)={\beta }_{i}P\left(t\right)+{ϵ}_i\left(t\right)$$

(7)

where $R_i$ (1 × time) are the residuals for voxel $i$ , $P$ (1 × time) is the predictor (initially constructed with equation (6)), ${\beta }_i$ is a weight from regression in units of % signal change (enforced by the residuals definition and normalization condition on $P$ ), and ${ϵ}_i$ contains error in the fit, ideally satisfying typical GLM requirements of $\text{mean}\left({ϵ}_i\right)=0$ , $\text{var}\left({ϵ}_i\right)={\sigma }^2$ .

Finally, to minimize bias and increase accuracy in our results, we shift the method to being more data-driven by iteratively (i) redefining baseline TACs to remove possible contamination from DA release-related signals, and (ii) redefining predictors to reduce possible bias introduced by an incorrectly assumed response function (Figure 3). Baseline TACs, and thus residuals, are iteratively redefined following an initial determination of $\beta$ : each voxel acquires a weight dependent on its $\beta$ value, rewarding low $\beta$ “baseline” voxels and penalizing high $\beta$ “release” voxels (equation (8)),

$$w\left(\beta \right)=\frac{1}{1+\exp \left(2z\left(\beta \right)\right)}$$

(8)

where $z\left(\beta \right)$ is the z-score across all voxel $\beta$ values. Taking a weighted average of voxel TACs in each k-region with this weighting scheme results in new baseline TACs ${\overline{C}}_{s,k}^{\prime }\left(t\right)$ that are less contaminated by DA release-related signal, helping to improve detection of low-magnitude signals.

Data-driven predictors are then defined using $\beta$ values: voxel residuals are normalized to a peak value of 1, weighted by their ${\beta }_i$ value, and averaged across voxels. The result is low-pass filtered to reduce the impact of noise then renormalized to a peak of 1, providing a new predictor $P^{\prime }\left(t\right)$ (Figure 3).

From simulations it was observed that only one iterative update of this protocol was required to achieve predictions of baseline voxel $\beta$ values near-zero, indicating adequate removal of contaminant DA-related signal from the baseline prediction, thus better defining true-release voxels (Supplemental Figure 7).

Methods and materials

Simulation design

To test our method in comparison to the lp-ntPET method, we simulated 79-minute single-task RAC scans on the GE SIGNA PET/MR system (GE Healthcare, Chicago, IL, USA). A T1 MRI from a healthy human subject scanned on the SIGNA was segmented using FreeSurfer ³⁷ to generate an anatomical reference for simulating regions with pre-determined RAC kinetics; the segmented image was then resampled to the PET voxel dimensions (1.39 × 1.39 × 2.78 mm) using nearest-neighbor interpolation. The striatum was defined using bilateral ROIs from the putamen, caudate, and nucleus accumbens.

Noiseless TACs were simulated following the ntPET model, ⁹ using the same approach as Wang et al. ¹⁵ Briefly, the arterial input function was modeled as in literature ³⁸ to match concentration levels observed in 15 mCi RAC injections on the GE SIGNA PET/MR scanner at our centre. Kinetic parameters and modeling from literature ³⁹ were used alongside pre-specified DA concentration increases and timings to obtain simulated TACs.

DA release clusters were placed in the left caudate (cluster size: 118 voxels), in the right putamen (160 voxels), and in the right caudate (104 voxels) to simulate a possible response to a single 10-minute functional task block starting at 32 minutes, involving motor and cognitive processing. The peak simulated DA release magnitude was set to 100% increase from basal DA concentration (for full details see Supplemental Material and Supplemental Figure 1); this is relatively low compared to previous studies which simulated peak DA increases of 200% or more relative to baseline.^14,15

To ensure our method applies more generally, we then performed identical simulations in the presence of spatially heterogeneous RAC binding to assess the performance of the k-means methodology. Motivated by observed binding heterogeneities in subjects with PD (Supplemental Figure 2) we simulated a gradient in the binding potential (BP), accomplished by modulating B_max, where

$$BP=\frac{B_{\text{max}}}{K_D}$$

(9)

and $K_D$ is the equilibrium dissociation constant. ³⁹ The gradient was designed to approximate the direction and magnitude of observed binding gradients in human data (Supplemental Figure 2): in the putamen the simulated B_max ranged from 30 pmol/ml in the anterior and inferior regions to 70 pmol/ml in the superior and posterior regions, while B_max in the caudate ranged from 30–54 pmol/ml with similar directionality. In the homogenous binding simulation, a constant striatal value of ${\mathrm{B}}_{\max }=44\hspace{0.17em}$ pmol/ml was used.

The TACs were framed 4 × 1, 3 × 2, 4 × 5 minutes for estimation of the baseline portion of the scan, then 20 × 2 minutes for higher temporal sampling during DA release, and finally 1 × 4, and 2 × 5 minutes at the end of the scan for better counting statistics. The dynamic images were then spatially filtered with a 1.8 × 1.8 × 3.6 mm FWHM Gaussian to simulate the SIGNA point spread function (PSF). ⁴⁰ We generated 100 noisy realizations (NRs) using scaled Gaussian noise, ¹⁴ before post-processing with IHYPR4D (7.2 mm, two frames spatiotemporal kernel, one iteration) ³⁰ ; the scaling coefficient was set so that the noise level after applying IHYPR4D matched human RAC scan noise levels reconstructed with PSF-HYPR4D-K-TOFOSEM (5.6 mm, ten frames spatiotemporal kernel, 10 iterations) ³¹ and post-processed with IHYPR4D.

Implementation of RSD

$\beta$ maps output from RSD can be used in a workflow similar to fMRI, ³⁶ however in this work we performed a direct comparison with lp-ntPET by first creating binary detection maps via thresholding, and then comparing these across NRs. First, RSD was performed as described in the Theory section, and in-plane edge voxels were removed from the maps to minimize false positives caused by partial volume effects. To obtain binary activation maps, $\beta$ maps were variably thresholded at a common value across NRs; specifically, the variable threshold ranged between the maximum and minimum $\beta$ value across all NRs, yielding sensitivity results at all possible values of FPR. These results were then summarized using cluster-level receiver operating characteristic (ROC) curves.

Cluster size thresholding

To control for the multiple comparisons problem, clusters of predicted DA release voxels in the binary maps were subject to a cluster size threshold (CST), where clusters smaller than the CST were rejected as DA release clusters. ¹⁵ In our simulations the left putamen was kept devoid of release, therefore we used this region to define the CST. The CST was set for each simulated dataset and methodology independently and defined so that 90% of noisy realizations resulted in the left putamen being completely free of false positives.

Comparison to model-based lp-ntPET method

We compare our proposed RSD method to the existing model-based lp-ntPET method as described in literature. ¹⁵ Both MRTM and lp-ntPET (equations (1) and (2) were fit to each voxel TAC in the striatum. Parameters for lp-ntPET basis functions, $h\left(t\right)$ (equation (3)), were discretized as follows: (1) $t_D\in [t_{\mathit{start}}-5\hspace{0.17em}\min , t_{\mathit{start}}+5\hspace{0.17em}\min ]$ in one-minute intervals, where $t_{\mathit{start}}$ is the known task start time; (2) $t_p\in [t_D+1\hspace{0.17em}\min , t_{\mathit{end}}]$ in one-minute intervals, where $t_{\mathit{end}}$ is the known task end time; (3) $\alpha \in \{0.\mathrm{25}\text{,}1,4\}$ .

Weighted residual sum of squares (WRSS) was calculated for lp-ntPET and MRTM fits, which permitted the calculation of the F-statistic. Voxels with an F-statistic greater than a critical value were accepted as release voxels. The critical value was varied similarly to the thresholding performed in RSD to generate ROC curves for each simulated release cluster, and eventually chosen to have a matching FPR to that of RSD. Note that this is in contrast to previous studies with lp-ntPET that used a fixed threshold for the F-statistic,^14,15 corresponding to a theoretical p-value < 0.05 as defined by the degrees of freedom in the models. The purpose of variable thresholding employed in this work is to enable a comparison of detection sensitivity between the two methods at matched FPRs. F-statistic thresholding defined binary maps for each noisy realization, which then went through the same in-plane edge voxel removal and CST processing as RSD.

Performance metrics

To determine the performance of our method in comparison with the lp-ntPET method we used two sensitivity measures as defined in literature. ¹⁵ (1) Voxelwise sensitivity (VWS): for a given true positive voxel, this is the percentage of noisy realizations in which the voxel was correctly classified as a release voxel after CST was applied. Mean voxelwise sensitivity is defined as the mean of this quantity across all voxels within regions where DA release was simulated. This measure can be interpreted as the mean true positive rate, and spatial maps of voxelwise sensitivity were used to compare the ability of a given method to determine the accurate shape and extent of release within a region for a single scan. (2) Binary sensitivity (BS): for a given true positive region, this is the percentage of noisy realizations in which at least one voxel in the region survived the CST. This measure can be interpreted as the likelihood of the method to correctly classify at least some portion of a given true positive region in a single scan.

Human scans

Scan protocol

As proof-of-concept, we present results from a simple finger tapping motor task study. Nine healthy controls (mean age = 64.8, SD = 7.6, range = 47–73) were scanned on the GE SIGNA PET/MR scanner. All participants were right-handed. The UBC Research Ethics Board (following ethical standards of the Helsinki Declaration of 1975 and as revised in 1983) approved this study (certificate number: H19-03166), and written informed consent was obtained for each subject prior to all study procedures.

At 36 minutes post bolus injection of RAC (550 MBq; 14.86 mCi), subjects tapped fingers of their dominant hand for 10 minutes, in blocks of 2 minutes with 15 seconds of rest between blocks. The task involved serial tapping of the index finger, middle finger, ring finger, and little finger, then reversed (‘1-2-3-4, 4-3-2-1’…), at the subject’s own pace.

Image processing

Raw PET data were reconstructed using PSF-HYPR4D-K-TOFOSEM, ³¹ with attenuation, scatter, randoms, and normalization corrections, after which PET images underwent frame realignment (i.e. inter-frame motion correction). Dynamic PET images were denoised using IHYPR4D, with the same reconstruction and post-processing parameters as previously mentioned. Masks of the striatum and cerebellum were extracted from subject-specific FreeSurfer segmentations of structural T1 images for use in RSD and lp-ntPET. To create group-level images, subject T1 images were normalized to the Neurodevelopmental MRI Database (age 65-69) atlas with a 2 mm isotropic voxel size, ⁴¹ using ANTs normalization software. ⁴²

Image analysis

RSD and lp-ntPET were separately performed on subject-space dynamic PET images. RSD $\beta$ maps (units of % signal change) indicate relative deflections of voxel TACs from baseline. Lp-ntPET $\gamma$ maps (equation (2), units of s⁻¹) can be interpreted as the magnitude of the effect of DA on the apparent tissue efflux rate. ¹⁰ Results were then brought to a common space using the forward transformation from the normalization step. One sample t-tests were applied to RSD $\beta$ values and lp-ntPET F-statistics, providing group level statistical parametric maps for both methods.

Results

Simulations

Release in the presence of homogeneous binding

Figure 4 displays voxelwise prediction maps calculated from the set of binary classification maps (thresholded at 5% FPR) across 100 NRs for the lp-ntPET + F-test method vs. the proposed RSD method, in addition to ROC curves for the cognitive clusters (bilateral caudate) and motor cluster (right putamen). CST was determined for each method from the set of binary maps and was subsequently applied to each individual binary map; the resulting voxelwise sensitivity values after CST application are given in Table 1.

Figure 4.

Simulated single-task release – homogeneous RAC binding. Left: Simulated ground truth clusters are placed bilaterally in the caudate and in the right putamen. Center: Voxelwise prediction maps for the lp-ntPET + F-test method and RSD method applied to single-task simulations of DA release in the presence of homogeneous binding. Color bar represents percentage of noisy realizations in which the voxel was predicted as release, determined using 100 noisy realizations described in Methods. Right: ROC curves are shown for each type of cluster for the lp-ntPET + F-test method and RSD, prior to CST application. The ROC abscissa range was limited to clearly demonstrate differences among method results.

Table 1.

Comparison between RSD and lp-ntPET, with sensitivities and CST for both simulation sets.

Method sensitivities (5% FPR)	lpnt-PET + F-test					RSD
	CST	Cognitive		Motor		CST	Cognitive		Motor
		Mean VWS (%)	BS (%)	Mean VWS (%)	BS (%)		Mean VWS (%)	BS (%)	Mean VWS (%)	BS (%)
Single-Task (homo. binding)	30	24.4	79	37.8	80	30	51.0	99	77.0	99
Single-task (hetero. binding)	39	18.5	69	31.9	67	45	36.5	92	57.9	96

Binary maps were produced from $\beta$ (RSD) and F-statistic (lp-ntPET) maps by thresholding at 5% FPR over 100 noisy realizations to determine the CST. Mean VWS was then determined after application of the CST. CST: cluster size threshold; VWS: voxelwise sensitivity; FPR: false positive rate; BS: binary sensitivity.

As seen in the ROC curves, RSD drastically outperforms lp-ntPET across all FPRs for this low-magnitude single-task release ( $\sim 100\hspace{0.17em}$ pmol/ml peak increase). Binary sensitivity (BS) for lp-ntPET was 79% and 80% for the bilateral caudate and putamen clusters after application of CST, respectively; RSD had a BS of 99% for both clusters after CST.

Release in the presence of heterogeneous binding

Results for the heterogeneous binding simulations prior to application of CST are displayed in Figure 5, while Table 1 summarizes the heterogeneous binding simulation results after application of the CST. Voxelwise sensitivity is reduced in both methods and clusters when compared with the results for homogeneous binding, with the RSD method still considerably outperforming lp-ntPET. Binary sensitivity of RSD remains quite high in these clusters, resulting from relatively high peak $\beta$ values within the caudate visible in Figure 5.

Figure 5.

Simulated single-task release – heterogeneous RAC binding. Left: Simulated ground truth clusters are placed bilaterally in the caudate and in the right putamen. Center: Voxelwise prediction maps for the lp-ntPET + F-test method and RSD method applied to single-task simulations of DA release in the presence of heterogeneous binding. Color bar represents percentage of noisy realizations in which the voxel was predicted as release, determined using 100 noisy realizations described in Methods. Right: ROC curves are shown for each type of cluster for the lp-ntPET + F-test method and RSD, prior to CST application. The ROC abscissa range was limited to clearly demonstrate differences among method results.

Human scans

Results from the proof-of-concept study with 9 healthy control subjects are shown in Figure 6, displaying statistical parametric maps from RSD and lp-ntPET. Due to the intrinsically high variability of lp-ntPET $\gamma$ values (see Discussion and Supplemental Figure 10), we used voxel F-statistics for group statistical testing. This helped to provide voxels that passed low t-statistic thresholding (p < 0.01, uncorrected) from lp-ntPET. RSD on the other hand shows significant clusters at much higher t-statistic thresholding (p < 0.001, uncorrected) that appear physiologically plausible, which is consistent with the findings of our simulations.

Figure 6.

Motor task human results – statistical parametric maps. Subject-level $\beta$ values (RSD) and F-statistics (lp-ntPET) were submitted to one-sample t-tests, with maps displaying raw and thresholded t-statistics. P-values are uncorrected for multiple comparisons, similar to other work in the literature. ⁴⁴

Discussion

This work presents a novel detection framework for PET-based DA release detection, most advantageous for spatially localized and low-magnitude responses which are inherently difficult to detect. Such release may be expected in motor, cognitive, and/or neuropsychological reward-based tasks rather than pharmacological interventions,^11,27,28 where the spatial extent and magnitude of release is much greater. ⁴³ Given the increasing incidence of addiction and neurodegenerative diseases, which alter dopaminergic function, the ability to detect dopamine release in response to non-pharmacological stimuli is becoming increasingly more important.

Our detection framework has two major differences compared to existing methods which primarily incorporate significance testing between models to detect release:

Moving detection to a residual space: we leverage the condition of spatially localized, low-magnitude release in behavioral tasks to construct voxelwise hypothetical baseline TACs from regional averages of the striatum. Regressing these against measured voxel TACs prior to task initiation and computing residuals moves the analysis to a space where structured temporal signals are most likely indicative of DA release.

Decoupling kinetic modeling from detection: The transformation to residual space removes the effect of baseline kinetics; in TAC space both release and non-release TACs have structured signals, whereas only DA-release TACs should have structure in residual space, so modeling the residuals with a kinetic model is unnecessary for detection. Rather, using GLM to identify similarities between task-informed predictors and the residuals is sufficient.

We have demonstrated that the RSD method improves detection sensitivity of single-task protocols in the presence of both homogeneous and heterogeneous RAC binding. While we did not test these methods on non-denoised data, we can confidently claim that the denoising algorithms employed significantly aided DA release detection, as has been shown in previous work with HYPR-based denoising applied to voxelwise DA release detection. ¹⁵

Simulations

In this work we simulated two separate datasets: single-task release with homogeneous RAC binding and with heterogeneous RAC binding. In both simulations, RSD outperformed existing model-based methods employing MRTM/lp-ntPET fitting alongside an F-test. RSD provides higher binary and voxelwise sensitivities for clusters of low-amplitude DA release when compared to lp-ntPET at equivalent FPR, exhibiting improved robustness to noise and greater ability to distinguish release from non-release.

In the presence of positive noise spikes prior to task initiation and/or negative noise spikes during the task, lp-ntPET is especially prone to predicting release when it does not exist. In these scenarios lp-ntPET overfits to the noise-altered TAC by tweaking any of its four free parameters ( $R_1, k_2,k_{2a},\gamma$ ) and three discrete parameters ( $t_D,t_P,\alpha$ ) simultaneously to achieve the best fit to a given voxel. Although RSD can still overfit to noise, it performs its analysis in two steps which acts to de-couple several of its degrees of freedom. First, a voxel baseline TAC is predicted from the k-means regional mean TACs alongside the cerebellum TAC, dictated by four free parameters ( $w_{s,1},w_{s,2},w_{s,3},w_R$ ); though for a given voxel only one of these parameters typically contributes to its baseline prediction (see Supplemental Figure 5). Lastly GLM is performed on voxel residuals, governed initially by two fixed parameters defining predictors ( $t_D,L$ ) that are iteratively replaced by data-driven predictors, and one free parameter ( $\beta$ ). In total, RSD effectively uses two free parameters ( $w_{s,k},\beta$ ) for a given voxel in region $k$ , in a sequential manner, making it more robust to noise. Bias is indeed observed when the baseline fit is erroneous (e.g. when drastic changes in binding occur within a k-region; see Supplemental Figure 3), however these instances are readily identified and can be accounted for by increasing the number of $k$ -regions from $k=3$ or by comparison to the lp-ntPET result. To demonstrate differences in methodology between RSD and lp-ntPET, we provide a comparison of each method’s behavior for selected voxels in a single human scan in Supplemental Figure 6.

Heterogeneous binding

When observing scans of patients with PD, it was clear that the binding homogeneity assumption was not always satisfied (see Supplemental Material), and thus using a single mean striatal TAC as regressor in RSD was insufficient; when the mean striatal TAC (mid-BP) is fit pre-task to a relatively low-BP voxel TAC, the resulting baseline prediction will consistently overshoot the measured voxel TAC, resulting in large positive residuals post-task that are consistently predicted as DA release. To remedy this, we introduced a k-means step that supplies further baseline regressors in a data-driven manner, which partially accounted for this issue by removing these structured biases.

Results from RSD for single-task simulations with heterogeneous binding provided moderate to high detection sensitivity as shown in Figure 5 and Table 1 at 5% FPR, albeit with diminished performance relative to the homogeneous binding simulations. This was to be expected to some degree; as implemented, RSD made use of k = 3 basis functions to model a continuous gradient of binding, which results in a slight intra-regional bias shifted toward the lower binding voxels. We tested different values of k to provide further basis functions, however increasing k did not improve performance (see Supplemental Material). Overall, RSD still showed a 26% and 18% absolute increase in voxelwise sensitivity in the putamen and caudate clusters respectively, compared with lp-ntPET at 5% FPR.

Human scans

To test the performance of RSD relative to lp-ntPET, we performed scanning and analysis on a cohort of nine healthy controls. RSD found significantly more consistent results across subjects as demonstrated in Figure 6, while lp-ntPET showed highly variable output, with several subjects’ $\gamma$ values being contained near edge voxels of the striatum. The $\gamma$ parameter of lp-ntPET has been noted to suffer from several modelling issues, ²⁴ including over-determination of solutions resulting in highly variable $\gamma$ values and dependence of parameter magnitude on the choice of basis functions. These reasons likely contribute to the low statistical power observed in $\gamma$ maps (see Supplemental Figure 10), thus requiring use of the F-statistic when performing group comparisons for lp-ntPET.

Statistical parametric maps were computed by means of one-sample t-test (Figure 6), using the $\beta$ values from RSD and the F-statistics of lp-ntPET. The release pattern observed with RSD appears plausible, with a hand representation in the left putamen similar to prior work done with fMRI. ³⁴ Unlike the fMRI study, our results exhibit a bilateral release pattern (albeit with higher $\beta$ values in the left striatum). Bilateral release was also found in a finger tapping study using simultaneous PET/fMRI as analyzed by lp-ntPET along with a t-test, ⁴⁴ however the significance level of these results was limited to p < 0.05 (uncorrected), compared with RSD showing results at a more stringent p < 0.001 (uncorrected) in this study. This finding is consistent with the observed higher voxelwise and binary sensitivities from simulations, indicating that RSD finds and characterizes release clusters more consistently, leading to more robust group-level detection of low-magnitude release compared with lp-ntPET. Finally, all simulation and human scanning work made use of a single bolus injection of RAC, showing that the combination of advanced reconstruction and denoising methods alongside RSD enhances the clinical ease of investigations into dopamine function.

Limitations

Limitations of RSD include the possible presence of TAC shapes that cannot be described accurately by the data-derived baseline TAC regressors. This may occur in the case of extreme binding heterogeneities; though if such heterogeneities are observed, the k-means protocol may be adapted to attempt to account for its presence. Another possibility is that the localization assumption is violated for physiological tasks, and the majority of the striatum is activated by DA release such that the baseline regressors are unavoidably contaminated. In this scenario RSD may fail to detect the full extent of release and may be limited to detecting the regions of highest TAC deflection.

It has been demonstrated that [¹¹C]raclopride TACs may also be sensitive to task-induced blood flow changes, though this perturbation has opposite polarity and a different temporal shape to that from task-induced DA release, ³⁹ which would result in negative $\beta$ values. Our initial predictors were designed based on realistic simulations of DA-induced perturbations, while thresholding above a positive critical $\beta$ excludes voxels with negative residuals induced by task-induced blood flow changes.

Additionally, the predictors presented in this work are by construct only sensitive to task-induced DA release close to task initiation. This may be a limitation if there is a reason to suspect significant DA release before task initiation, such as potentially present in a task anticipatory stage. We simulated cases where the predictor timing is partially offset from the true release time (up to ±5 minutes) and found that RSD sensitivity diminishes, but still outperforms lp-ntPET + F-test (shown in Supplemental Figure 8). RSD somewhat accounts for this possibility in its iteratively updated predictor by incorporating residuals up to 5 minutes prior to task initiation, though careful evaluation of the TAC behavior and parameter choice needs to be performed to fully describe such behavior.

From our simulations, the accuracy of the shape of the predictor is less important for detection than accurate quantification; RSD assigns $\beta$  > 0 for any positively-structured residuals—so long as the predictor has reasonable overlap with the residuals. Simulated results demonstrate even a simple boxcar predictor coincident with the task timing only slightly reduces RSD sensitivity (Supplemental Figure 9). The iteration step, where the predictor shape is fine-tuned, only offers a minor improvement in sensitivity. In human scanning, subject-level residual response to task was well approximated by the model-based predictor (Supplemental Figure 11).

Conclusion

In the presence of low-magnitude DA release in response to a task-based single-bolus PET protocol we have shown that our proposed method, Residual Space Detection (RSD), is capable of achieving high detection sensitivity at low FPR, outperforming the traditional lp-ntPET + F-test method for localized, low-amplitude release. As shown using simulations and proof-of-concept human scanning, RSD achieves superior performance as a result of its robustness to noise and data-driven methodology. Future work will explore the extension of RSD to two-task release in attempts of differentiating function-specific (e.g. cognitive, reward) release from confounding motor release in cohorts of healthy and diseased individuals.

Supplemental Material

sj-pdf-1-jcb-10.1177_0271678X231214823 - Supplemental material for Novel voxelwise residual analysis of [¹¹C]raclopride PET data improves detection of low-amplitude dopamine release

Supplemental material, sj-pdf-1-jcb-10.1177_0271678X231214823 for Novel voxelwise residual analysis of [¹¹C]raclopride PET data improves detection of low-amplitude dopamine release by Connor WJ Bevington, Jordan U Hanania, Giovanni Ferraresso, Ju-Chieh (Kevin) Cheng, Alexandra Pavel, Dongning Su, A Jon Stoessl and Vesna Sossi in Journal of Cerebral Blood Flow & Metabolism