Incorporation of anatomical MRI knowledge for enhanced mapping of brain metabolism using functional PET

Functional positron emission tomography (fPET) imaging using continuous infusion of [18F]-fluorodeoxyglucose (FDG) is a novel neuroimaging technique to track dynamic glucose utilization in the brain. In comparison to conventional static or dynamic bolus PET, fPET maintains a sustained supply of glucose in the blood plasma which improves sensitivity to measure dynamic glucose changes in the brain, and enables mapping of dynamic brain activity in task-based and resting-state fPET studies. However, there is a trade-off between temporal resolution and spatial noise due to the low concentration of FDG and the limited sensitivity of multi-ring PET scanners. Images from fPET studies suffer from partial volume errors and residual scatter noise that may cause the cerebral metabolic functional maps to be biased. Gaussian smoothing filters used to denoise the fPET images are suboptimal, as they introduce additional partial volume errors. In this work, a post-processing framework based on a magnetic resonance (MR) Bowsher-like prior was used to improve the spatial and temporal signal to noise characteristics of the fPET images. The performance of the MR guided method was compared with conventional Gaussian filtering using both simulated and in vivo task fPET datasets. The results demonstrate that the MR-guided fPET framework denoises the fPET images and improves the partial volume correction, consequently enhancing the sensitivity to identify brain activation, and improving the anatomical accuracy for mapping changes of brain metabolism in response to a visual stimulation task. The framework extends the use of functional PET to investigate the dynamics of brain metabolic responses for faster presentation of brain activation tasks, and for applications in low dose PET imaging.

is a novel neuroimaging technique to track dynamic glucose utilization in the brain. In comparison to conventional static or dynamic bolus PET, fPET maintains a sustained supply of glucose in the blood plasma which improves sensitivity to measure dynamic glucose changes in the brain, and enables mapping of dynamic brain activity in task-based and resting-state fPET studies. However, there is a trade-off between temporal resolution and spatial noise due to the low concentration of FDG and the limited sensitivity of multi-ring PET scanners. Images from fPET studies suffer from partial volume errors and residual scatter noise that may cause the cerebral metabolic functional maps to be biased. Gaussian smoothing filters used to denoise the fPET images are suboptimal, as they introduce additional partial volume errors. In this work, a post-processing framework based on a magnetic resonance (MR) Bowsher-like prior was used to improve the spatial and temporal signal to noise characteristics of the fPET images. The performance of the MR guided method was compared with conventional denosing methods using both simulated and in vivo task fPET datasets. The results demonstrate that the MR-guided fPET framework denoises the fPET images and improves the partial volume correction, consequently enhancing the sensitivity to identify brain activation, and improving the anatomical accuracy for mapping changes of brain metabolism in response to a visual stimulation task. The framework extends the use of functional PET to investigate the dynamics of brain metabolic responses for faster presentation of brain activation tasks, and for applications in low dose PET imaging.

