Spatiotemporal Kernel Reconstruction for Linear Parametric Neurotransmitter PET Kinetic Modeling in Motion Correction Brain PET of Awake Rats

The linear parametric neurotransmitter positron emission tomography (lp-ntPET) kinetic model can be used to detect transient changes (activation) in endogenous neurotransmitter levels. Preclinical PET scans in awake animals can be performed to investigate neurotransmitter transient changes. Here we use the spatiotemporal kernel reconstruction (Kernel) for noise reduction in dynamic PET, and lp-ntPET kinetic modeling. Kernel is adapted for motion correction reconstruction, applied in awake rat PET scans. We performed 2D rat brain phantom simulation using the ntPET model at 3 different noise levels. Data was reconstructed with independent frame reconstruction (IFR), IFR with HYPR denoising, and Kernel, and lp-ntPET kinetic parameters (k_2a: efflux rate, γ: activation magnitude, t_d: activation onset time, and t_p: activation peak time) were calculated. Additionally, significant activation magnitude (γ) difference with respect to a region with no activation (rest) was calculated. Finally, [¹¹C]raclopride experiments were performed in anesthetized and awake rats, injecting cold raclopride at 20 min after scan start to simulate endogenous neurotransmitter release. For simulated data at the regional level, IFR coefficient of variation (COV) of k_2a, γ, t_d and t_p was reduced with HYPR denoising, but Kernel showed the lowest COV (2 fold reduction compared with IFR). At the pixel level the same trend is observed for k_2a, γ, t_d and t_p COV, but reduction is larger with Kernel compared with IFR (10–14 fold). Bias in γ with respect with noise-free values was additionally reduced using Kernel (difference of 292, 72.4, and −6.92% for IFR, IFR+KYPR, and Kernel, respectively). Significant difference in activation between the rest and active region could be detected at a simulated activation of 160% for IFR and IFR+HYPR, and of 120% for Kernel. In rat experiments, lp-ntPET parameters have better confidence intervals using Kernel. In the γ, and t_d parametric maps, the striatum structure can be identified with Kernel but not with IFR. Striatum voxel-wise γ, t_d and t_p values have lower variability using Kernel compared with IFR and IFR+HYPR. The spatiotemporal kernel reconstruction adapted for motion correction reconstruction allows to improve lp-ntPET kinetic modeling noise in awake rat studies, as well as detection of subtle neurotransmitter activations.

Introduction

Transient changes in brain neurotransmitter levels can be investigated with dynamic positron emission tomography (PET) using for example the linear parametric neurotransmitter PET kinetic model (lp-ntPET) (Morris et al., 2005; Normandin et al., 2012). Transient changes in dopamine levels due to a rewarded task (Pappata et al., 2002), a motor planning task (Alpert et al., 2003), and due to psychosocial stress (Lataster et al., 2011), have been investigated in human brain PET. More recently, the effect of smoking (Cosgrove et al., 2014), gambling (Bevington et al., 2021), and cannabis (Calakos et al., 2021) on transient dopamine release has been investigated with the lp-ntPET method. These type of studies are not possible to perform in typical preclinical PET scans in which the animal is anesthetized. Therefore, the use of motion correction techniques, allowing to scan animals in the awake state (Kyme et al., 2014; Spangler-Bickell et al., 2016; Miranda et al., 2017), would make possible to investigate transient changes in neurotransmitter levels caused by a task or external stimuli in preclinical PET. Methods that involve head motion tracking followed by motion correction have been developed and improved over the last years to perform scans in awake rodents (Kyme et al., 2014; Spangler-Bickell et al., 2016; Miranda et al., 2017).

Using the lp-ntPET model (Normandin et al., 2012), the transient activation of certain neurotransmitter receptors can be quantified using tracers targeting these receptors (e.g., [¹¹C]raclopride for dopamine D_2/3 receptors) (Kyme et al., 2019). By modeling the tracer efflux in compartment modeling as a time varying parameter (Normandin et al., 2012), transient changes in endogenous neurotransmitter concentrations can be inferred by transient changes in tracer binding. For instance, the lp-ntPET has been used to quantify the striatal transient dopamine activation profile in awake rats following an amphetamine challenge (Kyme et al., 2019).

In order to perform kinetic modeling, dynamic PET reconstruction is necessary to determine the tracer concentration over time. Independent reconstruction of every time frame is the straightforward method to perform dynamic PET, but frame images, and therefore kinetic parameters, usually have high noise level due to the small number of events in each frame. To reduce noise in dynamic PET and kinetic modeling parameters, a wide variety of methods can be implemented (Reader and Verhaeghe, 2014; Wang et al., 2020), such as post-processing using the highly constrained backprojection method (Christian et al., 2010), or using machine learning denoising (Reader et al., 2021). Particularly for the case of dynamic PET for kinetic modeling, direct reconstruction has been developed to reduce noise (Matthews et al., 2010). In this method, the kinetic model is fitted to every voxel after every reconstruction iteration and therefore parametric images can be calculated during reconstruction. This method has been applied using the lp-ntPET kinetic model for noise reduction (Angelis et al., 2019). Another reconstruction developed for noise reduction in dynamic PET is the kernel method (Wang and Qi, 2015; Wang, 2019; Miranda et al., 2021). This method makes use of spatial and temporal correlations in the data to reduce noise in the iterative reconstruction (Wang and Qi, 2015; Novosad and Reader, 2016; Wang, 2019).

In this work, the lp-ntPET kinetic model was used to quantify transient dopamine changes in the rat striatum. As a first objective, we validated the spatiotemporal kernel reconstruction for lp-ntPET kinetic modeling in a 2D simulation and compare it with independent frame reconstruction. Then, we adapted the spatiotemporal kernel reconstruction for motion correction reconstruction to enable it in awake small animal scans. The method was used to perform a [¹¹C]raclopride scan, using cold raclopride as challenge, in an awake freely-moving rat using the point source tracking method (Miranda et al., 2017).

Materials and Methods

Motion Tracking and Independent Frame Motion Correction Reconstruction

The rat head motion in awake rat scans was tracked using the point source tracking method (Miranda et al., 2017). Four point sources prepared with [¹⁸F]FDG were attached on the rat head. Two point sources were attached below each ear, one on the nose bridge, and one in between the right ear and nose. Each point source was prepared with [¹⁸F]FDG and had an activity in the range of 222-370 kBq.

