MRI-guided PET Reconstruction with Adaptive Prior Strength: Study on Image Quality at Different Levels of Statistics


 BackgroundThe combination of magnetic resonance imaging (MRI)-anatomical information with positron emission tomography (PET) image reconstruction has been shown to improve PET image quality in terms of spatial resolution and image noise, especially in brain PET imaging. There are different approaches to combine MRI and PET available, being the use of a Bayesian framework the most extended. Generally, the strength of the prior is controlled by a hyperparameter that needs to be tuned depending on the acquired statistics/counts and the desired image quality in the resulting PET image. However, comparisons between methods is scant, and it is not clear how sensitive they are to the different levels of statistics that can be measured in a PET scan. MethodsIn this work we employed maximum a posteriori (MAP) reconstruction with MRI information to guide the PET reconstruction, and evaluated the performance of several prior models and optimization methods with a fixed and adaptive hyperparameter. Different simulated scenarios, and measured data using different radiotracers at different levels of statistics were employed for the evaluation. Comparisons in image quality and quanti cation between methods were performed in different brain cortical and sub-cortical regions. ResultsSimulated data showed that an adaptive hyperparameter consistently outperformed a fixed hyperparameter for every image reconstruction algorithm implemented. The best performance was achieved with a model combining the Bowsher prior weighted by similarity coefficients based on joint entropy between PET and MRI. One-step-late (OSL) and preconditioned gradient ascent (PGA) optimization methods performed similarly at any level of statistics and number of iterations, so long as the hyperparameter was adaptive. ConclusionsResults with simulated and measured data agreed that MAP reconstruction outperformed OSEM reconstruction, especially at low level of statistics, without any need of tuning. High resolution and low noise images were obtained using MAP reconstruction for 5{30 minutes scan times, showing negligible image quality difference for different radiotracers.


Background
Positron emission tomography (PET) image quality is mainly degraded by two factors: statistical noise and spatial resolution, mainly due to the low amount of radiotracer injected in subjects, limited scan time, and scintillator crystal size.
Using statistical iterative reconstruction methods, as more iterations are computed, the reconstructed image converges to the most likely result, at the cost of high frequency image noise [1]. A reconstructed image at a high number of iterations has the best achievable spatial resolution (limited by the scanner physical properties) and high quantification accuracy (as long as the optimizer does not introduce bias). However, the noise present in the data gets amplified and propagates to the image as more iterations are calculated, resulting in high frequency noise in the image, hindering its use for any imaging purpose.
In practice, few iterations are usually calculated followed by a low-pass filter, trying to reach a trade-off between noise, spatial resolution, quantification accuracy and reconstruction times acceptable in clinical practice.
PET image reconstruction is an ill-posed mathematical inverse problem, where including additional information to the measured data, using Bayesian maximum a posteriori (MAP) image reconstruction, is a way to regularize the reconstruction process [2,3].
The combination of magnetic resonance imaging (MRI), either acquired simultaneously or sequentially (shortly after) with PET, may potentially improve the spatial resolution and reduce image noise by guiding the PET image reconstruction. This is especially true in brain imaging, where MRI usually produces high isotropic spatial resolution, low noise and high contrast imagery.
The introduction of simultaneous PET/MRI scanners facilitated this approach since the PET and MRI scans can be taken at the same time point, providing an optimum match between anatomy and functionality, as opposed to scans taken at different time points where the disease can have evolved between both scans. T1and T2-weighted 3D MRI sequences are usually employed in brain imaging to obtain the best spatial resolution and noise properties. These sequences take between 2-5 minutes to be acquired, while PET scans are generally longer.
The combination of PET and MRI reconstruction has been previously explored [4,5,6,7,8,9,10]. However, there are no published works addressing clinical benefits or showing statistically differential outcomes between healthy volunteers and patients, from these improved PET images.
Two of the premises of combining MRI to guide the PET reconstruction are that if a set of voxels have similar intensity values in the anatomical image, they should have similar intensity values in the functional image, and that neighboring voxels have similar intensity values. A widely used model to describe a prior that fulfills this behavior is the Markov random field (MRF). However, there can be cases where these two premises may not be fulfilled, as is the case of heterogeneous tumors [11].
To combine the MRI information in the reconstruction process, a segmentationfree approach is generally desirable, to avoid segmentation errors propagating into the PET image. Following the above mentioned premises, Bowsher et al. defined a similarity coefficient constrained to those voxels that are similar in the anatomical image [12]. This method has been used and extended in several of the MRI-guided PET reconstruction approaches developed to date [6,7,8,10].
Usually, several MRI sequences are acquired in any PET/MRI scan to study the different properties of specific tissue characteristics. This results in observing different tissue features through different image contrasts (different MRI images).
Mehranian et al. [7] introduced the possibility of combining more than one MRI image in the MAP reconstruction framework, showing with simulated data that if a lesion is absent in one MRI sequence but present on a different MRI sequence, the lesion can be well recovered in the PET reconstructed image. This approach was termed as multiparametric PET/MRI reconstruction, given the multiple MRI sequences employed here.
One typical problem that most methods face is how to balance the influence of the anatomical information over the measured information. This critical part of the algorithm is generally done with a hyperparameter (β) which requires to be tuned depending on the amount of measured data (net true counts) and to control the smoothness and image sharpness of the reconstructed image.
Few works have been devoted to minimize or eliminate the tuning process from the reconstruction process, making the hyperparameter value adaptive during the reconstruction depending on some measure in the reconstructed image. Making the hyperparameter value dependent on the counts in each voxel at each iteration has been explored [13], in an effort to eliminate false positives in PET/CT MAP reconstruction. A more recent approach proposed to use intermediate regularized bootstrap replicates of the image estimate to automatically tune the hyperparameter value at each iteration [14].
In this work different reconstruction approaches were implemented to assess the impact of a fixed and an adaptive hyperparameter on the image quality at different levels of statistics. Different optimization algorithms, potential functions and similarity coefficients were explored. The evaluation consisted of a combination of Monte Carlo simulations based on the Biograph Vision scanner (Siemens Healthineers), and measured data acquired with the Biograph mMR scanner (Siemens Healthineers). Monte Carlo simulations of brain PET scans, describing realistic and challenging scenarios were explored. Image quality was further evaluated in measured datasets obtained with different radiotracers: [ 18 F]-Fluorodeoxyglucose (FDG) and [ 18 F]-Florbetapir.

Reconstruction Methods
Being y the set of measurements and u the measured activity distribution, MAP reconstruction aims at maximizing the logarithm posterior probability of û u = arg max{log p(y|u) + log p(u)}, where p(u) can be modeled as a Gibbs distribution exp −βR(u) , leading tô being R(u) the prior model, and β the hyperparameter that controls the strength of the additional information over the acquired data. R(u) can be expressed in a general form as an MRF [7]: where j is the voxel index, b is the index for the voxels in the neighborhood N j , N is the number of voxels, φ is a potential function operating in N j , ξ is the proximity coefficient, ω is the similarity coefficient, and ψ is the potential function that defines the relation between each voxel and its neighbors.
We explored different potential functions, weights on the potential function, also called similarity coefficient, which define the strength of each two pairs of voxels, and different optimization algorithms.

Prior Model: Potential Function and Similarity Coefficient
With respect to R(u), in this work we investigated two potential functions (ψ(u) eq. (3)): the quadratic difference The PFQ is known for smoothing edges when high β values are used. That undesirable effect has motivated the development of alternative potential functions, such as the RD. However, there are no comparisons between PFQ and RD using anatomical information in the prior model.
The similarity coefficient (ω in eq. (3)) weighs the potential between two voxels based on their intensity similarities. Bowsher et al. introduced the concept of determining this similarity based on an underlying anatomical image [12], setting ω jb to 1 if those voxels were similar in the MRI images, and 0 otherwise. In its original version, how similar two voxels are, is determined by including a fixed number of voxel-pairs, those with the lowest intensity difference in the anatomical image, hence several studies have implemented the Bowsher method [15,7,8]. The method does not require segmentation of the anatomical image. However, a region of interest (ROI) size and a fixed number of voxel-pairs were pre-determined.
The approach used in this work is a modification from the formulation introduced by Kazantsev et al. [16], where an intensity difference histogram of u j and each voxel in N j is calculated, and only those voxels u b that correspond to the first bin of the histogram are considered in the calculation of the prior. This approach prevents from fixing the number of voxels in N j , making it adaptive to edges and uniform regions. The initial size of the ROI, where the intensity difference histogram is calculated from, is fixed to 7×7×7 voxels.
In the case of the similarity coefficient, the work by Mehranian et al. [7] explored a number of coefficients, concluding that the similarity coefficient based on Bowsher and the Burg joint entropy outperformed the rest of the investigated coefficients, especially when recovering lesions present in the PET data but not in the anatomical data. Additionally, the model developed in [7], allowed the inclusion of more than one MRI sequence in the calculation of the similarity coefficient. This similarity coefficient is defined as where v is the MRI image, k is the index of the MRI sequence from the multiparametric approach, σ is the standard deviation in N j in u and v , and p(u j , v j ) is defined as a non-parametric Parzen window using Gaussian kernels: ). (5) In this work we explored the traditional Bowsher weight (Bw) and its extension including the Burg joint entropy weight (JEw).