Introduction
Brain imaging using positron emission tomography (PET) can provide unique insights into brain function in both healthy individuals and individuals with neuropathological conditions ( Nasrallah and Dubroff, 2013 ). [18F]-fluorodeoxyglucose (FDG)-PET imaging has long been a proxy for regional and global brain metabolism, as glucose uptake is closely correlated with the underlying neuronal activity ( Figley and Stroman, 2011 ;Phelps et al., 1979 ;Reivich et al., 1985 ). Conventional static FDG-PET based on a bolus injection of the radiotracer provides a snapshot of glucose metabolism over a long time-window (equal to the scan duration, usually 10-30 minutes). Dynamic PET imaging using a bolus administration of radiotracer provides an opportunity to of fPET remains substantially lower than that of fMRI, which is in the order of seconds or even sub seconds. The current temporal resolution of fPET (around 20-60 seconds) limits the opportunity to use fPET for detailed investigations of brain metabolic responses to rapidly switching tasks and brain stimulation paradigms.
Analysis of fPET data is challenging because of the relatively poor signal to noise ratio (SNR) and partial volume errors in the reconstructed PET images . Recent work has improved the SNR in fPET by applying a combined bolus and continuous infusion of radiotracer during experiments ( Jamadar et al., 2019 ;Rischka et al., 2018 ). However, the statistical power of these experimental approaches is still relatively low when compared with fMRI. To mitigate this issue, spatial smoothing of the reconstructed PET images is performed prior to functional analysis of the brain using techniques such as independent component analysis (ICA). Gaussian smoothing is widely used as a postreconstruction spatial and temporal smoothing operation for functional neuroimaging analyses ( Chen and Calhoun, 2018 ;Hahn et al., 2018 ;Jamadar et al., 2019 ;Pignat et al., 2013 ;Villien et al., 2014 ). However, the Gaussian kernel acts as a low-pass filter, and therefore, further worsens the partial volume errors in fPET images; this can cause errors in the localisation and quantification of brain functional activations and at high-temporal resolution fPET imaging. MRI-based PET reconstruction methods have shown substantial improvement in PET image quality compared to conventional methods Sudarshan et al., 2020Sudarshan et al., , 2018. For instance, several studies have explored post-reconstruction PET image enhancement using anatomical information from structural MRI ( Bousse et al., 2012 ;Hutton et al., 2013 ;Schramm et al., 2018 ) to perform partial volume correction and image deblurring ( Dutta et al., 2013 ;Song et al., 2019 ). These methods propose a regularized deconvolution problem formulation using MRI information to address the PET denoising and partial volume error problems. Works involving anatomical information as regularization terms include penalties that model the statistical dependencies across the PET and MRI images, such as joint entropy regularization. ( Vunckx and Nuyts, 2010 ) proposed an asymmetrical variant with a reduced neighbourhood of the original Bowsher prior ( Bowsher et al., 1996 ). The Bowsher prior, in principle, is a weighted Markov random field (MRF) model which provides an edge-preserving denoising technique for enhancing PET images based on the intensities in the spatially co-registered MRI image. The weights are computed based on a similarity metric (e.g. absolute difference) evaluated on the structural image. The weights assume discrete binary values (0 or 1) based on the similarity metric. Recently, ( Schramm et al., 2018 ) evaluated the performance of the asymmetrical Bowsher prior with respect to other image gradientbased priors such as parallel level sets ( Ehrhardt et al., 2015 ) for reconstruction of static PET images. They demonstrated that the asymmetrical version yielded PET image reconstruction with improved bias-variance trade-off in comparison to the original Bowsher prior and the parallel level sets prior. Another work in ( Loeb et al., 2015 ) proposed a variant of the well-known Bowsher prior, modelled as prior information in the reconstruction process for reconstruction of dynamic PET images. All these modifications focused on improved image quality with respect to minimization of certain image quality metric such as the root mean squared error between the reconstructed PET image and the reference PET image. This work differs from the earlier works mentioned here in that we wish to estimate the underlying brain activation by improving the fPET images using anatomical information from the static/fixed MRI structural image (e.g., T1 weighted or T2 weighted MRI images).
In the current study, we hypothesized that accurate identification of brain metabolic activations could be obtained by filtering the fPET images using knowledge from the anatomical MRI image. We propose a slightly modified version of the Bowsher prior that was proposed in ( Bowsher et al., 1996 ) and further modified in ( Vunckx and Nuyts, 2010 ) within a Bayesian image denoising context. The modified Bowsher prior is a weighted quadratic MRF prior with the PET image modelled as an MRF and the weights of the quadratic penalty obtained from the MRI image. The proposed anatomical prior improves the identification of independent signal components from the fPET data by improving the PET SNR and reducing partial volume errors, as demonstrated in prior works Song et al., 2019 ;Zhu et al., 2019 ). The proposed modification corresponds to selection of the weights from the static MRI image in a model that uses patch-based similarity in the MRI image rather than a similarity based on voxel-intensity alone. Furthermore, the weights within an MRF neighbourhood are no longer binary and can have continuous values between 0 and 1. As a consequence, this work demonstrates improved identification of brain activation in response to external stimulation for both two-fold and three-fold reductions in the duration of the experiment duration to acquire the fPET data. The formulation of the prior model in this paper differs from the one proposed in ( Loeb et al., 2015 ;Vunckx and Nuyts, 2010 ), in that the current model uses a smoothly-decaying function that in turn depends on the neighbourhood voxel-locations in the static MRI image. Additionally, we incorporate patch-level differences (as opposed to voxellevel differences) to estimate the weights within the neighbourhood of a voxel. The method is henceforth referred to as MRI-MRF prior and was validated using both simulated and in vivo visual task fPET datasets. The accuracy of the method was compared with conventional smoothing methods at both the subject and group level ICA, and the in vivo fPET dynamic data were downsampled to verify the robustness of the proposed method in response to reduced task stimulation durations.