Animals were scanned on an Inveon PET scanner (Siemens Medical Solutions, Inc., Knoxville, United States). Images are reconstructed in a grid of 128 × 128 × 159 voxels with a size of 0.776 × 0.776 × 0.796 mm along the x, y and z directions, respectively. Independent frame motion correction reconstruction was calculated using list-mode event-by-event motion correction (LMMC) with 16 subsets and 8 iterations (Rahmim et al., 2008). The sensitivity image for motion correction was calculated by interpolation in the image space (Rahmim et al., 2008). The attenuation map was calculated using the binary image of the activity body outline with an uniform attenuation factor for soft tissue (0.096 cm^–1) (Angelis et al., 2013). Motion dependent and spatially variant resolution modeling was implemented as well (Miranda et al., 2020). Dynamic images were reconstructed with a framing of 12 frames × 10 s, 6 ×20 s, 2 ×60 s, and 27 ×120 s.

Spatiotemporal Kernel Reconstruction for Motion Correction

The spatiotemporal kernel method (Wang and Qi, 2015; Wang, 2019) has been adapted for the case of PET rigid motion correction reconstruction. Briefly, the original method to calculate the spatial kernel matrix consists of dividing the entire PET scan in 3 frames and use the voxel intensity values of these 3 frames as the feature of the corresponding voxel. Using a Gaussian radial kernel, the correlation between voxels is calculated, which serve as the values of the spatial kernel matrix elements. For the case of small animal brain motion correction, the region of the image outside the head of the animal is not corrected for motion, and therefore it can be affected by blurring motion. Therefore, we define a rectangular region enclosing the animal head to calculate the spatial correlations for voxels only in that region. We use a neighborhood of 9 ×9 ×9 voxels, and a threshold of 0.8 in the radial Gaussian kernel value, to calculate the spatial kernel matrix. Only the 48 closest nearest neighbors were considered to create the sparse spatial kernel matrix.

To calculate the temporal kernel matrix, originally it was proposed to use the sinogram as the frame feature to calculate the correlation between frames (Wang, 2019). Although is possible to perform sinogram rebinning to calculate the motion corrected sinogram (Rahmim et al., 2008), these sinograms often present gaps due to the position of the detectors after motion correction which do not overlap the sinogram space. These gaps can differ between frames, and the effect can be pronounced in small animal brain scans, in which the animal head can have a wide range of orientations. For this reason, we replaced the sinogram with an approximate LMMC reconstruction as the frame feature. LMMC reconstruction allows to consider all events after motion correction for reconstruction (i.e., no events are discarded due to falling out of the sinogram space), and the motion corrected sensitivity image corrects for non-uniformities due to motion compensation (Rahmim et al., 2008). We calculate the approximate LMMC of every frame considering 16 subsets and only one iteration, without attenuation correction or resolution modeling. This smooth reconstruction also allows to reduce the difference between voxel kinetics, and therefore improve correlation between frames to produce temporal basis functions that can model the different kinetics present in the image. The point sources are masked from every frame image before calculating the frames features correlation with the radial Gaussian kernel (Wang, 2019). To reduce noise in the temporal basis functions calculated from the correlation between frames, we filter every temporal basis function using a Gaussian filter with σ = ts/100, were ts = 15 is the size of the frame neighborhood to calculate the correlation with other frames.

Highly Constrained Backprojection Denoising

Since the highly constrained backprojection (HYPR) denoising (Christian et al., 2010) has been shown to improve parameter estimation in lp-ntPET kinetic modeling (Wang et al., 2017), we applied HYPR denoising to independent frame reconstruction dynamic images. The time averaged sum of dynamic frames was used as the composite image, and a 3 ×3 ×3 boxcar filter was used to perform HYPR filtering in dynamic frames (Christian et al., 2010).

2D ntPET Simulation

In order to validate the spatiotemporal reconstruction for kinetic modeling using lp-ntPET, we performed a 2D simulation of a brain phantom using the ntPET model (Morris et al., 2005; Normandin and Morris, 2008). Similar to Angelis et al. (2019), we simulated a rat brain phantom with the striatum structure, where the left striatum did not present endogenous neurotransmitter activation (rest region), while the right striatum was activated (active region). A reference region, necessary to perform lp-ntPET kinetic modeling, was considered as well. The time activity curves (TACs) of the reference, rest, and active regions were generated for a 60 min scan with the same parameters as in Angelis et al. (2019). These parameters consider an activation profile peak of 200% the basal dopamine level (reducing dopamine binding by 10%), an onset activation time of 20 min, with a peak time at 25 min. Figure 1 shows the phantom image, and the TACs from the different regions. The simulation was performed considering [¹¹C] raclopride to simulate tracer decay, and incorporating photon attenuation with a uniform attenuation factor for soft tissue for the entire head. The image had a size of 128 ×128 pixels with a pixel size of 0.776 × 0.776 mm. List-mode data frames were generated with the same framing used for dynamic image reconstruction (section motion tracking and independent frame motion correction reconstruction). Simulations with 10, 40, and 80 million counts were generated, considering 30 realizations per count level. Data was reconstructed with independent frame reconstruction and spatiotemporal kernel reconstruction with 300 iterations in both cases.

FIGURE 1

Figure 1. 2D phantom regions and time activity curves used to generate the dynamic events data in the respective regions.

A second set of simulations were performed with the same previously described phantom and ntPET model parameters, but at 5 different peak levels of activation: 120, 140, 160, 180, and 200% the basal level (Figure 2). Ten realizations were calculated per activation level, with 40 million counts in all cases.

FIGURE 2

Figure 2. Active region (and rest region for reference) time activity curves with (A) 120 %, (B) 140 %, (C) 160%, (D) 180 %, and (E) 200 % percent peak dopamine release with respect to baseline. (F) Activation profiles for every respective level of activation.

Awake Rat Cold Raclopride Brain Scans

In order to perform the injection of the tracer in the awake state, a catheter was initially implanted in the jugular vein (Feng et al., 2015) in 2 Wistar female rats (Janvier Labs). Surgery was performed under isoflurane anesthesia (5% for induction, 1.5% for maintenance). After surgery, rats were left to rest during one week, followed by 3 days of acclimatization inside the holder used to maintain the rats inside the scanner field of view (Figure 3B). Catheter was flushed with heparin solution for maintenance every day for one week, and 2–3 times per week afterwards. The experiments followed the European Ethics Committee recommendations (decree 86/609/CEE) and were approved by the Animal Experimental Ethical Committee of the University of Antwerp, Antwerp, Belgium (ECD 2016-89).

