Hybrid PET/MR Kernelised Expectation Maximisation Reconstruction for Improved Image-Derived Estimation of the Input Function from the Aorta of Rabbits

Positron emission tomography (PET) provides simple noninvasive imaging biomarkers for multiple human diseases which can be used to produce quantitative information from single static images or to monitor dynamic processes. Such kinetic studies often require the tracer input function (IF) to be measured but, in contrast to direct blood sampling, the image-derived input function (IDIF) provides a noninvasive alternative technique to estimate the IF. Accurate estimation can, in general, be challenging due to the partial volume effect (PVE), which is particularly important in preclinical work on small animals. The recently proposed hybrid kernelised ordered subsets expectation maximisation (HKEM) method has been shown to improve accuracy and contrast across a range of different datasets and count levels and can be used on PET/MR or PET/CT data. In this work, we apply the method with the purpose of providing accurate estimates of the aorta IDIF for rabbit PET studies. In addition, we proposed a method for the extraction of the aorta region of interest (ROI) using the MR and the HKEM image, to minimise the PVE within the rabbit aortic region—a method which can be directly transferred to the clinical setting. A realistic simulation study was performed with ten independent noise realisations while two, real data, rabbit datasets, acquired with the Biograph Siemens mMR PET/MR scanner, were also considered. For reference and comparison, the data were reconstructed using OSEM, OSEM with Gaussian postfilter and KEM, as well as HKEM. The results across the simulated datasets and different time frames show reduced PVE and accurate IDIF values for the proposed method, with 5% average bias (0.8% minimum and 16% maximum bias). Consistent results were obtained with the real datasets. The results of this study demonstrate that HKEM can be used to accurately estimate the IDIF in preclinical PET/MR studies, such as rabbit mMR data, as well as in clinical human studies. The proposed algorithm is made available as part of an open software library, and it can be used equally successfully on human or animal data acquired from a variety of PET/MR or PET/CT scanners.