Optimization Algorithm
The optimization algorithm employed to maximize the cost function, should in principle not have an impact of the resulting image. However, in practice, different optimization algorithms converge at different speeds, hence potentially producing different images when stopped at a given iteration number. Furthermore, some optimization algorithms may run into numerical problems and instabilities.
In this work, the one-step-late (OSL) algorithm [17], which has been used in some of the works described in this study [7,4], was first implemented. OSL is known for having numerical problems when the β value is so high, that the prior dominates over the sensitivity, resulting in unstable or negative values in the reconstructed image. To avoid negative values, β is set to 0 in this work (effectively turning into OSEM reconstruction) when the prior dominates over the sensitivity in a per-voxel basis. OSL is defined as i being the line-of-response (LOR) index, M the number of LORs, r i is scatter and random coincidences, n i is the normalization, a i is the attenuation factors, and g ij is the system matrix.
MLEM has been re-writen by other authors as a preconditioned gradient ascent (PGA) manner [3], to avoid the instabilities above mentioned with OSL. This approach has been used and modified by other authors [8,15], and is defined as when ignoring the second derivative of R(u) from the original formulation [3]. Finally, we also implemented the method introduced by Wang et al. [18], which is based on constructing the penalized likelihood based on a separable surrogate function (PLSS) based on De Pierro's decoupling rule [19], where the update formula is where u n+1 j,EM is the EM estimate of u n j , followed by a smoothing step defined as This optimization method was used by one of the scant MRI-guided PET reconstruction methods published do far, focused on adaptively calculating the hyperparameter between iterations [14].
Inspired by the work in Salomon et al [13], we also explored the impact of an adaptive hyperparameter β value on the image quality at different levels of statistics. In [13], the β value is re-calculated at each iteration based on the noise in each image voxel, described as follows: where α is a constant. All these different reconstruction alternatives were implemented based on the original reconstruction code from the scanner, which uses a 2.09 mm voxel size for the Biograph mMR and a 1.65 mm for the Biograph Vision. Prior to the MAP reconstruction, the MRI sequence was registered to a PET image reconstructed with OSEM using SPM12 [1] . The registered MRI image was then used to guide the PET reconstruction.