FIGURE 3

Figure 3. (A) Scanning time schedules for the awake cold raclopride scan. (B) Diagram of the rat inside the holder tube connected to the line swivel and (C) scanner setup showing the injection pump and swivel.

Two rats underwent a cold raclopride challenge scan, one under anesthesia (210 g) and the other in the awake state (197 g). For the scans under anesthesia, the rat was initially administered with isoflurane (5% for induction, 1.5% for maintenance) and placed on the scanner bed. At the start of the 60 min PET scan, the rat was administered with [¹¹C]raclopride (11.6 MBq, Molar activity, MA: 38.2 MBq/nmol) through the jugular vein catheter. Twenty minutes after the start of the scan, cold raclopride in 0.2 mL saline (1 mg/kg) was administered. This dose was chosen to observe a clear displacement of [¹¹C]raclopride (Wadenberg et al., 2000; Kyme et al., 2019). For the awake scan (Figure 3A), 20 min before the start of the scan, four point sources were attached on the rat head in the awake state. At the onset of the PET scan, [¹¹C]raclopride was administered through the jugular vein catheter (12.4 MBq, MA: 45.3 MBq/nmol). Twenty minutes after the start of the scan, cold raclopride (1 mg/kg) in 0.2 mL saline was administered through the jugular vein catheter.

Kinetic Modeling

In both simulations and experimental data, the lp-ntPET kinetic model (Normandin et al., 2012) was used to calculate the magnitude and time of the activation profile. The lp-ntPET kinetic model represents the tissue TAC as a function of the reference region TAC as (Alpert et al., 2003; Normandin et al., 2012):

$C_T\left(t\right)=R_1{C}_R\left(t\right)+k_2{\displaystyle {\int }_0^t}{C}_R\left(u\right)\mathit{d}u\mathrm{\hspace{1em}}-k_{2a}{\displaystyle {\int }_0^t}{C}_T\left(u\right)\mathit{d}u-\mathrm{\gamma }{B}_i(t)(1)$

where C_T is the activity in tissue, C_R the activity in the reference tissue (cerebellum in our case), R₁ the ratio of the delivery in tissue compared to the reference tissue, k₂ is the rate constant transfer from free compartment to plasma, k_2a is the apparent rate constant transfer from specific compartment to plasma, and γ is the magnitude of the activation response modeled with basis functions B_i(t):

$B_i(t)={\int }_0^t{C}_T(u)h_i(u)\mathit{d}u(2)$

where h_i(t) is modeled with a gamma variate function:

$h_i\left(t\right)={\left(\frac{t-t_d}{t_p-t_d}\right)}^{\mathrm{\alpha }}\text{exp}\left(\mathrm{\alpha }\left[1-\frac{t-t_d}{t_p-t_d}\right]\right)u(t-t_d)(3)$

where t_d is the delay time (from injection onset) at which the activation starts, t_p is the peak time of the activation, α determines the skewness of the activation, and u(t) is the Heaviside function. Every i-th h_i(t) function is calculated with a different combination of t_d, t_p and α parameters, with the following ranges: t_d ranged from 10 to 40 min, in intervals of 1.5 min, t_p depended on t_d and ranged from t_d to t_end− 5 min (t_end:scan end time) in intervals of 1.5 min, and α ranged from 0.5 to 3, in intervals of 0.5. A total of 2,394 basis functions were calculated using (2) and (3). The activation response profile (ARP) is reported as the percentage of change in baseline dopamine efflux (k_2a):

$\text{ARP}=100\times \frac{\mathrm{\gamma }{h}_i(t)}{k_{2a}}(4)$

Since either in simulation experiments or in animal scans, a decrease in [¹¹C] raclopride binding is expected, non-negative linear least-squares was used to calculate the set of parameters [R₁k₂k_2aγ] with all basis functions and selecting the solution with the minimum least squared error from all basis functions. Parameters [t_dt_pα] are obtained from the basis function which results in the minimum least squared error. If no prior information about the decrease/increase of dopamine is known, linear least squares (i.e., allowing positive and negative γ magnitude) should be used.

Data Analysis

From simulation data, the mean and standard deviation (SD) of the activation parameters of interest, i.e., k_2a, γ, t_d, and t_p, over all realizations at the 3 different count levels, was calculated for independent frame reconstruction (IFR), spatiotemporal kernel reconstruction (Kernel), and IFR with HYPR denoising (IFR+HYPR). The relative difference with respect to the noise-free parameters is calculated as well. Mean and SD parametric γ, t_d, and t_p maps are calculated for the different count levels.

For the simulation with different activation profile levels, the relative magnitude of γ, i.e., the ratio γ/k_2a, which indicates the activation magnitude relative to the baseline washout, was calculated for the active and rest region. A paired t-test was calculated between rest and active regions γ/k_2a ratios at the 5 different activation levels for IFR, Kernel, and IFR+HYPR.

For the experimental data cold raclopride scans, a single frame reconstruction was used to manually align the rat brain to an MRI template with delineated striatum and cerebellum regions. TACs were extracted from these regions using PMOD 3.6 (Pmod technologies, Zurich, Switzerland). Kinetic modeling was performed and regional striatum k_2a, γ, t_d, and t_p parameters were calculated. In addition, the approximate Bayesian computation (ABC) framework (Toni et al., 2009; Kyme et al., 2019; Fan et al., 2021) was used to calculate the confidence intervals of the parameters k_2a, γ, t_d, and t_p in experimental data. For the prior, using the parameters obtained from the solution of least squares error, we considered a uniform distribution, with limits at +-100% the best fit values, considering 10 million sampling trials. The tolerance, obtained by trial and error as in Fan et al. (2021) was set at 2 times the least square error fit. Finally, t_d, t_p, and γ and its t-statistic (γ/SE(γ)) voxel-wise parametric brain maps, were calculated.

Results

2D ntPET Simulation