Accurate and precise quantitative biomarkers can be obtained by exploiting the pharmacokinetic information in the measured data [13]. is requires the estimation of the radiotracer concentration in the arterial blood plasma (input function). e gold standard for such measurement is blood sampling during the PET acquisition, via arterial cannulation [14]. Unfortunately, this technique is invasive and can be complicated, as it requires arterial blood samples in specific quantities and at precise times with corrections for delay and dispersion to account for the distance between the sampling site and the regions of interest (ROIs) [15].
A noninvasive technique is the image-derived input function (IDIF) [16] which uses a region of interest (ROI) to measure the uptake in the vessel over time. e IDIF is a simple way to calculate activity over time; however, it is challenging due to image-related issues. Firstly, the choice of the ROI has a very important impact, and nonaccurate ROIs will affect the measurement [17,18]. Other challenges are related to the use of MR images to extract the ROI because a potentially inaccurate registration between PET and MR images can lead to erroneous estimates of the activity in the chosen arterial ROI. With a hybrid PET/MR scanner, the problem of coregistration is expected to be minimised. e aforementioned problems are mostly related to the ordered subsets expectation maximisation (OSEM) method [19] which is usually followed by postreconstruction Gaussian filtering due to the high noise levels expected for the very short-time frames used for the IDIF estimation. OSEM with or without postfiltering has been shown to produce inaccurate values of IDIF with bias up to 30% propagating through the kinetic constant calculations [20]. In preclinical experiments, these issues can be even more challenging [13,21] because of the smaller size of animal vessel tissue, such as rabbit aortas, especially when they are performed with clinical scanners designed for larger human subjects. In this case, the PVE can be significant, as the diameter of the rabbit aorta is about 5 mm which is the same order of magnitude as the PET resolution.
Different studies have proposed methods for the use of IDIF by correcting or avoiding PVE [22][23][24][25][26]. Zanotti-Fregonella et al. [16] have shown in their comparison between cannulation-based and image-derived input functions that the use of high-resolution PET images is often not sufficient to avoid the need of cannulation to obtain a reliable IDIF. Moreover, the accuracy of the IDIF may vary between radiotracers and scanners. MR-guided techniques have been proposed and discussed [15], showing that erroneous registration between the PET and the MR images, as well as erroneous MR segmentation, can introduce an error in the IDIF estimation. e problem of PET/MR misalignment has been discussed for the kernel by Deidda et al. [27]. In this study, we apply a PET/MR-guided image reconstruction algorithm, hybrid kernelised expectation maximisation (HKEM) [28], to minimise PVE during the reconstruction step so that we can obtain more accurate IDIF estimates. In addition, to minimise the PET/MR misalignment, the HKEM-reconstructed image at the peak activity frame was used together with the MR image to extract the ROI to be used for the estimation of the input function. In this way, only a percentage of the maximum value is included in the ROI avoiding low-value voxels outside the carotid in case of PET/MR misalignment. e kernel method [29], which was first introduced in PET image reconstruction by Wang and Qi [30] and Hutchcroft et al. [31,32], makes use of only one prior information image, MR or PET, respectively. Furthermore, many other studies showing promising performances have appeared in the literature [33][34][35][36][37][38][39]. In contrast, the HKEM method, which we recently developed in the open-source STIR library [40], exploits both the PET and the MR coregistered images to derive PET information iteration after iteration.
e HKEM algorithm was introduced by Deidda et al. [28] as a method for improving PET image resolution and uptake recovery in PET/MR phantom experiments, as well as contrast and quantification of atherosclerotic plaque lesions in carotid arteries in clinical PET/MR studies-which could also be applied in PET/CT studies. In addition, it is a robust and stable method which gives consistent results across different datasets using the same parameter settings. In this paper, we focus on the quantification of the aorta IDIF of rabbits using 18 F-based radiotracers such as [ 18 F]-FDG and [ 18 F]-NaF, to extend the applicability and usefulness of our novel reconstruction algorithm. Here, we assume that if HKEM can recover the uptake while retaining satisfactory noise suppression for low-count PET acquisitions, it will also be capable of providing accurate IDIF estimates using a wide range of dynamic PET frame durations. e paper is structured as follows: Section 2 describes the datasets used to study image reconstruction, list mode (LM) subsampling, and the experimental methodology. Section 3 presents the results of the proposed method and comparison with different standard algorithms. e results are discussed in Section 4, and conclusions are drawn in Section 5.

Methods and Materials
2.1. Simulation. A realistic simulation was created using a model derived from real [ 18 F]-NaF rabbit data [41] and utilities implemented in the STIR library. e real data were acquired with the Siemens Biograph mMR scanner at Mount Sinai Hospital, NY, USA. e voxel size for the simulated image was 1.56 × 1.56 × 2.031 mm. e rabbit was a healthy subject and was anaesthetised before the scan. It was injected with [ 18 F]-NaF 170 MBq and scanned for 90 minutes. Different organs and tissues were segmented from the acquired MR UTE sequence, using 0.07 ms echo time. e original voxel size is 1.56 × 1.56 × 1.56 mm. It is then aligned to the PET field of view (FOV) and resliced to match the PET native voxel size, 1.56 × 1.56 × 2.031 mm 3 , and FOV size, 344 × 344 × 127 voxels. e same image is also used for the calculation of the kernel matrix. In particular, the abdominal aorta, kidneys, bladder, myocardium, lungs, stomach, and background were extracted as independent images. Each tissue type was segmented using a semiautomatic segmentation method in ITK-SNAP based on thresholding [39], and it was then used as a ROI in the real PET data to estimate the activity concentration over 45 frames organised as follows: 17 × 6 s, 4 × 15 s, 4 × 30 s, 4 × 60 s, 4 × 180 s, and 12 × 300 s. e measured values were then assigned to every tissue in the simulation.
In order to create the projection data, each simulated image is forward projected into the sinogram space. e attenuation sinogram is estimated using the attenuation coefficient, μ, map obtained from a Dixon MR sequence [42][43][44][45], and the precalculated hardware attenuation coefficients for the bed and coils. e projection data containing random events were estimated as a uniform sinogram containing 20% of the total number of events in the simulated acquisition sinogram. In order to estimate the scattered events, the Watson single scatter simulation was applied [46], and a mask obtained from the μ map was used for the tail fitting. At this point, the random and scatter sinograms were combined as an additive term in the emission sinogram to create the modelled prompts projection data. e final step was the simulation of Poisson noise from the prompts events. e above steps were repeated for each simulated frame image, and 10 noise realisations were created.