Simulated Data
The GATE Monte Carlo package was used to simulate one of the brain subjects from the brainweb database [2] . This database contains 20 models of normal brains divided in different tissues [20]. By downloading the individual tissues from a subject, one can compose different images, like update distribution, attenuation map and different MRI contrasts. A normal subject was downloaded per tissues with a 1 mm isotropic resolution, so the activity image, attenuation map, and T1w and T2w MRI images were later composed. Figure 1a shows the ground truth of a simulated [ 18 F]-fluodeoxyglucose (FDG) scan, together with the T1w and T2w MRI sequences. [1] https://www.fil.ion.ucl.ac.uk/spm/software/spm12/ [2] https://brainweb.bic.mni.mcgill.ca/ The scanner model that we used for the simulations was the Biograph Vision 600. The singles from the output root file were processed to generate a listmode file using the PETLINK format [3] , compatible with the rebinner and the reconstruction tools from the vendor. Coincidences and random events were calculated the same way as the coincidence processor in the scanner does.
The activity in the simulated brain was 22.6 MBq, resulting in 152×10 6 net true counts in a 5 minutes scan, with a 7.2% random fraction (RF) and a scatter fraction (SF) of 29.0%. The listmode file was additionally reframed to 1 minute (30×10 6 net true counts). Figure 1b shows a measured brain in the Biograph Vision for 5 minutes (136 MBq injected -101×10 6 net true counts, 32% SF, 29% RF), and a simulated brain for 5 minutes, both reconstructed using OSEM with TOF and PSF for 5 iterations and 5 subsets, voxel size 1.65×1.65×1.65 mm 3 in a 440×440×159 matrix. It is worth mentioning that the simulated brain did not contain the bed in the field of view (FOV), and there was no activity outside of the FOV, which explains the differences in RF and SF. Figure 1: Simulated FDG scan (top row), MRI-T1w (middle row) and MRI-T2w (bottom row) from brainweb database (a). PET brain scan measured in the Biograph Vision (top) and reconstructed simulated brain scan for 5 minutes from the brainweb database (bottom) using the vendor reconstruction and correction tools (b).
The same brain phantom was simulated adding three lesions, two in the white matter and one in the gray matter. The lesions in the white matter had three times the activity in the background, one had a diameter of 11.6 mm (0.8 mL volume) and the other one 5.0 mm (0.07 mL volume). The lesion in the gray matter had four times the activity in the background and a diameter of 5.0 mm (0.07 mL volume). There are two versions of the T1w image which was used to guide the PET reconstruction: one containing the lesions and another one not containing any lesion.
The same subject was also simulated by introducing an artificial reduction of 20% metabolism in the frontal cortex, to further evaluate the impact of MAP reconstruction in quantification, and compare with the regular OSEM reconstruction. [3] https://www.siemens-healthineers.com/en-us/molecular-imaging/petlinkdocuments