Figure 4 shows the regional and a pixel-wise TACs from the active and rest regions, and the ARP from the respective TACs, for noise realizations with 80 million counts. For IFR at the regional level, TACs present little variation (mean coefficient of variation, COV, along the TAC: 1.20%), however, the ARP present larger variation (COV: 2.59%). The mean magnitude of the ARP is higher in the active region compared to the rest region and has a sharper profile. At the pixel level, TACs present high noise (COV: 10%), and ARP have large magnitude differences across realizations (COV: 22.8%). However, mean active ARP is higher than in the rest region, and peak time in active ARP coincides with the noise-free ARP. Denoising IFR with HYPR reduces variation in active region TACs and ARP, both at the regional (TAC: 0.94%, ARP: 1.72%) and pixel level (4.01%, 5.99%), with better resemblance in ARP shape compared to noise-free ARP, although with reduced magnitude. For the Kernel reconstruction at the regional level, both TAC and ARP variation is further reduced (0.40 and 0.90%, respectively) compared to IFR and IFR+HYPR. Moreover, magnitude of the rest ARP is reduced using Kernel compared to IFR. For Kernel at the pixel level, variation in the active TAC (2.2%) is similar to IFR variation at the regional level (1.19%). Noise in ARP is greatly reduced at the pixel level using Kernel (1.46%) compared to IFR (22.8%), and IFR+HYPR (5.99%), with lower variation at the regional level than IFR (2.59%) and IFR+HYPR (1.72%). Corresponding plots for 40 and 10 million counts simulations are shown in Supplementary Figures 1, 2, respectively. Noise increases with lower counts for all methods. At the pixel level for 10 million counts, the shape of the ARP is greatly distorted for IFR and improved in IFR+HYPR, but the mean peak time (32 min) still is close to the noise-free value. On the other hand, Kernel ARP are less distorted and with good peak time (26 min) correspondence with the noise-free ARP.

FIGURE 4

Figure 4. Individual and mean TACs for every realization, at the regional (A,B,E,F,I,J) and pixel level (C,D,G,H,K,L), and activation response profiles (ARP) for the active and rest regions, using independent frame reconstruction (IFR), independent frame reconstruction with HYPR denoising (IFR+HYPR), and spatiotemporal kernel reconstruction (Kernel), for simulations with 80 million counts. Dotted lines show 2 SD.

Table 1 show the lp-ntPET kinetic modeling parameters of interest for the simulations with all 3 count levels, for IFR, IFR+HYPR and Kernel. At the regional level, coefficient of variation is improved in IFR when HYPR denoising is used, with Kernel further reducing variability. On the other hand, relative difference of k_2a and γ with respect to the noise-free value is larger using IFR+HYPR and Kernel compared with IFR, but t_d and t_p relative difference is lower using Kernel compared to IFR and IFR+HYPR.

TABLE 1

Table 1. Regional lp-ntPET kinetic modeling k_2a, γ, t_d and t_p mean, coefficient of variation (COV), and difference with respect to noise-free value, for independent frame reconstruction (IFR), independent frame reconstruction with HYPR denoising (IFR+HYPR), and spatiotemporal kernel reconstructions (Kernel), in simulations with 80, 40 and 10 million counts.

Similar to Tables 1, 2 shows kinetic modeling statistics, but at the pixel level. As in the regional analysis, at the pixel level for all count levels, coefficient of variation is larger for IFR, with HYPR denoising reducing variability. Kernel shows the lowest COV for all parameters at all noise levels. Except for k_2a at 80 and 40 million counts, Kernel shows smaller difference with respect to noise-free values than IFR. Particularly, γ present large differences with respect to noise-free values using IFR.

TABLE 2

Table 2. Pixel-wise lp-ntPET kinetic modeling k_2a, γ, t_d and t_p mean, coefficient of variation (COV), and difference with respect to noise-free value, for independent frame reconstruction (IFR), independent frame reconstruction with HYPR denoising (IFR+HYPR), and spatiotemporal kernel reconstructions (Kernel), in simulations with 80, 40 and 10 million counts.

Figure 5 shows the mean and SD γ parametric maps for the different count levels using IFR and Kernel. For the rest region, as also shown in Table 2, magnitude of γ becomes larger using IFR and IFR+HYPR at increasing noise levels, while using Kernel, magnitude is similar across noise levels. Standard deviation is larger in the active region compared to the rest region pixels using IFR and IFR+HYPR, but lower using Kernel. For all count levels and both regions, standard deviation is lower using Kernel. Supplementary Figures 3, 4 show the mean and SD t_d and t_p parametric maps for the different count levels using IFR, IFR+HYPR, and Kernel. As with γ, t_d Kernel parametric maps are similar across noise levels and with lower standard deviation compared with IFR and IFR+HYPR in the active region, but IFR+HYPR reduce noise to a level similar to that in the Kernel reconstruction. On the other hand t_p parametric maps are similar using IFR, IFR+HYPR and Kernel, with Kernel presenting lower standard deviation in 80 and 40 million counts, but IFR+HYPR showing lower SD in 10 million counts in the active region. Indeed, t_p was the parameter with lowest COV for IFR and IFR+HYPR as shown in Table 2.

FIGURE 5

Figure 5. Rest and active parametric γ mean (top row) and standard deviation (bottom row) images, for independent frame reconstruction (IFR), independent frame reconstruction with HYPR denoising (IFR+HYPR), and spatiotemporal kernel reconstructions (Kernel).

Table 3 shows the relative magnitude of γ (γ/k_2a) for different peak activation levels, and the difference between active and rest region relative magnitude. Using IFR, no significant difference in γ/k_2a was found for peak activation levels of 120 and 140%, but with HYPR denoising difference at 120% becomes significant. Significant difference in IFR is found for 160, 180, and 200% peak activation levels, with highest significance in the 200% level. For Kernel, difference between active and rest γ/k_2a is significant for all activation peak levels (p**< 0.01), reaching a significance of p****< 0.0001 for the activation levels higher or equal to 140%. In addition, γ/k_2a increase in proportion to the activation peak level (120%: 0.088, 140%: 0.142, 160%: 0.199, 180%: 0.250, 200%: 0.295) using Kernel, while this is not observed using IFR (120%: 0.263, 140%: 0.159, 160%: 0.223, 180%: 0.282, 200%: 0.330) or IFR+HYPR.

TABLE 3