Real Rabbit Data.
e acquisition was carried out using the Siemens Biograph mMR at Mount Sinai Hospital, NY, USA. e rabbit was a healthy subject and was anaesthetised during the scan. It was injected with [ 18 F]-NaF 170 MBq for the first study and [ 18 F]-FDG 133 MBq for the second, both scanned for 90 minutes. e attenuation images were obtained from the scanner, which included the attenuation coefficient for bed and coils. e LM data were divided into smaller frames, to permit calculation of the input function. e tracer was injected during the first seconds of the scan. e MR part of the kernel matrix was obtained from a MR UTE sequence with 0.07 ms echo time, and the original voxel size was 1.56 × 1.56 × 1.56 mm. It was then aligned to the PET field of view (FOV) and resliced to match the PETnative z voxel size, 1.56 × 1.56 × 2.031 mm 3 , and FOV size, 344 × 344 × 127 voxels.

Reconstruction Setup.
All the datasets were reconstructed using HKEM with 21 subsets and 10 iterations. e PET image voxel, λ j , using the HKEM can be written as where k fj is the fj th element of the kernel, N j is the number of feature vectors related to voxel j, and α f is the kernel coefficient to be estimated iteratively as follows: with p ij being the system matrix and s i the additive term. e fj th element of the kernel consists in two components, and it can be written as follows: where is the kernel derived from the MR image and is the kernel component derived from the updated PET image. e Gaussian kernel functions have been modulated by the distance between voxels in the image space. e quantity x j is the coordinate of the j th voxel, n is the subiteration number, z (n) j and v j are the feature vectors that are calculated from the n th updated PET image and the MR image, respectively, and σ m , σ p , σ dm , and σ dp are the scaling parameters for the distances in (4) and (5). Note that the HKEM uses a voxel-wise kernel. is means that the feature vector assigned for each voxel contains only one nonzero element with the same voxel value. e kernel parameters were chosen in order to obtain the minimum RMSE in the aorta. e values of the kernel parameters were set as follows: N � 27, σ m � 1, σ dm � 3, σ p �1, and σ dp � 3 (the last two are only used by HKEM).
For comparison, the same datasets have been reconstructed also with 21 subsets and 10 iterations of OSEM with and without 3 mm FWHM Gaussian postfilter. ese methods are denoted as OSEM+G and OSEM, respectively, in this study. e selected number of subsets and the application of the Gaussian post-filter are considered as standard settings in clinical routine. All datasets were reconstructed using span 1.
Scatter correction was performed with the method described by Tsoumpas et al. [47] and Polycarpou et al. [48]. Randoms were estimated from singles, which were calculated from delayed events [49]. e procedures for these evaluations, including attenuation and normalisation corrections [50], make use of STIR.