Measured Data
We used different brain scans acquired in the Biograph mMR, to test and evaluate the robustness and impact on quantification of the different methods presented here. The measured data were reframed in order to test the robustness of the reconstruction algorithms against different levels of statistics.
First, we used a static FDG scan of a diagnosed Epilepsy patient. FDG is used in Epilepsy to locate the epileptogenic region, which tends to show hypometabolism as compared to healthy tissue. The patient was injected 200 MBq and scanned for 30 minutes, 42 minutes post-injection, collecting 535×10 6 net true counts. The listmode file was reframed to 5 minutes and 1 minute, corresponding to 96×10 6 and 19×10 6 net true counts respectively.
In addition, we used two scans with [ 18 F]-Florbetapir, used to study the β-amyloid burden in Alzheimer's disease. [ 18 F]-Florbetapir shows non-specific uptake in the white matter, and low uptake in the gray matter in a healthy brain, showing high contrast between gray and white matter. However, a brain containing β-amyloid plaques will show high [ 18 F]-Florbetapir uptake in the gray matter also, reducing the contrast between gray and white matter.
The [ 18 F]-Florbetapir patients were injected 377.4 MBq and 373.7 MBq, and were scanned for 20 minutes, 50 minutes post-injection. The listmode files were reframed from the original 20 minutes, down to 5 minutes and 1 min.
All reconstructions were performed using OSEM with 3 iterations 21 subsets, 0.8×0.8×2.0 mm 3 voxel size (zoom 2.5) in a 344×344×127 matrix size. PSF modeling was used in all the reconstructions shown in this work.

Image Analysis
To evaluate the image quality obtained in the gray matter, in the case of the simulated and measured PET scans, we performed a volume of interest (VOI) analysis using SPM12, and calculated different figures of merit (FoM) in several cortical and sub-cortical regions. The following steps describe the VOI-analysis: • MRI volume is rigidly registered to PET volume.
• Registered MRI volume in PET space is normalized to the Hammers atlas [4] (1.0×1.0×1.0 mm 3 voxel size). • Normalization parameters from previous step are applied to PET volume, obtaining the PET volume warped onto the Hammers atlas. • Normalized MRI volume is segmented into gray matter, white matter, and CSF. • The regions from the Hammers atlas are masked by the gray and white matter. Only voxels with above 90% probability of belonging to the region were considered. Several sub-cortical and cortical regions were selected for the VOI analysis. We focused the analysis on multiple regions for the FDG datasets, trying to cover a spectrum from small to large regions, including hot, warm and cold regions. In the case of [ 18 F]-Florbetapir datasets, we focused the analysis on those regions known to present high burden of β-amyloid plaques.

