Validation of a spatially variant resolution model for small animal brain PET studies

In small animal positron emission tomography (PET) studies, given the spatial resolution of preclinical PET scanners, quantification in small regions can be challenging. Moreover, in scans where animals are placed away from the center of the field of view (CFOV), e.g. in simultaneous scans of multiple animals, quantification accuracy can be compromised due to the loss of spatial resolution towards the edge of the FOV. Here, we implemented a spatially variant resolution model to improve quantification in small regions and to allow simultaneous scanning of multiple animals without compromising quantification accuracy. The scanner’s point spread function (PSF) was characterized across the FOV and modelled using a spatially variant and asymmetric Gaussian function. The spatially variant PSF (SVPSF) was then used for resolution modelling in the iterative reconstruction. To assess the image quality, a line source phantom in a cold and warm background, as well as mouse brain [18F]FDG scans, were performed. The SVPSF and the vendor’s maximum a posteriori (MAP3D) reconstructions produced uniform spatial resolution across the scanner FOV, but MAP3D resulted in lower spatial resolution. The line sources recovery coefficient using SVPSF was similar at the CFOV and at the edge of the FOV. In contrast, the other tested reconstructions produced lower recovery coefficient at the edge of the FOV. In mouse brain reconstructions, less spill-over from hot regions to cold regions, as well as more symmetric regional brain uptake was observed using SVPSF. The contrast in brain images was the highest using SVPSF, in mice scanned at the CFOV and off-center. Incorporation of a spatially variant resolution model for small animal brain PET improves quantification accuracy in small regions and produces consistent image spatial resolution across the FOV. Therefore, simultaneous scanning of multiple animals can benefit by using spatially variant resolution modelling.