Image Analysis.
e comparison was carried out in terms of the mean value for all of the short frames and datasets, and the bias was estimated for the simulation to assess the accuracy of the proposed method. e ROI was obtained using the HKEM-reconstructed image and the MR image as follows ( (iii) e ROI, A, is obtained by taking into account only the voxels with a value bigger than 75% of the maximum in order to optimize those affected by PVE where I is the index of the voxel. Quantitative comparison between algorithms was performed using the following figures of merit: bias where t k is the mean value of the target ROI at frame k, t jk is the value of voxel j within the ROI at frame k, and V is the number of voxels within the ROI. e ROIs obtained with the proposed method are shown for each dataset in Figure 2.

Simulation.
e IDIF estimates for the simulated rabbit data and the early and late frames for the IDIF are illustrated in Figure 3. In the same figure, the reconstructed images with OSEM, OSEM+G, KEM, and HKEM, at the peak frame (24-30 s), are shown. Figure 4 presents the line profile of the aorta estimated for the images, as reconstructed with all investigated methods, at two different positions (LP1 and LP2), while Figure 5 reports the median IDIF estimated over the ten noise realisations using the HKEM. e shaded region is the range of possible values over the 10 simulated datasets, and the dashed line is the true IDIF. Finally, Table 1 reports the percentage value of the mean, maximum, and minimum absolute bias over the frames and the noise realisations.
A voxel-wise analysis example is reported in Figure 6, where the 10 peak frame images were combined to extract the bias and the SD images for each algorithm. Figure 7 shows the comparison, on the bottom row, between the initial 200 s of the input function on the left, and the later section of the IDIF on the right. Moreover, to give an idea of the image quality, the reconstructed [ 18 F]-NaF images for the peak time are shown on the top. Figure 8 reports the line profile of the aorta in two different positions (LP1 and LP2) for the [ 18 F]-NaF peak images reconstructed with the investigated methods to illustrate in detail the differences between the images reconstructed with different techniques. Figure 9 gives an example of fused PET/MR image quality for all the reconstruction techniques.

FDG Study.
e IDIF was estimated for a [ 18 F]-FDG study in order to assess the method on a different tracer.

Discussion
In this study, we have proposed the use of our recently developed hybrid kernelised reconstruction algorithm HKEM, for the estimation of the IDIF in the aorta artery of rabbits having undergone [ 18 F]-FDG and [ 18 F]-NaF PET/ MR studies using a clinical PET/MR scanner. e study was driven by the fact that many applications, where dynamic PET is used to extract more accurate and precise kinetic imaging biomarkers, rely on the estimation of the IDIF which is problematic in preclinical studies due to extensive PVE. As a consequence, it is relevant to propose a method which provides accurate estimates of IDIF. e results in Figure 3 show that the proposed reconstruction method and ROI extraction provide accurate results for all time points. e mean, maximum, and minimum bias were also calculated over the frames and the ten noise realisations (Table 1).
We were able to obtain a mean bias of 5% using the HKEM with the maximum value being 16.1%. Note that due to the applied threshold in the definition of the ROI, the OSEM also provided accurate results although the dynamic PET image frames were very noisy, and thus it becomes challenging to accurately delineate the appropriate aortic input function ROI, which is crucial for the IDIF calculation. In addition, a 52% averaged CoV over noise realisations means that there is a probability of about 68% that the repeated measure will have a value within ± 52% around the mean. As a consequence, values with high bias are very likely with OSEM. e results suggest that MR information can provide substantial improvement in terms of PVE and noise suppression. Nevertheless, the inclusion of the PET functional information allows better accuracy at similarly low noise levels (Table 1), compared to KEM. Figure 4 shows the line profiles in two different points of the carotid for the image corresponding to the peak. Here, we can notice the better delineation of the aorta for both the KEM and HKEM MRguided techniques, thanks to the broader smoothness applied in the background tissue regions. It is also important to highlight that the extraction of the ROI from the OSEM image in Figure 1 would not be accurate, as the maximum value was very high due to noise. us, the 75% thresholding would only extract very few voxels, therefore causing up to 100% bias in the OSEM IDIF values despite being associated with high accuracy estimates. Figure 5 illustrates the median full IDIF estimate over the 10 realisations, and it is possible to notice the accuracy over time compared to the true values. A voxel-wise analysis example is reported in Figure 6, where it can be seen the better image quality of KEM and HKEM, with lower bias in the aorta and low SD overall. e ROI analysis was also performed on this image. e results reported in Table 2 Figure 6: Voxel-wise image analysis over the ten noise realisations using the peak frame. e top row shows the average absolute bias, and the bottom row shows the SD over the ten noise realisations. frame; however, the repeatability of the measure is around 3 times worse. When Gaussian filter is applied, the value is extremely biased with similar CoV to HKEM and KEM. e same analysis was applied to two real PET/MR rabbit datasets acquired with the Biograph mMR scanner, using [ 18 F]-NaF or [ 18 F]-FDG radiotracers. Figure 7 shows consistent results for the IDIF plots. e reconstructed images using the real data show regions of high uptake only in some places of the aorta, thereby demonstrating the benefit in contrast and resolution of exploiting a hybrid PET/MR kernel matrix. Figure 8 presents the line profiles obtained with the different methods, showing the good resolution of the aorta when using the HKEM method and the poor quality of the postfiltered OSEM which is highly affected by the PVE. In Figure 9, the fused PET/MR image is illustrated for each technique, confirming the better alignment of the aorta region between the PET and the MR images and the resulting higher PET image resolution and aortic contrast. Moreover, the comparison between the [ 18 F]-FDG and [ 18 F]-NaF PET/MR studies allowed to assess the feasibility and performance of HKEM in estimating the aorta IDIF for two of the most commonly employed radiotracers in oncology and cardiology. From the results in Figure 10, the benefit of the synergistic PET/MR information encoded in the kernel matrix is visible especially in the IDIF plot. ese results are also supported by the line profiles in Figure 11 showing a clear definition of the aorta for the proposed method and minimum spill-out of activity from the aorta. It is worth noticing that, for the real data, there are two peaks in the early frames IDIF; this is probably due to the fact that the injection was not continuous during the scan but there was a sudden stop making the uptake rate drop down in that specific time frame. We could show the IDIF with one peak by summing the frame associated with the first peak and the frame having low uptake; however, we think it is interesting to show the effect of a noncontinuous injection on the IDIF estimation. e input function represents a very crucial data component when estimating kinetic parameters, and its accurate estimation can become extremely challenging for small animal imaging due to the very small sizes of the associated aortic vessels. In this study, we proposed the use of PET/MR synergistic information for the more accurate  and precise extraction of the aortic ROI and IDIF estimation in the framework of the HKEM method. We demonstrated that, despite the small size of the rabbit aorta, it is feasible and promising to employ the HKEM method for the extraction of an aorta IDIF estimate of improved accuracy and reduced PVE even when using a clinical PET/MR scanner. In addition, the method described to extract the ROI is easy to use and implement as it only involves trivial mathematics    between matrices. It is worth mentioning that, although this study was performed with PET/MR data, it could also work with PET-CT data especially if the CT image to use as anatomical information is a CT angiography image.

Conclusion
In this investigation, we demonstrated that the HKEM method allows the more accurate extraction of the aortic ROI for improved IDIF estimation even when using a human hybrid scanner, compared to conventional OSEM or anatomically guided KEM reconstruction. Our findings were validated with both 10 simulated [ 18 F]-NaF PET/MR datasets as well as 2 real rabbit PET/MR studies. Further, the methodology can be applied to most of the available radiotracers and with PET-CT without any major modification. is technique can enhance the use of dynamic PET in the context of imaging biomarkers with direct pharmacokinetic information.

Data Availability
A demonstrative code for the creation of the simulated study, reconstruction, and ROI extraction is available in CODE OCEAN at https://doi.org/10.24433/CO.bde84e0c-4c73-47fa-8ba5-81fb8bd2af77. e real rabbit data used to support the findings of this study, however, have not been made available because the Translational and Molecular Imaging Institute Group, who provided the data, retains the right to publish the data before making them generally available.

Conflicts of Interest
e authors declare that there are no conflicts of interest regarding the publication of this paper.