Simulated Data Image Analysis
Since one of the objectives in this work is to find an MRI-guided PET reconstruction method that is robust to datasets with different levels of statistics, one way the methods were compared was by measuring the difference between the 5 minutes and 1 minute scans. The analysis was focused on the hippocampus, amygdala, putamen, accumbens, caudate, the frontal, occipital, parietal and temporal lobes, the cerebellum, and the ventricles. The difference between the 5 minutes and 1 minute scans was measured through the mean relative difference in bias (∆b) for all the brain regions, defined as where N ROIs is the number of brain regions, r is the region index, the subscripts 5m and 1m refer to the 5 minutes and 1 minute scans, µ r,x is the mean intensity value in the brain region r for the 5 minutes or 1 minute scans, and µ r,ref is the mean intensity value in the brain region of the ground truth used for the simulation (fig. 1a).
We further compared those methods that did not show important differences (<10%) between the 5 minutes and 1 minute scans with OSEM, by looking at the bias versus noise from 1-30 iterations (21 subsets) in all the cortical regions.
Regarding the simulated brain with lesions, we focused the quantitative analysis on the largest lesion (11.6 mm diameter -0.8 mL). The quality to recover and detect the lesion was evaluated by automatically calculating the lesion size and the ratio between the activity in the lesion and a ring surrounding the lesion, as a measure of lesion contrast/detectability.
The lesion size was obtained by segmenting those connected voxels with intensity values greater than half the maximum intensity value in the lesion.
The activity in the lesion and a ring around the lesion were obtained by applying a template on the lesion with the known lesion size and shape, and a three voxel ring around the lesion respectively. The simulated lesion-to-outer-contour ratio was 3. However, the phantom used for the simulations has a 1.0 mm 3 voxel size, while the image analysis was performed in the reconstructed voxel size: 1.65×1.65×1.65 mm 3 . This loss of resolution resulted in a degradation of the lesion-to-outer-contour ratio of 2.25, which was considered as the reference of the analysis.
The simulated brain with the reduced activity concentration in the frontal cortex was analysed by filtering each reconstructed image with the mask of the region where the activity was reduced, and comparing the resulting mean activity concentration reduction in that region between reconstruction methods.

Measured Data Image Analysis
The method used to obtain the VOIs in the FDG measured data was the same as for simulated data. The analysis was focused on the caudate (small region), frontal lobe (large region), cerebellum (reference region in semi-quantitative analysis), and ventricle (cold region). We measured the mean intensity and the noise, measured as coefficient of variation (CoV), for each region for 1-10 iterations.
We compared the OSEM algorithm only with those methods that performed the best with the simulated data. The comparison was done with the original 30 minutes dataset (usual clinical scan time), and the decimated 5 minutes and 1 minute scans.
The clinical [ 18 F]-Florbetapir datasets were analysed in the same way as the FDG data, but the focus was on a different set of regions. Studies have shown that the parietal, temporal and frontal lobes, and posterior and anterior cyngulate gyrus, are regions that can potentially show β-amyloid plaques, with the cerebellum used as reference [21].
The clinical scan time was 20 minutes. We decimated the datasets to 5 minutes and 1 min. Taking the 20 minutes scan as reference, we calculated the relative difference for each region in both patients for the 5 minutes and 1 minute datasets. Figure 2 shows a central axial slice of the 5 minutes and 1 minute simulated brain dataset, reconstructed with all the explored reconstruction algorithm variations considered in this work: optimization algorithms (OSL, PGA, PLSS), potential functions (PFQ and RD), and similarity coefficients (JEw and Bw), for fixed (HPf) and adaptive (HPa) hyperparameter value.