Table 3. Mean and standard deviation of the relative magnitude of γ in the active and rest regions, as well as their difference and significance, for peak levels of activation of 120, 140, 160, 180, and 200% the baseline level.

Awake Rat Cold Raclopride Brain Scans

Figure 6 shows the regional striatum and cerebellum TACs as well as the kinetic modeling fit and ARP, in both anesthetized and awake rats using IFR, IFR+HYPR, and Kernel. Noise is greater in Awake IFR, compared to Anesthesia IFR TACs, with both IFR+HYPR and Kernel reducing noise in both cases. ARP calculated from IFR, IFR+HYPR, and Kernel fits have very good correspondence, but IFR+HYPR shows reduced magnitude compared to IFR and Kernel. The ARP from awake rats have a sharper profile and larger relative magnitude compared with anesthesia ARP. Table 4 shows the parameters of the lp-ntPET kinetic modeling. Anesthesia scan parameters k_2a, and γ have good correspondence between IFR, IFR+HYPR and Kernel, but confidence intervals (95%) are smaller for values calculated using Kernel. Difference in awake scans k_2a, and γ are larger between IFR and Kernel (12 and 15% difference, respectively), with IFR+HYPR showing intermediate values. Timing parameters t_d, and t_p are the same using IFR, IFR+HYPR but slightly change using Kernel, for which confidence intervals are smaller.

FIGURE 6

Figure 6. Regional striatum and cerebellum TACs calculated from independent frame reconstruction (IFR), independent frame reconstruction with HYPR denoising (IFR+HYPR), and spatiotemporal kernel reconstructions (Kernel), in cold raclopride challenge scans of anesthetized and awake rats. Activation response profiles calculated from the lp-ntPET modeling are shown below. Cold raclopride injection at 20 min (red dotted line).

TABLE 4

Table 4. Regional striatum k_2a, γ, t_d, and t_p parameters calculated with the lp-ntPET kinetic model, using independent frame reconstruction (IFR), independent frame reconstruction with HYPR denoising (IFR+HYPR), and spatiotemporal kernel reconstructions (Kernel), for the anesthetized and awake rats.

Supplementary Figures 5, 6 show the posterior distribution of parameters estimated with ABC, in anesthesia and awake scans, respectively. In the anesthesia scan, compared to IFR, IFR+HYPR slightly decreases distribution spread, with Kernel showing distributions with the smallest spread in all parameters. In the awake scan the same trend is observed, where the peak of the distributions in timing parameters is more clear with the Kernel method.

Figure 7 shows the γ, γ/SE(γ) (γ t-statistic), and t_d parametric maps for the anesthesia scan, as well as striatum voxel-wise ARP, calculated using IFR, IFR+HYPR, and Kernel reconstructions. Figure 8 shows the equivalent figure for the awake scan. The shape of the striatum is not visible from the IFR anesthesia γ parametric map, showing high intensity values throughout the brain, while HYPR denoising slightly improves striatum structure. With Kernel the striatum structure is visible in the γ map, but also showing high intensity values outside the striatum. The IFR γ t-statistic parameter map shows lower intensity values outside the striatum, and some regions of high intensity within the striatum. HYPR denoising further improves striatum structure and reduces noise outside striatum. The same effect is observed in the Kernel γ t-statistic map, but showing a better striatum structure shape and slightly higher intensity values outside the striatum. The onset time t_d parametric map shows no striatum structure at the cold raclopride injection time value (20min, assigned to white color) using IFR and IFR+HYPR, while t_d maps using Kernel show good correspondence at the striatum region with the cold raclopride injection time. This is also shown in the ARP curves calculated from striatum voxels TACs, with IFR showing large variation across voxels, which is reduced using IFR+HYPR, and further reduced using Kernel. Supplementary Figure 7 shows histograms of γ, t_d, and t_p values considering striatum voxels. Distribution of values is larger in IFR histograms, especially in the t_d histogram where no clear distribution centered at a single value can be discerned. HYPR denoising reduces spread, and Kernel histograms show distributions with the lowest spread, with clear peaks in γ and t_d histograms (0.16 min^–1, and 20 min, respectively), close to the regional analysis values (0.16 min^–1, and 22 min).

FIGURE 7

Figure 7. Parametric γ, γ/SE(γ) (γ t-statistic), and t_dbrain maps, as well as activation response profiles (ARP) for striatum voxels, calculated using independent frame reconstruction (IFR, first row), independent frame reconstruction with HYPR denoising (IFR+HYPR, middle row), and spatiotemporal kernel reconstruction (Kernel, third row), for the anesthetized rat scan. White regions in t_d maps correspond to cold raclopride injection time (20 min). Striatum delineated from MRI template shown in red.

FIGURE 8

Figure 8. Parametric γ, γ/SE(γ) (γ t-statistic), and t_dbrain maps, as well as activation response profiles (ARP) for striatum voxels, calculated using independent frame reconstruction (IFR, first row), independent frame reconstruction with HYPR denoising (IFR+HYPR, middle row), and spatiotemporal kernel reconstruction (Kernel, third row), for the awake rat scan. White regions in t_d maps correspond to cold raclopride injection time (20 min). Striatum delineated from MRI template shown in red.

As in the anesthesia scan, in awake γ parametric maps using IFR the striatum structure cannot be identified, with HYPR denoising reducing noise outside the striatum. Using Kernel, the striatum structure is better defined, but also showing high values outside the structure. Calculating the γ t-statistic map reduces high intensity values outside the striatum for all methods, and the striatum structure is visible using IFR and IFR+HYPR, but shape is better defined using Kernel. The IFR and IFR+HYPR t_d parametric map show regions within the striatum with values close to the true cold raclopride injection time (white regions), but striatum shape is not well defined. On the other hand, Kernel t_d map white regions are well defined within the striatum. Voxel-wise striatum ARP variation is reduced using Kernel compared with IFR and IFR+HYPR. Supplementary Figure 8 shows γ, t_d, and t_p histograms considering striatum voxels. Histogram distributions are improved with HYPR denoising, but have less spread using Kernel, particularly in γ and t_p histograms. Peak value of t_p in the IFR and IFR+HYPR histograms (30 min) shows a large difference with respect to the region analysis value (38.5 min), while Kernel t_p histogram peak value has a closer value (38 min).

Discussion

