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 (k2a: efflux rate, γ: activation magnitude, td: activation onset time, and tp: activation peak time) were calculated. Additionally, significant activation magnitude (γ) difference with respect to a region with no activation (rest) was calculated. Finally, [11C]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 k2a, γ, td and tp 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 k2a, γ, td and tp 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 td parametric maps, the striatum structure can be identified with Kernel but not with IFR. Striatum voxel-wise γ, td and tp 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.

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, [ 11 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 noisefree 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 γ,

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 , 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., [ 11 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 postprocessing 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 . 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 [ 11 C]raclopride scan, using cold raclopride as challenge, in an awake freely-moving rat using the point source tracking method (Miranda et al., 2017).

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 [ 18 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 [ 18 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 nonuniformities 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 [ 11 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.
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.

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).
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 [ 11 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 [ 11 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, [ 11 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): where C T is the activity in tissue, C R the activity in the reference tissue (cerebellum in our case), R 1 the ratio of the delivery in tissue compared to the reference tissue, k 2 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): where h i (t) is modeled with a gamma variate function: 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 ): Since either in simulation experiments or in animal scans, a decrease in [ 11 C] raclopride binding is expected, nonnegative linear least-squares was used to calculate the set of parameters R 1 k 2 k 2a γ with all basis functions and selecting the solution with the minimum least squared error from all basis functions. Parameters t d t 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 noisefree 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. 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 noisefree ARP. 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.

2D ntPET Simulation
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. 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. 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.  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.   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.  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.
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).
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).
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), Calculated from simulations with 40 million counts and considering 10 realizations per activation level. Reconstructed with independent frame reconstruction (IFR), independent frame reconstruction with HYPR denoising (IFR+HYPR), and spatiotemporal kernel reconstructions (Kernel). n.s: not significant.
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. Uncertainty is calculated using the ABC algorithm.
Frontiers in Neuroscience | www.frontiersin.org FIGURE 7 | Parametric γ, γ/SE(γ) (γ t-statistic), and t d brain 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.
but striatum shape is not well defined. On the other hand, Kernel t d map white regions are well defined within the striatum. Voxelwise 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 FIGURE 8 | Parametric γ, γ/SE(γ) (γ t-statistic), and t d brain 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.
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 3 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 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.