Simulated Data
The hyperparameter value was set to obtain optimum image quality for the 5 minutes scan, and it was kept constant for the 1 minute scan. The hyperparameter value used in each case was 10 2 for PFQ-JEw, 0.1 for PFQ-Bw, 5.0 for RD-JEw, 0.005 for RD-Bw (same values for OSL and PGA), 10 5 for PLSS-JEw and 10 2 for PLSS-Bw. In the case of the adaptive hyperparameter value, α (eq. 10) was set to 10 2 to obtain equivalent image quality for the 5 minutes scan between HPf and HPa.
From visual inspection we observed that in some cases there was an image degradation comparing the 5 minutes and 1 minute datasets. This image degradation was clear in most cases where Bw similarity coefficient was employed, and also where RD potential function was used. The image degradation was mitigated in some cases when using HPa instead of HPf.
The mean ∆b between the 5 minutes and 1 minute scans, for all the analysed brain regions, is shown in figure 3. The numerical value for all bars with height <0.1 are shown numerically in the figure.
One consistent observation is that ∆b is lower for all methods with adaptive β value (HPa) than with fixed β value (HPf), indicating that HPa introduces less bias and makes the reconstruction algorithm more robust to different levels of statistics than HPf, which agrees with figure 2.
PLSS shows ∆b > 0.1 for all methods, indicating that its performance is highly dependent on the β value against the level of statistics.
We further analysed those methods that showed the smallest difference between the 5 minutes and 1 minute scans. These methods were OSL-PFQ-JEw, PGA-PFQ-Bw, PGA-PFQ-JEw and PGA-RD-JEw, all with HPa. Figure 4 shows the mean bias vs mean noise for 1-30 iterations, for the cortical regions.
All methods behave similarly at low number of iterations. However, at large number of iterations OSEM shows consistently higher bias and noise, while OSL-PFQ-  With regards to the lesion evaluation, figure 5 shows the resulting FoM for the largest lesion (11 mm diameter) for the 5 minutes and 1 minute scans. The axial slices of the reconstructed 5 minutes scans, centered where the lesions are located, together with the ground truth are also included in figure 5. The cases where the lesions were included in the T1w-MRI (T1w-les), absent in the T1w-MRI (T1w-noles), and present in the T1w-MRI and absent in the T2w-MRI (T1w-les T2w-no-les) were considered.
The lesion size error shows a consistent higher underestimation of the lesion size for the 1 minute scan, compared to the 5 minutes scan. As expected, the best results were obtained when the lesions were present in the T1w-MRI image, followed by the case where the lesions were present in the T1w-MRI and absent in the T2w-MRI, which resulted in similar lesions sizes as OSEM. In the scenario where the lesions are absent from the MRI images, even though the lesions are still clearly visible as observed in the axial slices in figure 5 (other similarity weights have shown to largely blur the lesions [7]), the lesion size is underestimated compared to OSEM. Lesion detectability, measured as the lesion-to-outer-contour activity concentration ratio, was degraded from the reference 2.25, down to 1.7 for the 5 minutes scan with OSEM (1.6 for 1 minute scan). The detectability with all the methods with the lesion present in the MRI was in general slightly higher as the one obtained with OSEM for the 5 minutes and 1 minute scans, except for the 1 minute scan using PGA-RD-JEw. The detectability was degraded when the lesion was absent in the MRI images compared to OSEM. These small differences are attributed to the PSF modeling included in all the reconstructions. The PSF introduces an overshoot in the lesions, but the MRIguided PET reconstruction partly corrects for the PSF overshoot when the lesion is present.
Regarding the simulated brain with the reduced activity concentration in the frontal cortex, figure 6 shows a reconstructed axial slice with each of the evaluated methods.
The simulated mean activity concentration reduction, compared to the rest of the gray matter, was 20%. The resulting mean activity concentration reduction was 22.4% for OSEM, 20.7 % for PGA-PFQ-JEw, 21.9% for OSL-PFQ-JEw, and 20.7% for PGA-RD-JEw.