Theory
Let { } =1 represent the uncorrupted (unobserved) dynamic sequence comprising fPET images, each containing voxels. Let { } =1 represent the corresponding observed low-quality dynamic sequence comprising fPET images. In this work, we consider the sequence { } =1 as that output by the PET-MRI scanner using conventional model-based iterative methods such as maximum likelihood expectation-maximization (MLEM) ( Shepp and Vardi, 1982 ). Typically, the PET detector PSF models cross-crystal effects during the gamma-ray coincidence event, which is leveraged by the scanner during PET image reconstruction resulting in reduced partial volume errors. However, in comparison to MRI, the PET resolution is still much lower due to the fundamental limitation in the current PET detector technology. Therefore, in this work, we propose to denoise and improve the anatomical accuracy of the reconstructed PET images by relying on the MRI anatomical information modelled as a weighted MRF prior detailed below. The proposed method serves as a post-reconstruction image enhancement procedure that removes noise and reduces the partial volume errors by removing edge-artifacts.
Let represent the fixed MRI image geometrically aligned with the PET image sequence, i.e., each individual PET image is spatially aligned to the MRI image . We model the PET images as an MRF consisting of a neighbourhood system  : = {  } =1 where  denotes the neighbours of voxel . In the MRF model, we consider square-shaped neighbourhoods with the square width corresponding/equivalent to L mm. We model the prior distribution on given the fixed MRI image , as an MRF ( | ) . Because our goal is to perform image enhancement of the PET image, we consider a Gaussian likelihood model that penalizes the squared difference of the intensities in the uncorrupted (unobserved) image { } and the corrupted observed image { } . Thus, given the sequence of acquired and reconstructed fPET images, { } =1 , and the fixed MRI image of the same subject, we formulate the restoration of the fPET image sequence { } =1 as maximum-a-posteriori (MAP) estimation: Thus, the negative-log posterior with the mean-squared error data fidelity term and the MRF prior becomes Here ( ⋅) acts the MRI-guided MRF (MRI-MRF) regularization term which incorporates the anatomical information from MRI image, .
( ⋅) is the negative logarithm of the probability distribution which is the Gibbs energy associated with the MRF in terms of the weighted squared difference between the voxel intensities within the neighbourhood, where the weights are obtained from the static MRI image .The parameter determines the strength of the regularization ( ⋅) . The formulation in Eq. (2) is generic and allows incorporation of arbitrary prior models that enforce certain type of regularity, e.g. piecewise smoothness, on the fPET images. The work in ( Vunckx and Nuyts, 2010 ) proposed a modified Bowsher prior, called asymmetric Bowsher prior (Bow-Asym). The modification involved choosing a subset of cliques from the MRF neighbourhood and assigning binary (0 or 1) weights to the neighbours. For the task of PET image reconstruction, the parameters were optimized to obtain a lower root mean square error which resulted in choosing either 4 out of 18 closest voxels or 9 out of 80 closest voxels within a neighbourhood of 3 × 3 × 3 and 5 × 5 × 5, respectively. Focusing on the identification of brain activation maps, we further modify Bowsher prior model for improved estimation of the activation maps from the fPET images. Specifically, instead of involving only a few voxels within the MRF neighbourhood, we maintain a fully connected MRF with continuous weights (values in [0, 1]). In this paper, ( ⋅) is modelled as a weighted quadratic MRF function defined as, Here the weights are computed based on the intensity values from the co-registered MRI image, , as where the operator ( . ) extracts a vectorized isotropic 3D patch of desired width, centred around voxel , and the parameter determines the spatial pattern of weights within the MRF neighbourhood of voxel i . The strategy of determining the weights in the MRF-based regularization term by relying on patch-difference norms has been used within the literature on patch-based denoising methods, first proposed on natural images in the works of ( Awate and Whitaker, 2005a ;Buades et al., 2005 ) and on MRI images in the works of Whitaker, 2007 , 2005b ;Coupé et al., 2012 ). While a high value of leads to weights that are similar for all the neighbouring voxels, a low value of assigns higher weights to a few selected voxels in the neighbourhood. The latter scenario leads to an extension of the strategy in for Bow-Asym ( Vunckx and Nuyts, 2010 ) that (i) enforces neighbourhood weights to be binary (1 or 0) and (ii) designs weights only based on voxel-intensity differences (instead of patch differences). The patch size considered in the MRI image can be different from the MRF neighbourhood considered in the PET image. In this paper, we used patches of size 3 × 3 × 3 in the MRI image to compute the weights in Eq. (3) . Our proposed strategy of using patch-based differences can provide additional robustness to noise, artefacts, and minor misalignments between the MRI and PET images while leading to better structure preservation, in ways that are similar to those studied in general image denoising ( Milanfar, 2012 ). While iterative denoising algorithms, as in Whitaker, 2006 , 2005a ), offer algorithms for data-driven tuning of the parameter to improve performance, in the application in this manuscript, where the weights only need to be precomputed once (for a particular subject, due to the MRI image being static), we tune the parameter based on validation on simulated data ( Section 2.2.1 ). The Bayesian optimization problem with the MRI-MRF prior in Eq. (2) was solved using the limited memory BFGS method (L-BFGS) ( Byrd et al., 1995 ), with positivity constraints.
To perform spatial ICA, we construct a spatiotemporal data matrix, , using { } =1 , such that the dimension of is x . ICA models as a linear combination of the underlying independent components: = , where contains the independent components and is the mixing matrix. In the context of PET imaging, the measured PET data is affected by the blurring matrix, ( Bousse et al., 2012 ;Zhu et al., 2019 ), and the ICA model becomes where 0 represents the spatiotemporal matrix constructed from the true PET signals, and 0 models the true underlying independent components of 0 . The matrix, , acting on the spatial dimension, models the partial volume errors in PET measurements, and hence, the resultant independent components though the mixing operation, .
The goal of fPET data analysis is to identify 0 from Eq. (4) . Image denoising in the spatial domain is an important pre-processing step prior to application of the ICA algorithm. The characteristics of an ideal filter for estimation of the source components, 0 , would be to recover the signal without compromising the independence of the underlying true components. Typically, a Gaussian smoothing filter with a suitable width, specified by its full width at half maximum (FWHM) is used to reduce spatial noise, for example, during fMRI data analysis. However, performing a Gaussian smoothing can introduce additional bias in fPET images and the corresponding independent components, due to worsening of the partial volume errors ( Chen and Calhoun, 2018 ;Pignat et al., 2013 ). Hence, this work proposes an MRI guided filtering scheme that can (i) perform denoising, as well as (ii) reduce partial volume errors, to provide an improved estimation of the underlying source components, 0 .