The spatiotemporal kernel method has been validated for lp-ntPET kinetic modeling and adapted in conjunction with motion correction reconstruction, reducing noise in dynamic reconstructions and in kinetic parameters. Depending on the level of noise, subtle transient changes in neurotransmitter levels might be difficult to detect using regular independent frame reconstruction. The spatiotemporal kernel reconstruction reduces noise in the dynamic PET data, therefore allowing one to detect changes in the time activity curves and calculate kinetic parameters with less uncertainty. This was validated in simulation experiments, and then applied in a real data awake rat experiment.

For simulations data, at the pixel level Kernel produce reconstructions TACs with similar noise to IFR TACs at the regional level. Some negative bias is present in γ, at the regional and pixel level using IFR+HYPR and Kernel, which can be observed in the ARP lower magnitude compared to noise-free ARP (Figure 4). This could be caused by some smoothing effect of the HYPR filtering and Kernel reconstruction in the TAC. For HYPR this could be caused by the mismatch between the composite image and the lower intensity voxels in the release frames, while for Kernel the contrast can be adjusted by fine-tuning the spatial kernel matrix threshold (Wang and Qi, 2015). However, timing parameters (t_d and t_p) have excellent statistics (low coefficient of variation and bias) in both regional and pixel-wise data using Kernel. Although variability in the TAC magnitude increases with noise level using Kernel, the overall shape is well preserved, which might be an important factor in calculating accurate lp-ntPET kinetic modeling timing parameters. On the other hand, noise in IFR TACs can be overfitted using lp-ntPET and wrongly interpreted as an activation by the model. This is shown in the high variability in the IFR ARP parameters at the pixel level.

Similarly, parametric γ, t_d and t_p maps present lower noise using Kernel compared with IFR and IFR+HYPR, and have less variation across noise levels. But HYPR denoising presents closer performance to Kernel. At reduced count levels magnitude of γ parametric maps increase with respect to high count levels, but this effect is less pronounced using Kernel compared with IFR and IFR+HYPR. Activation onset time t_d parametric maps calculated using IFR present larger bias with increasing noise, but have stable values across different noise levels using Kernel. For peak time t_p parametric maps, mean values are similar between IFR, IFR+HYPR and Kernel maps, but noise is lower using IFR+HYPR and Kernel. As also observed in the pixel-wise analysis, this indicates that peak time t_p might be the most robust parameter calculated with lp-ntPET kinetic modeling.

The spatiotemporal kernel reconstruction also allows to detect more subtle activation profiles, as observed in the simulations with different ARP magnitudes. At the lowest activation of 120% the baseline level, kinetic modeling with lp-ntPET using IFR+HYPR and Kernel produced relative activation magnitudes (γ/k_2a) in the active region significantly different from the rest region, but IFR+HYPR failed to show significant differences at 140% the baseline level. The significance and magnitude of the difference increased with largest activation peak value using Kernel. Using IFR on the other hand, only for activations of 160% baseline and larger, a significant difference in the relative activation magnitude was observed between rest and active regions. Reduced noise in the Kernel reconstruction would therefore allow to detect neurotransmitter release of less intensity compared with IFR.

In anesthesia and awake cold raclopride scans, there is good correspondence between IFR, IFR+HYPR and Kernel lp-ntPET parameters at the regional level, with less noisy TACs calculated using Kernel. k_2a and γ parameters confidence intervals calculated using ABC (Fan et al., 2021) are smaller using Kernel compared with IFR and IFR+HYPR due to the TACs lower noise. For all methods, the activation due to the cold raclopride challenge has larger relative magnitude and faster rise in the awake compared to the anesthesia scan. This can also be observed in the striatum TACs where the activity level reaches the cerebellum level at an earlier time in the awake scan compared to the anesthesia scan. Due to the relatively large volume of the rat striatum structure (0.043 cm³ single side), regional analysis performs well using IFR. However, for smaller structures, or for studies with lower activity injection, the spatiotemporal kernel method would show larger differences compared with independent frame reconstruction (Miranda et al., 2021).

At the voxel level, performance of IFR is suboptimal, with HYPR denoising improving performance. The striatum structure cannot be identified in the IFR and IFR+HYPR activation magnitude γ parametric maps, but it is visible in the Kernel γ parametric maps, although large intensity values are still present outside the striatum, in both anesthesia and awake scans. Calculating the γ t-statistic helps to better identify the striatum structure in parametric maps using all methods, reducing high intensity regions outside the striatum, but striatum shape is better defined using Kernel. Activation onset time t_d parametric maps also show no striatum structure using IFR and IFR+HYPR, but it is well identified using Kernel by looking at the voxels with values close to the true cold raclopride injection time (20 min). Looking at the striatum voxels ARP, variation is large across voxels using IFR, and reduced with HYPR denoising, but more consistent ARP are obtained using Kernel, and consequently more consistent kinetic parameters are obtained. This is also shown in the γ, t_d, and t_p histogram plots with less spread values using Kernel compared with IFR and IFR+HYPR.

The spatiotemporal kernel reconstruction benefits from the relatively large rat striatum size spanning several voxels (about 100), which allows calculation of spatial basis functions also spanning several voxels. In addition, the temporal basis functions spanning only a finite time interval allow to preserve the shape of the TAC, including the transient changes. Temporal noise in the IFR voxel-wise TACs is detrimental for the lp-ntPET kinetic modeling since noise can be overfitted and wrongly interpreted as an activation in the TAC. Additional incorporation of HYPR denoising in calculation of the kernel matrix could further improve performance of the Kernel method (Cheng et al., 2021; Miranda et al., 2021).

Future work involves using the spatiotemporal kernel reconstruction in studies that could present subtle activations, for example in drug challenge (Kyme et al., 2019) or behavioral studies (Koepp et al., 1998; Lataster et al., 2011).

Conclusion

The spatiotemporal kernel reconstruction (Kernel) has been validated for lp-ntPET kinetic modeling and adapted for PET brain motion correction reconstruction in freely moving animals. In simulation experiments, Kernel improves kinetic parameters noise and bias, and allows to detect neurotransmitter activations of lower magnitude, compare with independent frame reconstruction (IFR) and HYPR denoising. In anesthetized and awake rat experiments, Kernel produce lp-ntPET parametric maps with better definition of the striatum structure and with more consistent parameter values across striatum voxels compared with IFR. Noise reduction using Kernel allows to perform neurotransmitter activation studies with lower parameters noise, and with detection of subtle neurotransmitter activations.