Measured FDG Data
Through different simulated scenarios we determined that OSL-PFQ-JEw and PGA-PFQ-JEw performed similarly, while PGA-RD-JEw performed slightly differ- ently in the scenario with lesions. OSL-PFQ-JEw and PGA-RD-JEw were further compared with OSEM using measured FDG data.
The example shown in figure 7 is an FDG patient diagnosed with Epilepsy, scanned for 30 minutes in the Biograph mMR. We compared OSL-PFQ-JEw and PGA-RD-JEw with OSEM, in all cases with PSF. The dataset was reframed from 30 minutes, down to 5 minutes and 1 min. The original and decimated datasets were reconstructed using 1-10 iterations, to see the level of degradation at a high number of iterations. Subfigure 7a shows the central axial slices of the patient for the 30 minutes, 5 minutes and 1 minute, at 3 iterations (clinical standard) and 10 iterations.
The quantitative evaluation was performed similarly to the simulated data, described in section 1.4, to obtain the mean intensity value and noise of the VOIs. Subfigure 7b shows the mean intensity versus the noise, measured as the CoV, for OSEM, OSL-PFQ-JEw and PGA-RD-JEw. We show some representative VOIs such as the caudate, which is a small sub-cortical region part of the striatum, the frontal cortex which is a large region, the cerebellum which is often used as reference in brain studies, and the ventricles, which are a cold region filled with cerebrospinal fluid.
From visual observation we can determine that the PGA-RD-JEw method shows image degradation for the 5 minutes and 1 minute scans. In addition, the MRI information is more enforced on the PET image (for the 5 minutes and 1 minute scans) for the PGA-RD-JEw than for the OSL-PFQ-JEw.
On the other hand, OSL-PFQ-JEw shows almost no visual difference between the 30 minutes and 5 minutes images at 3 and 10 iterations. For the 1 minute scan the image is degraded, but most of the structure is visible, as compared to OSEM, where the image is too noisy to perform any kind of analysis. From figure 7b we observe that PGA-RD-JEw shows non-monotonic trends for the 1 minute scan for every region, which is in agreement with the degraded images observed in figure 7a, while for the 30 minutes scan, the behavior is similar to OSL-PFQ-JEw but with slightly higher noise.
The noise with OSL-PFQ-JEw is consistently lower than OSEM for every scan time and region, with the exception of the cerebellum where they are comparable for the 30 minutes scan, and the caudate and the frontal cortex, which is slightly lower for OSEM for the 30 minutes scan. These small differences can be attributed to the 90% confidence of the voxels belonging to each mask, which in turn reduces the impact of partial volumes.

Measured Florbetapir Data
The reconstruction methods used for the [ 18 F]-Florbetapir datasets were OSEM and OSL-PFQ-JEw, since PGA-RD-JEw showed image degradation with the FDG dataset, and PGA-PFQ-JEw showed similar performance compared to OSL-PFQ-JEw. Figure 8 shows the axial and sagittal views of two [ 18 F]-Florbetapir cases.
From visual inspection we can observe that the reconstructed images with OSL-PFQ-JEw between the 20 minutes and 5 minutes show high resemblance, while OSEM shows a clear noise increment in the 5 minutes scan. However, the 1 minute scan shows an important image degradation with both reconstruction methods. There are some visible hot spots/clusters in the OSL-PFQ-JEw reconstructed images, which are also visible in the OSEM reconstructed images. However, because the noise is in general lower with OSL-PFQ-JEw, the hot spots/clusters stand out. The quantitative comparison per regions between the 20 minutes scans and the decimated 5 minutes and 1 minutes scans shown in figure 8b, shows that OSL-PFQ-JEw produces overall more similar quantitative values among scan times than OSEM. The differences are more prevalent for the comparison between 20 minutes and 1 minute scans than for the 20 minutes and 5 minutes scan. The overall small differences were partly attributed to the large size of every region, being the smallest 7.9 mL for the frontal lobe.