Data and experiments
Both simulated and in vivo fPET and MRI data were used to validate the performance of the MRI-MRF prior. For comparison, the MRI-MRF prior processed fPET images were compared with those obtained using Gaussian smoothing with varying kernel sizes (specified by FWHM).

Experiments on simulated data
Continuous infusion of FDG PET activity was simulated for 60 minutes using a two-tissue compartment model involving the three kinetic parameters 1 , 2 , and 3 and a fitted arterial input function from intravenous blood samples collected from our previous in vivo experimental data ( Jamadar et al., 2019 ) to yield a total of 60 frames. PET image was simulated based on the T1 weighted MNI template MRI image with a voxel resolution of 8 mm 3 isotropic. The simulated FDG activity was corrected by the blood partition fraction and haematocrit using the same procedure as in our previous work ( Li et al., 2020 ). Brain regions were segmented into grey matter, white matter and the occipital cortices using the MNI structural atlas ( Mazziotta et al., 2001 ) using FSL ( Diedrichsen et al., 2009 ). The MRI and PET images were simulated with an isotropic spatial resolution of 2 mm in the MNI space. The regional specific metabolic kinetic parameters used for the simulated dataset were 1= 0.101, 2 = 0 . 071 , 3 = 0 . 042 for grey matter and 1= 0.047, 2 = 0 . 070 , 3 = 0 . 035 for white matter, respectively ( Lucignani et al., 1993 ). A visual task stimulus was simulated between 20 to 30 minutes in the visual cortex region similar to the in vivo experimental paradigm. During the visual stimulation period, the parameter 3 in the occipital cortex was simulated to have a 20% increment. The tomographic iterative GPU-based reconstruction toolbox (TI-GRE) was used for PET image reconstruction ( Biguri et al., 2016 ). The PET images were forward projected, and Poisson noise was applied in the measurement space, to generate a high-dose dataset. Subsequently, we simulated dynamic low-dose PET data using the Poisson thinning approach ( Kim et al., 2018 ) such that the low-dose data had a dose reduction factor (DRF) of 100 compared to that of the high-dose data. The generated PET images were smoothed using a 3D Gaussian filter with FWHM 2.35 mm to simulate the partial volume effect produced by the system. Finally, MLEM algorithm was used to reconstruct the PET images to yield 60 frames for both low and high dose datasets.
The ICA-specific pre-processing steps including spatial normalization and dimensionality reduction were performed as described in detail by Li et al. (2020) on the post-reconstruction smoothed images. In this work, we performed both subject-level and group-level ICA on fPET data. For group analysis, the spatiotemporal matrix from each subject was concatenated along the temporal dimension before the application of ICA ( Calhoun et al., 2004 ). The pre-processed data, which was an estimate of 0 , was then decomposed using an ICA unmixing algorithm in the FastICA toolbox ( A. Hyvarinen and E. Oja, 2000 ;Hyvärinen and Oja, 1997 ).

In vivo experiments and data
A cohort of five healthy volunteers were scanned for a visual task stimulus experiment using a 3T Siemens Biograph mMR (Siemens Healthineers, Erlangen, Germany) PET-MRI scanner, approved by the institute human ethics committee. The overall stimulation protocol consisted of three visual stimulation periods consisting of alternating periods of rest and task blocks. A detailed description of the experiment is provided in our earlier work in ( Jamadar et al., 2019 ). The subjects rested for a period of 20 minutes to allow sufficient FDG accumulation in the brain, during which structural T1 MPRAGE MRI scans were acquired. Following this, the subjects viewed a circular flickering checkerboard stimulus for 10 minutes. The checkerboard was retained for a period of 120 seconds and subsequently, an intermittent 32 seconds on and 16 seconds off design was employed. Following the first task stimulation, which involved 3 blocks: rest, task, and rest, two other stimulation experiments, using the full checkerboard visual, were carried out. We used the PET data acquired during the first full checkerboard. Hence, the PET data for each subject was of 30-minute duration, including 10 minutes resting before the stimulation, 10 minutes of a full checkboard stimulation followed by another 10 minutes of rest ( Fig. 4 (a)). The average dose of FDG given to each subject was 95.9 ± 5.9 MBq which was infused at a constant rate of 36mL/hr over the 90-minute duration.
We reconstructed PET images from the list-mode data using two different values for the temporal bin-width (Tbin) of (i) 30 seconds for the low-dose PET images, and (ii) 3 minutes for high-dose PET images. The average dose for the corresponding low dose fPET images across the group of subjects was calculated to be 7.5 kBq/kg/frame. The PET data was corrected for attenuation using a pseudo-computed tomography (pCT) map Burgos et al., 2013 ). The corrected PET data sinograms were reconstructed using ordered subsets expectation maximization (OSEM) algorithm with 3 iterations and 21 subsets along with point spread function modelling. The PET 3D volumes were reconstructed with voxel sizes of 3 × 3 × 2.03 mm 3 . For standard analysis, all the images were registered to the MNI-152 template. The high-dose PET images from the 3-minute binned data were used to register and resample the low-dose PET images with the T1 weighted MRI (acquired at 1 mm 3 isotropic resolution) for each subject using ANTS ( Avants et al., 2011 ).
Additionally, we also undertook an empirical analysis when the duration of the task and resting blocks was reduced. This analysis was carried out by downsampling the total number of low-dose fPET images reconstructed from the list mode data. Functional PET analyses were computed at both the subject-level and group-level for downsampling factors (DF) of 2 and 3 to simulate fPET images of duration 30 secs but acquired at 1:00 minute and 1:30 minute intervals, respectively. The downsampled PET images correspond to reduced task duration with a lower number of temporal frames.