Data Availability Statement

The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

Ethics Statement

The animal study was reviewed and approved by Animal Experimental Ethical Committee of the University of Antwerp.

Author Contributions

AM was involved on the experimental design, software writing, data analysis and manuscript drafting and editing. JV was involved in the experimental design, data analysis, and manuscript drafting. DB, SiS, and StS were involved in drafting and editing the manuscript and figures. All authors approved the final manuscript and they are accountable for the content of the work.

Funding

AM was supported by the ERA-NET NEURON project SleepLess supported by BMBF (01EW1808) and FWO under the frame of Neuron Cofund. DB was supported by the Research Foundation Flanders (FWO) through a post-doctoral fellowship (1229721N). The University of Antwerp also funded the work through a partial assistant professor position for JV and a full professor position for StS.

Conflict of Interest

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.

Publisher’s Note

All claims expressed in this article are solely those of the authors and do not necessarily represent those of their affiliated organizations, or those of the publisher, the editors and the reviewers. Any product that may be evaluated in this article, or claim that may be made by its manufacturer, is not guaranteed or endorsed by the publisher.

Supplementary Material

The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fnins.2022.901091/full#supplementary-material

References

Alpert, N. M., Badgaiyan, R. D., Livni, E., and Fischman, A. J. (2003). A novel method for noninvasive detection of neuromodulatory changes in specific neurotransmitter systems. Neuroimage 19, 1049–1060. doi: 10.1016/s1053-8119(03)00186-1

PubMed Abstract | CrossRef Full Text | Google Scholar

Angelis, G., Bickell, M., Kyme, A., Ryder, W., Zhou, L., Nuyts, L., et al. (2013). “Calculated attenuation correction for awake small animal brain PET studies,” in Proceedings of the 2013 IEEE Nuclear Science Symposium and Medical Imaging Conference (Nss/Mic), (Piscataway, NJ: IEEE).

Google Scholar

Angelis, G. I., Gillam, J. E., Ryder, W. J., Fulton, R. R., and Meikle, S. R. (2019). Direct estimation of voxel-wise neurotransmitter response maps from dynamic PET data. IEEE Trans. Med. Imaging 38, 1371–1383. doi: 10.1109/TMI.2018.2883756

PubMed Abstract | CrossRef Full Text | Google Scholar

Bevington, C. W., Cheng, J. K., Klyuzhin, I. S., Cherkasova, M. V., Winstanley, C. A., and Sossi, V. (2021). A Monte Carlo approach for improving transient dopamine release detection sensitivity. J. Cereb. Blood Flow Metab. 41, 116–131. doi: 10.1177/0271678X20905613

PubMed Abstract | CrossRef Full Text | Google Scholar

Calakos, K. C., Liu, H., Lu, Y., Anderson, J. M., Matuskey, D., Nabulsi, N., et al. (2021). Assessment of transient dopamine responses to smoked cannabis. Drug Alcohol Depend. 227:108920. doi: 10.1016/j.drugalcdep.2021.108920

PubMed Abstract | CrossRef Full Text | Google Scholar

Cheng, J. K., Bevington, C., Rahmim, A., Klyuzhin, I., Matthews, J., Boellaard, R., et al. (2021). Dynamic PET image reconstruction utilizing intrinsic data-driven HYPR4D denoising kernel. Med. Phys. 48, 2230–2244. doi: 10.1002/mp.14751

PubMed Abstract | CrossRef Full Text | Google Scholar

Christian, B. T., Vandehey, N. T., Floberg, J. M., and Mistretta, C. A. (2010). Dynamic PET denoising with HYPR processing. J. Nucl. Med. 51, 1147–1154. doi: 10.2967/jnumed.109.073999

PubMed Abstract | CrossRef Full Text | Google Scholar

Cosgrove, K. P., Wang, S., Kim, S. J., McGovern, E., Nabulsi, N., Gao, H., et al. (2014). Sex differences in the brain’s dopamine signature of cigarette smoking. J. Neurosci. 340, 16851–16855. doi: 10.1523/JNEUROSCI.3661-14.2014

PubMed Abstract | CrossRef Full Text | Google Scholar

Fan, Y., Emvalomenos, G., Grazian, C., and Meikle, S. R. (2021). PET-ABC: fully Bayesian likelihood-free inference for kinetic models. Phys. Med. Biol. 66:115002. doi: 10.1088/1361-6560/abfa37

PubMed Abstract | CrossRef Full Text | Google Scholar

Feng, J., Fitz, Y., Li, Y., Fernandez, M., Puch, I. Cortes, Wang, D., et al. (2015). Catheterization of the carotid artery and jugular vein to perform hemodynamic measures, infusions and blood sampling in a conscious rat model. J. Vis. Exp. 5:51881. doi: 10.3791/51881

PubMed Abstract | CrossRef Full Text | Google Scholar

Koepp, M. J., Gunn, R. N., Lawrence, A. D., Cunningham, V. J., Dagher, A., Jones, T., et al. (1998). Evidence for striatal dopamine release during a video game. Nature 393, 266–268. doi: 10.1038/30498

PubMed Abstract | CrossRef Full Text | Google Scholar

Kyme, A., Se, S., Meikle, S., Angelis, G., Ryder, W., Popovic, K., et al. (2014). Markerless motion tracking of awake animals in positron emission tomography. IEEE Trans. Med. Imaging 331, 2180–2190. doi: 10.1109/TMI.2014.2332821

PubMed Abstract | CrossRef Full Text | Google Scholar

Kyme, A. Z., Angelis, G. I., Eisenhuth, J., Fulton, R. R., Zhou, V., Hart, G., et al. (2019). Open-field PET: simultaneous brain functional imaging and behavioural response measurements in freely moving small animals. Neuroimage 188, 92–101. doi: 10.1016/j.neuroimage.2018.11.051

PubMed Abstract | CrossRef Full Text | Google Scholar

Lataster, J., Collip, D., Ceccarini, J., Haas, D., Booij, L., van Os, J., et al. (2011). Psychosocial stress is associated with in vivo dopamine release in human ventromedial prefrontal cortex: a positron emission tomography study using [F-18]fallypride. Neuroimage 58, 1081–1089. doi: 10.1016/j.neuroimage.2011.07.030