Discussion
In this work different implementations of MAP MRI-guided PET reconstruction were explored, combining different potential functions, similarity coefficients and optimization algorithms. More importantly, each of these variations was implemented with a fixed and a variable hyperparameter value adapted to the image noise.
The prior model was constructed using the quadratic (PFQ) and relative differences (RD) potential functions, in combination with the traditional Bowsher similarity coefficient (Bw), and a variation based on joint entropy (JEw). These prior models were implemented within three different optimization algorithms: the one-step-late (OSL), a modified gradient ascent method (PGA), and an additional method based on a separable surrogate function to optimize (PLSS).
Performance achieved by each of these implementations was first evaluated based on different datasets obtained through Monte Carlo simulations. Results between the different reconstruction methods were compared using a 5 minutes and 1 minute simulated brain scan. The comparison was performed visually, and quantitatively focused on different cortical and sub-cortical regions from the Hammers atlas. In all cases, the hyperparameter was individually optimized to obtain optimum image quality for the 5 minutes scan.
To evaluate which method is less sensitive to the level of statistics we compared the mean intensity value of several brain regions between the 5 minutes and 1 minute scans. In the case of using the fixed hyperparameter value, it was not modified for each level of statistics.
The mean difference for OSEM was 0.2%. For every method with MAP, the difference between the 5 minutes and 1 minute scans was consistently lower for the adaptive hyperparameter value than for the fixed hyperparameter value. The smallest differences were measured in OSL-PFQ-JEw (-3%), PGA-PFQ-JEw (0.2%), PGA-PFQ-Bw (2%), and PGA-RD-JEw (0.8%). Out of these methods, noise and bias were evaluated for 1-30 iterations, concluding that the best results were achieved with OSL/PGA-PFQ-JEw and PGA-RD-JEw. Out of these methods, we could observe slightly noisier images using PGA optimization than OSL optimization at high iterations.
In the case were lesions were inserted in the simulated brain phantom, small differences were observed among MAP methods. Overall, having the lesions in the MRI improved lesion size and detectability for all methods. Having the lesions in one MRI sequence but missing in the other one, showed how none of the methods degraded the results compared to OSEM. When the lesions were missing in the MRI, all methods under-performed compared to OSEM, but PGA-PFQ-JEw showed to be the method that underestimated the lesion size the least.
The dataset where the frontal cortex had 20% reduced activity concentration showed lower bias with MAP reconstruction, without important differences between them, compared to OSEM.
Measured data were evaluated using OSL-PFQ-JEw and PGA-RD-JEw. PGA-PFQ-JEw was excluded because of its similar performance compared to OSL-PFQ-JEw. An FDG dataset, obtained from a 30 minutes acquisition, was decimated to 5 minutes and 1 min. Images were compared at 3 iterations (21 subsets), and at 10 iterations.
Results showed that PGA-RD-JEw is sensitive to the level of statistics. The datasets for the decimated 5 minutes and 1 minute showed artifacts in the PET image produced by the excessive enforcement of the MRI information on the PET image.
OSL-PFQ-JEw showed improved noise and spatial resolution compared to OSEM for the 30 minutes, 5 minutes and 1 minute datasets. It also showed high image quality at high number of iterations, demonstrating that OSL-PFQ-JEw can mitigate the noise as the number of iterations increase, still improving the spatial resolution.
Quantitative evaluation showed consistent significant lower noise using OSL-PFQ-JEw compared to OSEM (except the caudate for the 30 minutes scan), at comparable levels of mean intensity values in all analysed regions. The cerebellum showed the largest difference, with a difference of 2%.
The Florbetapir data showed that the observations made with FDG with regard to the improvements made by the MAP reconstruction compared to OSEM, can be extended to other radioisotopes presenting different uptake patterns. Also importantly, the clear positive case shows that there is no distinct pattern in the boundary between the gray and white matter in the PET image, which indicates that the MRI aids the PET reconstruction by reducing noise and enhancing the spatial resolution, without introducing any anatomical patterns that are not present in the PET image, which was one of the initial premises.
The quantitative evaluation showed that the mean intensity measured in the analysed regions from the datasets of a 5 minutes and 1 minute scan are more similar to the 20 minutes scan when using OSL-PFQ-JEw than with OSEM.
Even though the initial evaluation performed with simulated data was done based on the Biograph Vision, and the measured data were taken on a Biograph mMR, we think that the comparison between methods is directly related to the amount of net trues, not PET scanner.
One important remark to make is the importance of the need for an accurate registration between the PET and MRI. Given the voxel-to-voxel propagation of the MRI to the PET information, any misregistration due to patient motion should be below the voxel size (2.09 mm).
Potential applications of the presented MAP reconstruction expand from dose reduction, useful in pediatrics and repeated scans to study disease progression, to extraction of an image-derived input function for dynamic analysis of different radiotracers, where the image noise is usually high, especially in the first few frames.
It is also in the scope of our future work to determine the impact of dose/scan time reduction in clinical outcomes. It is worth exploring the ability to combine more than one MRI sequence in the MAP reconstruction, especially in neuro-oncology, where different sequences are usually acquired in order to study intratumor heterogeneity and surrounding tissue. Time-of-flight angiography can also be interesting in dynamic PET/MRI studies, to optimize the extraction of the image-derived input function.

Conclusions
In this work we explored different implementation alternatives for PET/MRI MAP MRI-guided PET reconstruction. We employed an adaptive hyperparameter value based on image noise which avoids the need for hyperparameter tuning. We conclude that OSL and PGA optimizations perform similarly when using PFQ as potential function and JEw as similarity coefficient, which were found to achieve the best image quality at any level of statistics used in this work. This conclusion was supported by a variety of simulated and measured datasets with different radioisotopes, for different scan times. We also conclude that the suggested MAP reconstruction algorithm reduces noise and enhances spatial resolution, without introducing quantitative bias.