Optimal kernel width selection
The optimal kernel sizes for all the methods including the optimal MRF neighbourhood for the MRI-MRF prior, for processing the in vivo data were selected and validated using simulated data. We optimized the parameters to achieve high sensitivity without substantial loss of specificity using ICA computed activation maps. For computing the sensitivity and specificity values, the region of interest (ROI), occipital cortex, was obtained using the segmentation procedure as described in Section 2.2.1 . The sensitivity and specificity performance metrics defined as follows: provide a quantitative assessment of the activation maps (z-score map) obtained from the different filtering operations. The metrics were computed by considering a voxel as activated if | | ≥ 1 . 6 and | | ≥ 2 . 1 , for subject-level and group-level analysis respectively ( Li et al., 2020 ). The parameter search-space for the MRI-MRF prior, includes varying values of the regularization parameter, MRF neighbourhood size ( , L, respectively). The parameter was determined empirically in this work. We found via simulation studies ( Section 3.1 and Supplementary material) that the sensitivity and specificity values do not change dramatically for minor perturbations in the parameter . The effect of varying the around the empirically tuned value has been demonstrated in the Supplementary material ( Supplementary Fig. 1).

Comparison methods
We evaluate our proposed method against several denoising approaches with and without anatomical information from MRI. In addition to traditional Gaussian filtering approach, because there is possibility of other edge-preserving filters for fPET processing, we evaluate our method against other denoising schemes such as 3D Bilateral filter, BM4D, and the original asymmetrical Bowsher prior (Bow-Asym) for the simulated data.
In the 3D bilateral filter, the weights for the neighbouring pixels are dependent on two Gaussian smoothing kernels, one in the spatial domain, i.e., and the other in the intensity domain, i.e., . The BM4D filter exploits the strategy similar to non-local means (NLM) filtering in identifying similar blocks of the image but differs from NLM in that it groups them to find an enhanced sparse representation in a transform domain (typically wavelet) of the grouped blocks. We used the implementation provided for ( Maggioni et al., 2013 ). Finally, we also compare against Bow-Asym ( Vunckx and Nuyts, 2010 ) which is also an MRF-based prior that uses binary weights (1 or 0) instead of the continuous weights used for the MRI-MRF prior. Table 1 compares the sensitivity and specificity for both denoising schemes at different parameter configurations. For the MRI-MRF prior, the MRF neighbourhood was varied based on cube-shaped patches of width from 10 mm to 18 mm which represented a width of 5 to 9 voxels in each direction, respectively. In the case of Gaussian filtering, the kernel size was determined by the full width at half maximum of the Gaussian function. The FWHM for Gaussian varied from 11 mm to 15 mm. It is to be noted that while the parameter L represents the width of the entire patch in the PET image (for the MRI-MRF prior), FWHM (for the Gaussian filter) represents approximately half of the kernel-width. The parameter range chosen for the Gaussian smoothing is consistent with the Gaussian kernel widths used in the prior work ( Li et al. 2020 ). The sensitivity values for the MRI-MRF processed image are dramatically higher than that of the Gaussian smoothed images, whereas the specificity values are comparable between the two methods. For the fPET data analysis, an MRF neighbourhood corresponding to a length of 14 mm (in all directions) was chosen for the MRI-MRF prior. However, both the Gaussian kernels with FWHM 11 mm and 13 mm show similar sensitivity and specificity values. Therefore, the analysis using the Gaussian-filtered in vivo fPET data was undertaken using both the 11 mm and 13 mm FWHM filters. Table 2 shows the effect of varying the free parameter for a given fixed on the sensitivity and specificity values for the visual-task component obtained from the ICA. The sensitivity and specificity values for varying prior-weight is shown for two different MRF neighbourhood L = 10 mm and 14 mm. The sensitivity and specificity values in Table 2 correspond to a z-score threshold value of 1.6 (i.e., | | ≥ 1 . 6 ).