PubMed Abstract | CrossRef Full Text | Google Scholar

Matthews, J. C., Angelis, G. I., Kotasidis, F. A., Markiewicz, P. J., and Reader, A. J. (2010). “Direct reconstruction of parametric images using any spatiotemporal 4D image based model and maximum likelihood expectation maximisation,” in Proceedings of the 2010 IEEE Nuclear Science Symposium Conference Record (Nss/Mic), (Piscataway, NJ: IEEE), 2435–2441.

Google Scholar

Miranda, A., Bertoglio, D., Stroobants, S., Staelens, S., and Verhaeghe, J. (2021). Low activity [(11)C]raclopride kinetic modeling in the mouse brain using the spatiotemporal kernel method. Phys. Med. Biol. 66:115005. doi: 10.1088/1361-6560/abfbf0

PubMed Abstract | CrossRef Full Text | Google Scholar

Miranda, A., Staelens, S., Stroobants, S., and Verhaeghe, J. (2017). Fast and accurate rat head motion tracking with point sources for awake brain PET. IEEE Trans. Med. Imaging 36, 1573–1582. doi: 10.1109/TMI.2017.2667889

PubMed Abstract | CrossRef Full Text | Google Scholar

Miranda, A., Staelens, S., Stroobants, S., and Verhaeghe, J. (2020). Motion dependent and spatially variant resolution modeling for PET rigid motion correction. IEEE Trans. Med. Imaging 39, 2518–2530. doi: 10.1109/TMI.2019.2962237

PubMed Abstract | CrossRef Full Text | Google Scholar

Morris, E. D., Yoder, K. K., Wang, C., Normandin, M. D., Zheng, Q. H., Mock, B., et al. (2005). ntPET: a new application of PET imaging for characterizing the kinetics of endogenous neurotransmitter release. Mol. Imaging 4, 473–489.

PubMed Abstract | Google Scholar

Normandin, M. D., and Morris, E. D. (2008). Estimating neurotransmitter kinetics with ntPET: a simulation study of temporal precision and effects of biased data. Neuroimage 39, 1162–1179. doi: 10.1016/j.neuroimage.2007.09.046

PubMed Abstract | CrossRef Full Text | Google Scholar

Normandin, M. D., Schiffer, W. K., and Morris, E. D. (2012). A linear model for estimation of neurotransmitter response profiles from dynamic PET data. Neuroimage 59, 2689–2699. doi: 10.1016/j.neuroimage.2011.07.002

PubMed Abstract | CrossRef Full Text | Google Scholar

Novosad, P., and Reader, A. J. (2016). MR-guided dynamic PET reconstruction with the kernel method and spectral temporal basis functions. Phys. Med. Biol. 612, 4624–4645. doi: 10.1088/0031-9155/61/12/4624

PubMed Abstract | CrossRef Full Text | Google Scholar

Pappata, S., Dehaene, S., Poline, J. B., Gregoire, M. C., Jobert, A., Delforge, J., et al. (2002). In vivo detection of striatal dopamine release during reward: a PET study with [(11)C]raclopride and a single dynamic scan approach. Neuroimage 16, 1015–1027. doi: 10.1006/nimg.2002.1121

PubMed Abstract | CrossRef Full Text | Google Scholar

Rahmim, A., Dinelle, K., Cheng, J. C., Shilov, M. A., Segars, W. P., Lidstone, S. C., et al. (2008). Accurate event-driven motion compensation in high-resolution PET incorporating scattered and random events. IEEE Trans. Med. Imaging 27, 1018–1033. doi: 10.1109/TMI.2008.917248

PubMed Abstract | CrossRef Full Text | Google Scholar

Reader, A. J., Corda, G., Mehranian, A., da Costa-Luis, C., Ellis, S., and Schnabel, J. A. (2021). Deep learning for PET image reconstruction. IEEE Trans. Radiat. Plasma Med. Sci. 5, 1–25.

Google Scholar

Reader, A. J., and Verhaeghe, J. (2014). 4D image reconstruction for emission tomography. Phys. Med. Biol. 592, R371–R418. doi: 10.1088/0031-9155/59/22/R371

PubMed Abstract | CrossRef Full Text | Google Scholar

Spangler-Bickell, M. G., Zhou, L., Kyme, A. Z., De Laat, B., Fulton, R. R., and Nuyts, J. (2016). Optimising rigid motion compensation for small animal brain PET imaging. Phys. Med. Biol. 619, 7074–7091. doi: 10.1088/0031-9155/61/19/7074

PubMed Abstract | CrossRef Full Text | Google Scholar

Toni, T., Welch, D., Strelkowa, N., Ipsen, A., and Stumpf, M. P. (2009). Approximate Bayesian computation scheme for parameter inference and model selection in dynamical systems. J. R. Soc. Interface 6, 187–202. doi: 10.1098/rsif.2008.0172

PubMed Abstract | CrossRef Full Text | Google Scholar

Wadenberg, M. L., Kapur, S., Soliman, A., Jones, C., and Vaccarino, F. (2000). Dopamine D2 receptor occupancy predicts catalepsy and the suppression of conditioned avoidance response behavior in rats. Psychopharmacology 150, 422–429. doi: 10.1007/s002130000466

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, G., Rahmim, A., and Gunn, R. N. (2020). PET parametric imaging: past, present, and future. IEEE Trans. Radiat. Plasma Med. Sci. 4, 663–675. doi: 10.1109/trpms.2020.3025086

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, G. B. (2019). High temporal-resolution dynamic PET image reconstruction using a new spatiotemporal kernel method. IEEE Trans. Med. Imaging 38, 664–674. doi: 10.1109/TMI.2018.2869868

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, G. B., and Qi, J. Y. (2015). PET image reconstruction using kernel method. IEEE Trans. Med. Imaging 34, 61–71.

Google Scholar

Wang, S., Kim, S., Cosgrove, K. P., and Morris, E. D. (2017). A framework for designing dynamic lp-ntPET studies to maximize the sensitivity to transient neurotransmitter responses to drugs: application to dopamine and smoking. Neuroimage 146, 701–714. doi: 10.1016/j.neuroimage.2016.10.019

PubMed Abstract | CrossRef Full Text | Google Scholar