Introduction
To improve the spatial resolution of positron emission tomography (PET) images, knowledge about the system point spread function (PSF) can be used in the reconstruction process. Using the maximum-likelihood expectation maximization (ML-EM) algorithm, the resolution model can be introduced in the detection forward model to account for the loss of spatial resolution (Qi et al 1998, Reader et al 2003.
A simple resolution model uses a spatially invariant and isotropic Gaussian, with a width approximating the PSF width at the center of the scanner field of view (FOV). However, since the PSF of PET scanners is usually spatially variant and its shape does not follow a symmetric Gaussian, the isotropic Gaussian PSF is an approximation. More realistic models, which capture the asymmetric and spatially variant nature of the PSF, have been proposed (Lee et al 2004, Cloquet et al 2010.
The determination of a realistic PSF can be performed using several methods. For example, radioactive point sources can be measured in several positions in the scanner FOV, and the PSF shape can be parametrized with respect to the spatial location Any further distribution of this work must maintain attribution to the author(s) and the title of the work, journal citation and DOI. (Panin et al 2006, Cloquet et al 2010. Monte Carlo simulations and analytical calculations of the detection process can also be performed to circumvent the experimental measurement of the PSF (Alessio et al 2006, Zhang et al 2010. The resolution model can be introduced in the system matrix either in the projection or image space. For the projection space approach, usually the sinogram is blurred with the projection space PSF. Since for listmode reconstruction no sinograms are calculated, an alternative to incorporate the resolution model in the projection space consists on repositioning the LORs according to the probable LOR position calculated from the resolution model. However, this approach has been shown to not perform well using ML-EM listmode reconstruction (Bickell et al 2016). On the other hand, resolution modelling in the image space can be easily adapted for ML-EM list-mode reconstruction, showing good image quality (Cloquet et al 2010), and in practice is easier to adapt than the projection space approach. It is considered that the image space approach is only an approximation of the true PET resolution model (Rahmim et al 2013). However, little difference has been reported in image quality between images reconstructed using the projection and the image space approach (Kotasidis et al 2011).
Although resolution modelling is relevant in preclinical PET due to the necessity to differentiate small structures present in small animals, most of the research in resolution modelling has been focused on clinical PET imaging (Alessio et al 2006, Panin et al 2006, Sureau et al 2008, Cloquet et al 2010, Rapisarda et al 2010, Kotasidis et al 2011, Ashrafinia et al 2017. Moreover, spatially variant resolution modelling would improve the image quality in scanner regions away from the center of the FOV, therefore allowing to perform scans of multiple animals in the FOV without compromising quantification accuracy. Spatially variant resolution modelling has been performed in small animal PET scanners (Qi et al 1998, Lee et al 2004, Bickell et al 2016 but not many studies demonstrate the advantages of resolution modeling for preclinical biological studies. A study by Aide et al (2010) explored the effect of spatially variant resolution modelling in oncological small animal studies, but using a clinical PET scanner. In another study by Spinelli et al (2009), simulations of a cardiac study were performed to study the effect of the resolution modeling in PET quantification.
In this work we adapted a state-of-the-art spatially variant resolution model to the Inveon microPET preclinical scanner in the context of mice brain scanning. Although mouse brain regions are small and some are on the limit of visualization due to the PET spatial resolution, to our knowledge, no previous studies have investigated the role of resolution modelling in mice brain scans. The performance of the method in terms of spatial resolution and quantification accuracy, was evaluated across the scanner FOV and compared to a spatially invariant Gaussian model, as well as the vendor's reconstructions. Scans of line source phantoms served to assess the image spatial resolution across the scanner FOV. Finally, [ 18 F]FDG mice brain scans of mice positioned at the center of the FOV were compared to scans of mice positioned in an off-center position, as encountered when scanning multiple animals.

Scanner
The scanner used for this work is the Inveon microPET scanner (Siemens Medical Solutions, Inc., Knoxville, USA) and consists of 80 rings of LSO crystals, with 320 crystals per ring (Bao et al 2009). The ring diameter is 16.1 cm and the axial extent is 12.7 cm. The detection area of each crystal is 1.51´1.51 mm. The transaxial FOV has a diameter of 9.9 cm. Images were reconstructed using a matrix size of 128×128´159 voxels with a size of 0.776´0.776´0.796 mm in the x, y and z directions respectively.

Experimental measurements
Twenty sodium polyacrylate point sources were soaked in [ 18 F]FDG (480 kBq) and pasted in a linear array on a foam holder. The grains had a size smaller than 1 mm. The point sources were placed at 5 mm intervals, covering the 9.9 cm of the transaxial FOV diameter along the x axis. The horizontal platform with the point sources was placed at the center of the FOV along the vertical and axial directions. The point sources were scanned for 10 min, acquiring in total 181 million coincidences.

Parametrization of the spatially variant PSF
We used the model proposed in Rapisarda et al (2010) where an asymmetric Gaussian, with different widths along the radial internal and radial external directions is used to model the asymmetric shape of the PSF observed towards the edge of the FOV (i.e. elongated tail towards the center of the FOV).The radial internal direction points towards the center of the FOV and the radial external in the opposite direction. To simplify the parametrization of the PSF, and due to the moderate change of the PSF width along the axial and tangential directions (Bao et al 2009), the PSF width parameters along the radial direction were the only parameters which varied with the radial distance. The PSF can then be parameterized using a radial (r), tangential (t ) and axial (a) coordinate system as: The PSF parameters ( ) s r i and ( ) s r e are the radial internal and external width parameters, s a and s t are the axial and tangential width parameters and m , r m a and m t are the mean radial, axial and tangential coordinates, respectively. (·) H is the Heaviside function. The normalization factor A is calculated such that the integral of the asymmetric Gaussian function is 1, given by: The asymmetric Gaussian model presented above is fitted to each of the 20 point sources in the image reconstructed without resolution modelling (see section 2.3). Initially, for every point source image, the voxel with maximum intensity is located and the 2D transaxial plane containing the maximum intensity voxel is extracted to fit a 2D Gaussian. From this Gaussian, the centroid of the point source in the xy plane is determined (figure 1a). Using the coordinates of the centroid (x, y), the radial coordinate is calculated as = + r x y 2 2 and the angle with respect to the x axis is calculated as . The image is then rotated about the z axis (using trilinear interpolation) to align the point source radial direction with the x axis and the tangential direction with the y axis (figure 1b): is the inverse rotation matrix about the z axis. The 3 orthogonal planes (xy, y-z and zx) are extracted from the point source image (figure 1b) and the model in (1) is fitted, considering only the 2 respective components of the model for each respective plane (radial-tangential, tangential-axial and axial-radial, respectively). The fitting was performed using the Levenberg-Marquardt algorithm of the nonlinear least-square curve fitting toolbox from Matlab (The MathWorks, Inc., Massachusetts, USA). Finally, a quadratic polynomial fit is performed on s i and s e in function of the radial coordinate r, and s a and s t is calculated as the average of all fitted values in all point sources.

Resolution model
It is assumed that the spatially variant PSF (SVPSF) acts as a blurring matrix on the unblurred image: jn is the blurring kernel for voxel j, N j is the neighborhood of the voxel j and l¢ is the unblurred image one wants to estimate. The values q jn of the blurring kernel can be calculated for any voxel using the model in (1) with the parameters calculated from the point source fit. Initially, a kernel size is defined and the coordinates of all voxels that belong to the kernel neighborhood N j are calculated according to the defined size (figure 1c). Then, from the central voxel coordinates the radial coordinate r j and the angle q j with respect to the x axis is calculated. Afterwards, all kernel voxel coordinates in N j are rotated using (4) to align the voxels' radial coordinate with the x axis and the tangential coordinate with the y axis (figure 1d).
Finally (1) is used to calculate the kernel values with m = r , e j s a and s t calculated form the fitted parameters.
Additionally, in order to compare the use of the SVPSF with a spatially invariant isotropic Gaussian kernel, a 3D isotropic Gaussian was fitted to a point source at the center of the FOV (CFOV) and the Gaussian width was calculated.

Image reconstruction
Resolution modelling was performed in the image space by using (5) in the forward model of the listmode reconstruction algorithm (Reader et al 2003): ,voxels) at iteration c, g lj is the intersection path of detected line of response (LOR) l with voxel j, s j is the sensitivity correction factor for voxel j, ℓ w is a weight factor for LOR ℓ (ℓ = ¼ 1 ,total number of LORs) to incorporate normalization correction, and the detected events are subdivided in subsets S m ( = ¼ m M 1, , , subsets). Attenuation correction is performed by precorrecting the LOR's with the attenuation correction factor A , l calculated from the mu-map of the object reconstructed from a CT scan.
Scans were reconstructed using algorithm (6) incorporating the spatially variant resolution model in (1) with the parameters calculated from the model fit to the point sources. In addition, list-mode reconstruction using (6) with a spatially invariant Gaussian resolution model with FWHM of 1.6 mm (width at the CFOV) and list-mode reconstruction using (6) without resolution modelling, were performed. These reconstructions are referred as SVPSF, Gauss16 and noRM, respectively and were performed with 16 subsets and 16 iterations. Finally, the vendor's software was used to perform 2 additional reconstructions, using the parameters suggested by the vendor. The first reconstruction consists on sinogram 3D histogramming followed by fourier rebinning and ordered subsets EM (OSEM) 2D reconstruction (16 subsets, 4 iterations). The second reconstruction consist on sinogram 3D histogramming followed by 3D MAP reconstruction (16 subsets, 18 iterations, target resolution of 1.5 mm FWHM). Standard parameters were used for both reconstructions. These reconstructions are referred as OSEM2D and MAP3D respectively.

Line sources
In order to validate the PSF model in (1), a line source phantom consisting of 9 crystal capillaries was used. Capillaries were placed parallel to each other and their long axis was aligned with the scanner axial axis. The spacing between capillaries was 1 cm along the x axis and the internal diameter of the capillaries was 1.5 mm. The capillaries were filled with [ 18 F]FDG (420 kBq/cc) and scanned during 10 min. The scan was reconstructed with the 5 methods detailed in the previous section.
The reconstructions were analyzed by first aligning the capillaries long axis with the axial direction.
Alignment was performed in the reconstruction with best spatial resolution (thinner line sources images) by manually positioning the capillary to a position in which all pixels with maximum intensity were located in the same transaxial pixel coordinate. This same transformation was applied to all reconstructions. The average of 35 transaxial planes was then calculated and the image was upsampled by a factor of 4 to a pixel size of 0.194 mm using trilinear interpolation. Then the FWHM for each of the 9 line sources was calculated in these images. For each line source the contour at 50% of the maximum intensity was obtained and the eigenvectors of the contour were calculated. Then, the chords of the contour along the 2 eigenvectors were obtained and the average of the 2 chords lengths was calculated as the line source FWHM. In addition, to calculate the deformation due to the parallax effect of the, ideally, circular shape of the line source in the transaxial plane, the ratio of the contour largest eigenvalue magnitude to the smallest was calculated. We refer to this metric as the eigenvalue ratio. Since the image deformation due to the parallax effect is seen as an elongation of circular objects along the radial direction (i.e. ellipsoidal deformation) for objects placed The image is rotated about the z axis to align the radial direction with the x axis and the tangential direction with the y axis. The z axis is aligned with the axial direction. From the rotated image, the 3 central orthogonal planes are extracted and a 2D fit of the PSF model is performed on each plane. This procedure is performed for every point source (different r ) and the PSF is parametrized with respect to the radial coordinate. (c) In a similar manner, the PSF kernels for every voxel C j are determined by first calculating the coordinate r j and q j of the central voxel. (d) All voxels in the neighborhood of C j are then rotated using ( ) q R . z j After rotation, the coordinates x, y and z of the voxels correspond to the radial, tangential and axial coordinates. Finally, the kernel is calculated with these coordinates evaluating the PSF model with parameters as shown in the figure. close to the edge of the FOV, eigenvalue ratio values larger than 1 reflect a more prominent deformation of the shape for round objects. Figure 2 details the analysis procedure.

Line source phantom recovery coefficient
A phantom with 5 capillaries embedded in a warm background was used to evaluate the activity quantification in small structures. Five capillaries, with internal diameters of 0.6, 0.8, 1.0, 1.5 and 2.0 mm, where placed in a plastic container with 18 mm diameter (Deleye et al 2013) with their long axis aligned with the scanner axial axis. The capillaries and background had an initial activity concentration of 291 kBq/cc and 76.6 kBq/cc respectively, i.e. an activity ratio of 3.8. The phantom was scanned at the center of the FOV (center), midway between the center and edge of the FOV (mid-off-center), and at the FOV edge (offcenter), during 10 min in each position.
The recovery coefficient (RC) was measured at each position following the NEMA protocol (Bao et al 2009). Transaxial planes along 1 cm were averaged to obtain a line sources image with less noise. In these images, a circular region of interest (ROI) with twice the diameter of the line source, was centered at each line source and the maximum activity concentration within the ROI was calculated. The ratio of the true (measured) line source activity concentration to the maximum within the ROI is reported as the RC for each line source. The RC is reported for the 5 reconstructions detailed in section 2.3.

Center and off-center mouse brain [ 18 F]FDG scans
In order to evaluate the discrepancy in brain regional quantification between mice positioned in a centered and off-center position, mice brain [ 18 F]FDG scans were performed. All animals were housed in a temperature-controlled room with a 12-hour light-dark cycle (food and water available ad libitum). The night before each scan, all animals were fasted for at least 12 h with free access to water. Isoflurane gas anesthesia was used for procedures (5% induction, 2% maintenance). The experiments followed the European Ethics Committee recommendations (Decree 86/609/CEE) and were approved by the Animal Experimental Ethical Committee of the University of Antwerp, Antwerp, Belgium (ECD 2018-37). Four mice were scanned twice, in a center and off-center position respectively. Mice were injected with [ 18 F]FDG through tail vein injection (14.0 ± 1.0 MBq). After an awake uptake period of 20 min, mice were anesthetized and scanned during 10 min in the center position. The 10 min PET scan was followed by an anatomical CT scan. Ten minutes after the end of the center scan, the mice were positioned in an off-center position (16 mm from the center of the FOV) and scanned during 15 min. The scan duration in the off-center position was prolonged to account for radiotracer decay. Mouse brain images were analyzed using the software PMOD 3.6 (PMOD technologies Ltd, Zurich, Switzerland). Initially, the image was upsampled to a pixel size of 0.2 mm and the brain was cropped. The brain image was manually registered to an [ 18 F]FDG  e 2 ) of the contour is calculated and the FWHM is calculated as the average of the cords length, while the eigenvalue ratio is calculated as the ratio of the eigenvectors magnitude (longest/shortest). brain template, followed by non-rigid registration. The average of the 4 mice brain images in the center position, as well as in the off-center position, is calculated for all 5 reconstructions. Additionally, using brain regions delineated in an MRI template, the uptake in the caudate putamen (CP, hot region) and in the septum (Sept, background region) was calculated, and the contrast recovery coefficient (CRC) was calculated as: where A CP and A Sept are the activity concentrations in the caudate putamen and septum respectively. For each mouse, the CRC was normalized with respect to the CRC of the SVPSF reconstruction in the center position. The volume of the CP and Sept regions is 29.9 and 4.84 μl respectively. Figure 3 shows the shape of 2 point sources in the transaxial plane for 2 different radial distances (20.5 and 45.7 mm) and the asymmetric Gaussian model (1) fitted on the point sources. The shape of the point source is more elongated along the radial direction at the radial distance of 45.7 mm compared to 20.5 mm, which the model can reproduce correctly. The difference between the experimental measurement and the model fit is between  5% at both radial distances. However, at r=45.7 mm the difference with the experimental measurement is higher than at r=20.5 mm. Figure 4 shows the radial internal and radial external width parameters from the model fit to the point sources as well as the quadratic fit performed on the values of the model fit. Data points follow the quadratic curves shape for both radial internal and external widths. The R-squared of the fit was 0.99 and 0.98 for s i and s e respectively.

Line sources
Reconstructions of the line source phantom using the 5 different reconstruction methods are shown in figure 5. Line sources show lower intensity towards the edge of the FOV using noRM, Gauss16 and OSEM2D reconstructions, and elongation of the rods shape in the xy plane is visible also at the edge of the FOV. In contrast, using SVPSF and MAP3D the intensity of the line sources looks more uniform along the x direction, and deformation at the edge of the FOV is not visible. Figure 6a shows the FWHM values of the line sources using the 5 different reconstructions. As it was qualitatively seen in figure 5, resolution (mean  STD FWHM) decreases towards the edge of the FOV using noRM (2.05  0.23 mm), Gauss16 (1.36  0.27 mm) and OSEM2D (1.86  0.27 mm), while more uniform resolution along the radial direction (x axis) is found using SVPSF and MAPD3D. However, the FWHM is larger using MAP3D (1.92  0.05 mm) than using SVPSF (1.19  0.06 mm). The deformation of the line sources shape towards the edge of the FOV using noRM, Gauss16 and OSEM2D is observed in the eigenvalue ratio plot (figure 6b). The minimum-maximum eigenvalue ratio values are 1.09-1.73, 1.13-2.12 and 1.13-1.82 for noRM, Gauss16 and OSEM2D. On the other hand, the eigenvalue ratio is smaller using MAP3D and SVPSF along the entire radial direction. The minimum-maximum eigenvalue ratio values are 1.01-1.25 and 1.04-1.29 for MAP3D and SVPSF respectively.  3.4. Center and off-center mouse brain [ 18 F]FDG scans Figure 8 shows the average of the 4 mice brain reconstructions using the different reconstruction methods, for mice scanned in the center and off-center positions. In all cases, reconstructions in the center position show brain structures with more detail than reconstructions in the off-center position. Reconstructions with OSEM2D and noRM show little structural detail and present large spill-over from high uptake regions to cold regions. Differentiation between hot and cold brain regions is more clearly observed in reconstructions using MAP3D, Gauss16 and SVPSF, with MAP3D presenting more spill-over. Symmetry between left and right brain uptake is improved using SVPSF, especially in the off-center position, in which, e.g., Gauss16 shows lower uptake   in the left hippocampus (outermost hippocampus) compared to the right hippocampus (indicated with white arrows).

Line source phantom recovery coefficient
With respect to the SVPSF reconstruction of the centered mouse, there was on average a loss in CRC of 39%, 81%, 13% and 78% for the MAP3D, OSEM2D, Gauss16 and noRM respectively (figure 9). For the SVPSF off-center reconstructions there was a reduction of 9% in CRC with respect to the center reconstructions. For the other off-center reconstructions, the reduction in CRC with respect to their center reconstruction was 10%, 25%, 38% and 33% for MAP3D, OSEM2D, Gauss16 and noRM respectively. For OSEM2D off-center and noRM off-center, even a negative CRC was obtained, indicating that the activity in the hot region (A CP ) was lower than in the cold region (A Sept ).

Discussion
A spatially variant resolution model was implemented on the Inveon microPET scanner. The asymmetric Gaussian model proposed in Rapisarda et al (2010) reproduced correctly the larger elongation of the PSF along the radial internal direction in comparison with the radial external direction. The measured PSF presented longer tails which were not reproduced by the Gaussian shape. This caused larger model errors away from the PSF center of about 5%. However, the model still performed well as demonstrated in phantom and animal experiments. Similar to what was found in Rapisarda et al (2010), a quadratic fit to the radial width parameters as a function of the radial distance describes the data accurately.
The line source phantom reconstructions clearly show the loss of spatial resolution towards the edges of the FOV when no resolution modelling is considered. For the list-mode reconstruction considering a spatially invariant Gaussian resolution model as well as with the vendor's OSEM2D reconstruction, the spatial resolution is lower at the edge of the FOV. For the listmode reconstruction using the spatially variant resolution model as well as for the vendor's MAP3D reconstruction, the spatial resolution is uniform along the radial direction. However, the MAP3D reconstruction produced a lower spatial resolution than the spatially variant resolution model.
The deformation of the circular cross-section of the line sources due to the parallax effect near the edge of the FOV is also corrected using the spatially variant resolution model and the MAP3D reconstruction, as can be seen in the small eigenvalue ratios at the edge of the FOV. The other reconstruction methods showed an increasing deformation (eigenvalue ratio) along the radial direction.
A challenging test, given the scanner spatial resolution (1.5 mm), was performed to quantify the recovery coefficient in small structures (lines sources) with high activity embedded in a warm background. Considering the center, mid-off-center and off-center positions, the RC is below 0.5 for the 0.6, 0.8 and 1 mm line sources using all reconstruction methods. The RC is more uniform in the 3 positions using MAP3D, SVPSF and noRM. However, the maximum RC in the 2 mm line sources using SVPSF is 0.89, while for MAP3D is 0.7, and noRM produces the smallest RC among all reconstruction methods. Gauss16 performs similar to SVPSF, but producing lower RC in the off-center position.
In the mouse brain reconstructions, improvement is observed using MAP3D, Gauss16 and SVPSF, compared to OSEM2D and noRM. Spill-over from hot regions to cold regions greatly deteriorates differentiation between brain structures in OSEM2D and noRM reconstructions. MAP3D shows less spill-over to cold regions compared to OSEM2D and noRM, but shows more spill-over than Gauss16 and SVPSF. Symmetry between left and right brain regions uptake is better Figure 8. Average of the 4 mice brain reconstructions using the 5 different methods, in the center (left) and off-center (right) positions. Relative standard uptake value (SUVr, normalized to whole brain) is used for visualization purposes. The [ 18 F]FDG PET and MRI templates used for image registration and brain region delineation are shown at the bottom. Caudate putamen and septum regions, used to calculate the CRC, are shown in magenta and blue respectively. using SVPSF than using Gauss16, observed for example as a more symmetric uptake in the left and right hippocampus.
The contrast recovery coefficient in the mouse brain is reduced in the off-center position compared to the center position in all 4 mice using OSEM2D, Gauss16 and noRM. Only for MAP3D (mouse 1) and SVPSF (mouse 1 and 2) the CRC increase in the offcenter position. The best CRC, and the smallest difference between center and off-center CRC, is observed using SVPSF for all mice.
The possibility to obtain reproducible image quality along the entire scanner FOV allows to scan multiple animals at the same time (e.g. side by side) without compromising accuracy in quantification. Particularly for mouse brain PET scans, quantification in small brain regions is challenging. Using the spatially variant resolution model implemented here, brain regional quantification is minimally affected by the loss of spatial resolution observed at the edge of the scanner FOV. Quantification in mouse brain studies is therefore expected to improve using the proposed reconstruction, in both single and multiple animal scans.

Conclusions
List-mode reconstruction with spatially-variant resolution modelling was implemented in the Inveon microPET scanner. Compared to reconstructions without resolution modelling, or using a spatially invariant Gaussian model, the spatially-variant resolution model produced consistent spatial resolution along the entire scanner field of view. Correction for the parallax effect was also achieved using the spatially variant resolution model. In mouse brain reconstructions, improved symmetry and contrast is observed in the mouse brain uptake at an off-center position using the spatially variant resolution model. The use of a spatially variant resolution model in small animal scanners allows to obtain an improved resolution in mouse brain scans, especially if the brain is positioned at an off-center location. This method therefore allows to scan 2 mice positioned side-by-side at the same time without major resolution loss.