Results for simulated data
We tune the parameters for the other denoising schemes following similar strategy as above, i.e., obtaining improved activation from ICA, The optimal values for bilateral filter were determined to be = 3 pixels (6 mm) and = 0 . 6 max , where max is the maximum intensity in the noisy PET image. For the BM4D filter, we tuned the noise standard deviation parameter BM 4D = 0 . 9 max . For the Bow-Asym, a slightly larger set of cliques provided an improved estimate of activation compared to using the standard configuration of 18 nearest voxels out of a 5 × 5 × 5 neighbourhood. For Bow-Asym, we used 50 voxels out of 125 voxels (5 × 5 × 5). Fig. 1 shows the denoised images obtained using different methods, for visual comparison. The difference between the 3D bilateral filter and the BM4D filter is negligible. However, both the 3D bilateral filter and the BM4D filter perform substantially better than the Gaussian filter in terms of preserving overall structural regularity without oversmoothing. The MRI-MRF method (fourth column) recovers edges and textures that are closer to the reference PET image (last column). Similarly, the Bow-Asym is very close to the image quality shown by MRI-MRF, with both of them showing substantially improved contrast compared to other methods. Although the anatomical delineations in both bilateral and BM4D filters are not as clear as the MRI-MRF method, they do show better contrast between different structures compared to the Gaussian filter. Fig. 2 shows the visual task specific activation for the reference noiseless fPET images, and for the five denoising schemes using the optimal parameters. The ICA activation map obtained using the noiseless images serves as the reference map ( Fig. 2 (a)). The activation map obtained from post-reconstruction filtered fPET images using the MRI-MRF prior ( Fig. 2 b) was closest to the reference activation map in the visual cortex. The corresponding activation maps using Bow-Asym ( Fig. 2 c) was able to recover most of the activation but resulted in reduced sensitivity (82.7) compared to the MRI-MRF prior (87.7). On the other hand, the activation maps obtained using Gaussian smoothing ( Fig. 2 d) yields suboptimal activation maps in the visual cortex with asymmetrical patterns. The maximum sensitivity obtained using the 3D bilateral ( Fig. 2 e) and   BM4D filter ( Fig. 2 f) were found to be similar to that of the Gaussian filters. The maximum sensitivity obtained by tuning the parameters for the 3D bilateral filter and the BM3D filters were 0.78 and 0.76 respectively for | | ≥ 1 . 6 .

Results for in vivo data
The results reported in this section are for the fPET images reconstructed using the list-mode data binned at Tbin = 30 s, and for DF = 1, 2 and 3. For the in vivo data, we demonstrate the efficacy of our method in comparison to Gaussian filters with varying FWHM values due to its popularity in functional neuroimaging processing. Fig. 3 shows the post-reconstruction filtered fPET images along with the subject's MRI image ( Fig. 3 , left column) and the corresponding vendor provided low dose fPET image ( Fig. 3 , second column). The denoised image using the MRI-MRF prior ( Fig. 3 , third column) shows superior recovery of PET signals in different regions of the brain while removing substantial amount of noise. Specifically, the white and grey matter boundaries are well delineated, the shape of the ventricles has been recovered (which is not evident in the low dose PET image), and anatomical features in the gyri, sulci and details of the cortical folding (refer Fig. 3 m) have been restored. On the other hand, the denoised images using both Gaussian kernels (FWHM 11 mm and 13 mm) are heavily blurred and show substantial loss of anatomical details due to the partial volume errors ( Fig. 3 , fourth and fifth columns). Fig. 4 shows the results of an individual subject-level fPET analysis obtained using different filtering techniques for a downsampling factor of one (i.e. DF = 1, includes all list-mode data). The ICA activation maps corresponding to the visual task component along with the timecourses are calculated for each method. The timecourses are derived from the normalized (z-scored) spatial components following the methods in Li et al. 2020 . The component maps for all sections of the brain are are noisy and do not closely follow the experimental task paradigm. Moreover, the shape of the region of brain activation does not follow the known anatomical structure of the primary visual cortex but extends into adjacent neuroanatomical structures including the white matter, likely due to large partial volume errors. Conversely, the activation map obtained using the MRI-MRF prior ( Fig. 4 c) shows localized activity near the visual cortex with a significantly higher z-score within the visual cortex compared to both Gaussian kernels. The ICA timecourse for the MRI-MRF prior ( Fig. 4 b) accords more closely with the experimental design with increased uptake during the visual task block. The comparison of the visual task components for the three methods for all brain sections is consistent with these observations (see Supplementary material).
The results for the group-level fPET analyses for the three filtering techniques using the complete list-mode dataset (DF = 1) are shown in Fig. 5 . A higher z-score range was observed for the group-level analyses compared to the single subject-level analysis. However, in contrast to the subject-level analysis, the timecourses estimated from all methods ( Figs. 5 b, 5 d & 5 f) at the group level recapitulated the experimental design paradigm. The activation map corresponding to the MRI-MRF prior followed the known neuroanatomical representation of the primary visual cortex and was consistent with the subject-level result. On the other hand, the activation maps using the two Gaussian kernels did not represent activation in the primary visual cortex and demonstrated diffuse cerebral metabolic activity into large adjacent anatomical regions including the white matter. Fig. 6 shows the subject-level fPET analyses for downsampling factors of two and three (DF = 2 & 3). Timepoints T1 and T2 represent the onset of the task and second resting block in the downsampled task paradigm. Plausible ICA activation maps were not generated using an 11 mm FWHM Gaussian kernel for both DFs and therefore no results are included. The ICA timecourses during the task-block for the MRI-MRF filter demonstrated a steadier gradual increase, in agreement with the task paradigm, in comparison to the 13 mm FWHM Gaussian kernel for DFs of 2 and 3 respectively ( Figs. 6 a & 6 e compared to Figs. 6 c & 6 g respectively). The activation map axial view for DF = 3 did not reveal activation in the left hemisphere as was expected for the visual task ( Fig. 6 h). However, for DF = 2 there was some activation in the left hemisphere visual cortex ( Fig. 6 d) although it was not as widespread as for the fully sampled dataset. On the other hand, the activation maps for the MRI-MRF prior ( Figs. 6 b & 6 f) showed spatial congruency across the three DFs, whilst the discrepancy between the z-scores for the MRI-MRF prior and the 13 mm FWHM Gaussian filter was largest for DF = 3 compared to DF = 2 and 1. Fig. 7 shows the group-level fPET analyses at DF = 2 and DF = 3. In contrast to the group-level analysis for the fully sampled dataset where there was little difference between the activation maps estimated by the MRI-MRF method and the Gaussian kernel with FWHM 13 mm ( Fig. 5 ), the activation maps estimated for the group-level analyses for DF = 2 and DF = 3 showed marked differences. For both DF = 2 and 3, the ICA timecourses for the MRI-MRF prior ( Figs. 7 a & 7 e) showed agreement with the task experimental design with higher z-scores than for the 13mm FWHM Gaussian filter timecourses ( Figs. 7 c & 7 g). The activation maps show that while the MRI-MRF prior was able to resolve brain activation that was consistent with activation of the visual cortex ( Figs. 7 b & 7 f), at both DF = 2 and 3 the 13mm FWHM Gaussian filter was unable to resolve extended activation throughout the primary visual cortex ( Figs. 7 d & 7 h) with no activation in the left hemisphere for DF = 3 ( Fig. 7 h). Conversely, the activation maps for the MRI-MRF prior were congruent across the subject-level and group-level analyses although greater consistency in the right hemisphere.

Discussion
We have proposed an MRI-assisted fPET processing framework for the analysis of task-related metabolic changes in the brain using high temporal resolution fPET data and for low-dose fPET brain mapping applications. We investigated the effect of using the anatomical information from a subject's MRI to denoise the fPET dataset to reduce the partial volume error in the PET images in order to increase the sensitivity of the ICA analysis. The PET image restoration problem was posed as a solution to a Bayesian optimization problem which was solved using L-BFGS due to its greater computational efficiency compared to gradientdescent based optimization techniques.
This study compared the post-reconstruction filtered images from the MRI-MRF method and Gaussian smoothing with varying kernel sizes as well as the ICA activation maps from the fPET dataset using a visual stimulation task. Visual assessment of the post-reconstruction smoothed images showed that the MRI-MRF processed PET images recovered many features which were not readily observed in the conventional low dose PET images. The MRI-MRF filtered PET images revealed localized tracer uptake in the sub-cortical nuclei adjacent to the lateral ventricles (e.g., Fig. 3 c) whereas little or no uptake was apparent in the comparable low-dose and Gaussian denoised PET images. Furthermore, the level of Gaussian smoothing required to obtain plausible activations in the visual cortex rendered the fPET image hard to interpret visually as there was a substantial loss of features. The MRI-MRF method provides a balance between visual interpretability of the fPET images together with improved resolution and sensitivity for functional analysis using ICA.
The task-based experimental design paradigm enabled meaningful comparison of the ICA timecourses obtained using the two filtering techniques, by inspection of the FDG uptake in the visual cortex during the visual stimulation task. The proposed methodology was able to achieve consistent activation maps at both the subject-level and group-level for DF = 1, 2 and 3. However, this was not true for the Gaussian smoothing kernels. Moreover, since the fPET data was acquired for an exogenous task-based stimulus, good correlation between the subject-level and group-level activation maps was expected. In particular, the improved results for the individual subject-level analysis demonstrates the benefit of the MRI-MRF method to enhance single subject-level analysis using low dose high temporal resolution fPET data with reduced task durations.
The study involving downsampling factors that demonstrated that the proposed processing pipeline could detect dynamic brain metabolic increases for visual task stimulation for as short as approximately three minutes. However, this interpretation does assume that the FDG dosage per frame in the fPET images is consistent for the different downsampling factors. In practice this would be achievable experimentally by altering the infusion protocol or slightly increasing the dosage of the radiotracer ( Verger and Guedj, 2018 ). The Gaussian smoothing technique failed to identify task related ICA components for the shorter task durations (i.e. at higher DFs) due to reduced sensitivity.
Unlike fMRI, fPET images suffer from very low SNR and hence the spatial denoising scheme must be carefully chosen to provide an optimal bias-variance trade-off. MRI-guided PET image denoising and deblurring has been extensively reported in the literature Song et al., 2019 ) with many solutions for post-reconstruction PET image enhancement. However, this paper is the first to demonstrate the effectiveness of the MRI-based spatial denoising technique for dynamic fPET imaging, such that fPET images are both visually interpretable and produce accurate functional maps with improved temporal resolution. The high specificity and sensitivity of the algorithm also enabled single subject-level analyses along with reasonable visualization of the fPET images without loss of anatomical details. Traditional methods such as Gaussian smoothing perform averaging without consideration of the anatomical boundaries and hence the quantitative accuracy of FDG signals is degraded by partial volume errors. This was reflected in the diffuse visual activation maps obtained with the Gaussian filtering. Using edge-preserving denoising techniques such as anisotropic filtering would also yield suboptimal performance because of the poor SNR of the fPET images and the difficulty to distinguish between tissue boundaries and noise. Fig. 6. Subject-level (representative) estimation of brain activations using the reduced task and resting blocks with DF = 2 and 3 at MNI coordinate (14, -94, -8). The independent components estimated from the filtered fPET images using different schemes are provided. DF = 2: MRI-MRF prior (b) and Gaussian kernel with FWHM 13 mm (d). DF = 3: MRI-MRF prior (f), and Gaussian smoothing with FWHM = 13 mm (h). The T1 and T2 represent the onsets of the task and second resting block, respectively.
The formulation of the MRI-MRF prior in this work is generic and allows for modelling of higher-level image features such as dictionary atoms. Nevertheless, the proposed filtering framework is efficient and computationally less expensive in comparison to other patch-based techniques and hence, the framework is easier to adapt for other research and clinical applications.
Although the MRI-MRF prior in this work was applied in the image domain on post-reconstructed fPET images, the same could be applied within the image reconstruction process provided the PET list-mode data was accessible. Research using image restoration techniques in a reconstruction framework generally employ a Poisson noise model for the sinogram data and a system matrix composed of several matrix operations representing the point spread function, forward projection, attenuation correction, scatter correction, and back projection. Our work solved the image restoration problem in the image domain and employed a least-squares-type data term rather than a fixed noise-model in the image space. This was because the noise characteristics of the reconstructed PET images inherently depend on the reconstruction algorithm. For example, noise characteristics during filtered back projection reconstruction depend upon the filter employed, whilst in maximum likelihood expectation maximisation reconstruction and its variants, the noise characteristics depend on the number of iterations as well as the strength of the prior function.
In this work, we have applied an MR assisted PET data processing method to identify brain metabolic activations using FDG PET. This concept can further be applied to other types of PET experiments such as neuroreceptor ( Ametamey and Schubiger, 2006 ) and neuroinflammation tracer studies ( Alauddin, 2012 ), where counts are statistically limited during dynamic imaging for kinetic modelling analyses. However, these tracers will have different binding targets and biodistribution in brain tissues compared with FDG PET. Therefore, it is important to consider appropriate MR contrasts when designing the MRI-MRF prior model. For example, ( Kaunzner et al., 2019 ) investigated the correspondence of quantitative susceptibility mapping (QSM) MRI and [11C]-PK11195 PET for the localization of tracer associated with neuroinflammation.
The current work has a number of limitations. One of the limitations is the small sample size. In this work, we show proof of the principle for utilizing anatomical information for fPET data processing. Advanced statistical image restoration models such as joint patch-based techniques and neural networks may further improve the image quality for shorter image acquisition durations and potentially in future approach the temporal resolution offered by fMRI. However, the proposed framework is readily adaptable to use these techniques in the research context although modelling higher statistical dependencies would increase the number of hyperparameters that were required to be tuned. A further limitation is that the MRI-MRF modelled as a Bowsher-like prior may be perceived as a technique that relies excessively on the anatomical modality. Although this may be relatively unimportant or in fact beneficial in the case of tracers like FDG that are widely distributed Fig. 7. Group-level estimation of brain activations using the reduced task and resting blocks with DF = 2 and 3 at MNI coordinate (14, -94, -8). The ICA components estimated from the filtered fPET images using different schemes are provided. DF = 2: MRI-MRF prior (b) and Gaussian kernel with FWHM 13 mm (d). DF = 3: MRI-MRF prior (f), and Gaussian smoothing with FWHM = 13 mm (h).
throughout the brain, this may not be the case for other heterogeneously distributed tracers such as for amyloid PET imaging. More sophisticated image restoration models which maintain a balance between the PET and MRI features for each tracer may need to be incorporated at the cost of more computational time. The use of spatial regularization could be carefully extended to include a temporal smoothing constraint governed by studies in tracer kinetics. A comprehensive study of several MRI-PET joint priors in the context of dynamic functional PET denoising and analytical techniques is an important direction for future studies.

Conclusion
We have presented a novel MRI-assisted fPET processing framework for functional analysis of fPET data at high temporal resolution and for low doses of radiotracer. Compared to traditional Gaussian smoothing, our approach yields visually interpretable PET images while increasing the sensitivity and anatomical accuracy of activation maps estimated using ICA. Through validation using simulated data, we have demonstrated that the MRI-MRF method is able to accurately estimate visual task related brain activation maps even under poor SNR conditions. The application to in vivo fPET data demonstrated that the MRI-MRF prior achieves detection of reduced task durations of approximately three minutes and provides an avenue for further increases in the temporal resolution and sensitivity of both single subject and group-level brain metabolic mapping studies.

Data statement
Data and code from this study are available on request from the corresponding author pending the institute Ethics approval.

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi:10.1016/j.neuroimage.2021.117928 .