Simulations of a micro-PET system based on liquid xenon

Figure 1 shows the planned configuration of the LXePET scanner consisting of 12 trapezoidal sectors arranged in a ring geometry. The inner bore has 10 cm diameter and 10 cm axial length. The LXe is contained in a stainless-steel vessel thermally insulated by a vacuum space. Each sector is a LXeTPC viewed by two arrays of LAAPDs. The anode and cathode areas are 10 × 9.2 cm² and 10 × 3.2 cm², respectively, and the drift length is 11.2 cm. Each APD array consists of seven APDs with 16 mm diameter, nine APDs with 10 mm diameter and eight APDs with 5 mm diameter. Smaller APDs are used in the inner region to enhance the reconstruction where most of the events occur. Figure 2 shows the APD layout in one of the sectors.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The simulation of the LXe prototype was carried out with the Geant4 simulation package (Agostinelli et al 2003). A positron emitter (¹⁸F or ²²Na depending on the study) was simulated. Following the decay of the radioisotope, positrons with energy sampled from a continuous distribution of the beta decay process were generated and tracked until annihilation. To simulate the non-colinearity of the annihilation photons, a new process was created and integrated in Geant4. The new process simulates the positron annihilation in-flight according to the Geant4 annihilation process and replaces the Geant4 annihilation at rest with a model where the non-zero momentum of the electron–positron pair is taken into account. The interactions of the annihilation photons with the phantom and PET scanner were simulated with the low-energy package of Geant4. Energy and 3D position of every photon interaction in the LXe detector were recorded. The numbers of ionization charges Nⁱ_{e −} and scintillation photons S_i created in the interactions were calculated as Miceli et al (2011):

$$\begin{equation} N^i_{e-} = \frac{(1- {\rm Fr}^*) \times E_i^{G4}}{15.6eV} \end{equation} \tag{ 1 }$$

$$\begin{equation} S^i = \frac{(\upsilon + {\rm Fr}^*) \times E_i^{G4}}{15.6eV} \end{equation} \tag{ 2 }$$

where Fr* is the electron–ion recombination fraction, E^G4_i is the energy deposited in the interaction i and υ = 0.2 (Aprile et al 2007) is the ratio of the number of excitons and ion pairs produced. The electron–ion recombination fraction Fr* varies on an event-by-event basis. It was modeled as a Gaussian function centered at Fr = 0.24 with width ΔFr = 0.032 (Amaudruz et al 2009). Electronics and photo-detectors were not simulated directly. Instead, instrumental responses were parameterized in subsequent analyses as described in Miceli et al (2011). The parameters used in the simulation are listed in table 2.

Coincidence events were selected using a two-step procedure. The first stage of the event selection simulated the response of the detector trigger using only the information from the scintillation light. Events producing less than 5000 scintillation photons (corresponding to approximately 180 keV) were rejected. For each photon of each annihilation pair passing the first selection stage, we calculated the energy from the scintillation light corrected for the solid angle using the information of the position from the charge measurement (Miceli et al 2011) and we used the resulting value to calculate the light-charge combined energy as described in Aprile et al (2007). Events with combined energy 450–600 keV were kept. The first interaction points defining the lines of response (LOR) of the selected events were stored in a list-mode format. The Compton reconstruction algorithm described in section 2.3 was used to find the first interaction point for multi-sites events.

When a photon interacts in the detector, it can Compton scatter multiple times before being photo-absorbed. A 511 keV photon is roughly three times more likely to Compton scatter than be photo-absorbed when it first interacts in LXe. The simplest interaction configuration is the 1–1 case in which the detector registers only one discernible interaction point for each of the two photons, corresponding to photo-absorption without scattering. Practically, however, multi-hit scenarios such as 1–2, 1–3, 2–2, etc are more common, and must be taken into account, as they contribute to blurring of the image due to ambiguity in the location of the first interaction point. The goal of the Compton reconstruction algorithm is to sort through all the possible scattering sequences, determine the path that is the most probable, and define the most likely first interaction point and its associated LOR.

For each pair of photons interacting M − N times in the detector, with M representing photon 1 and N photon 2, and M ⩽ N, there are M!N! number of possible interaction sequences. For each sequence, a LOR check is first performed, determining whether the trajectory passes through the phantom. Then, if the sequence was found to be viable, Compton kinematics were used to compute a test statistic score associated with the sequence.

The Klein–Nishina formula determines the scattering angle based on the energy deposited:

$$\begin{equation} \cos (\theta _{E}) = 1 + mc^2 \times \big(E^{-1}_{\gamma i} - E^{-1}_{\gamma i+1}\big) \end{equation} \tag{ 3 }$$

where E_γi is the photon energy before the ith step given by

$$\begin{equation} E_{\gamma i} = E_{\gamma 1} - \sum _{j=1}^{i-1} {\rm d}E_j, \end{equation} \tag{ 4 }$$

where m is the electron mass, θ_E is the Compton scattering angle, dE_j is the energy deposited at the jth step and E_γ1 = 511 keV is the energy of the photon before it reaches the detector. Alternatively, the scattering angle θ_G based on the position of the interaction site is calculated as

$$\begin{equation} \cos (\theta _{G}) = \frac{ \vec{u_{i}} \cdot \vec{u_{i+1}} }{ |\vec{u_{i}}||\vec{u_{i+1}}| }, \end{equation} \tag{ 5 }$$

where $\vec{u_{i}} = (x_{i}-x_{i-1},y_{i}-y_{i-1},z_{i}-z_{i-1})$.

For each candidate interaction site, we could, in principle, determine if the sequence was the correct one by comparing the scattering angles computed using the energy deposited (θ_E) with the observed scattering angles given the geometric distribution of interaction sites (θ_G). In the ideal situation, the difference would be zero.

The ability to resolve the correct sequence, however, depends on the position and energy resolution of the system. A statistical weighting was used to account for instrumental resolution limits:

$$\begin{equation} \chi ^{2} = \sum ^{N-1}_{i=1}\frac{ \left( \cos (\theta _{E})_{i}-\cos (\theta _{G})_{i} \right)^{2}}{\Delta \cos (\theta _{E})^{2}_{i}+ \Delta \cos (\theta _{G})^{2}_{i}}, \end{equation} \tag{ 6 }$$

where the error terms are defined as Aprile et al (2008):

$$\begin{equation} \Delta \cos(\theta _E)^2_i = m^2c^4 \times \left( \frac{\sigma _{{\rm d}E_i}^2}{E_{\gamma i}^4} + \sigma _{E \gamma i+1}^2 \times \big(E^{-2}_{\gamma i} - E^{-2}_{\gamma i+1}\big)^2\right) \end{equation} \tag{ 7 }$$

and

$$\begin{equation} \Delta \cos (\theta _G)^2_i = \sigma _{\rm pos}^2 \times (\sigma g_{i,x} + \sigma g_{i,y} + \sigma g_{i,z}), \end{equation} \tag{ 8 }$$

where

$$\begin{equation} \fl \vec{\sigma g_i}= \left( \frac{\vec{u_{i+1}}}{|\vec{u_{i}}| \cdot |\vec{u_{i+1}}|} - \frac{\vec{u_i} \times \cos(\theta _G)}{|\vec{u_{i}}|^2} \right)^2 + \left( \frac{\vec{u_i}}{|\vec{u_{i}}| \cdot |\vec{u_{i+1}}|} - \frac{\vec{u_{i+1}} \times \cos(\theta _G)}{|\vec{u_{i+1}}|^2} \right)^2 . \end{equation} \tag{ 9 }$$

The error on the energy deposited at the ith step, $\sigma _{{\rm d}E_i}$, and the error on the photon energy after the interaction step i, σ_{Eγi + 1}, are given by

$$\begin{equation} \sigma ^2_{{\rm d}E_i} ={\rm ENC}_Q^2 + \Delta\, {\rm Fr}^2 \times {\rm d}E_i^2 , \end{equation} \tag{ 10 }$$

$$\begin{equation} \sigma ^2_{E \gamma i+1} = {\rm i} \times {\rm ENC}_Q^2 + \Delta\, {\rm Fr}^2 \times \sum ^i_{j=1}{\rm d}E_j^2 . \end{equation} \tag{ 11 }$$

Finally, the viable sequence with the lowest test statistic score was chosen by the reconstruction algorithm, and the associated LOR was defined and recorded. If no suitable interaction sequence was found, then the event was discarded. This reconstruction technique is similar to the one used in Oberlack et al (2000), Aprile et al (2008) modified for PET applications.

At high rates, fast scintillation light signals are used to roughly (1 cm³) localize the event in order to match correctly the light signal with the slowly drifting charge. Pile-up of events can occur in the small volume determined by the light localization region and may contribute to the count losses. In order to improve the count rate capability of the LXePET system, a pile-up event recovery method based on energy balance and proximity to the light signal was developed. The efficiency of the algorithm was found to be 99%, 95% and 89% for two, three and four events, respectively, of simulated pile-up. The fraction of pile-up events was evaluated by simulating a mouse and a rat phantom filled with water and ¹⁸F. The time of each decay was simulated using a Poisson distribution. The count rate correction factor for the pile-up with the recovery method, _p, is given by

$$\begin{equation} \epsilon _p = f_s + \sum _{k=2}^{4} f_e^k \times \mu ^k , \end{equation} \tag{ 12 }$$

where f_s is the fraction of pile-up free events, f^k_e is the fraction of k-events pile-up and μ^k is the efficiency of the pile-up recovery method for k-events pile-up.

The output of the simulation consisted of interaction steps for two types of events: singles where only one of the two photons interacted with the detector and double events where both photons reached the detector. These data are source activity independent and do not contain random coincidence events. In order to simulate the count rate performance of the LXe detector, the detection rates at different source activities and instrumental parameters, such as dead time and coincidence window, had to be calculated. Also, count losses due to the pile-up of events in the TPC were taken into account.

The calculation of detection rates was done by Poisson statistical modeling, taking into account the probability of each interaction type, and assuming that only events with exactly two emitted photons detected were selected. Given the trigger probabilities of detecting zero (P₀), one (P₁), and two (P₂) photons from annihilation, the trigger rates for true and scatter events C_{2, 0} and for random events C_{2r, 0} for a given source activity, A, and coincidence window, Δt, can be computed:

$$\begin{equation} C_{2,0} (A)= \frac{1}{\Delta t}\sum _{k=1}^{\infty } \frac{{\rm e}^{-\lambda }\lambda ^{k}}{k!}P_{2}P_{0}^{k-1},\qquad {\rm where }\quad \lambda = A * \Delta t \end{equation} \tag{ 13 }$$

and

$$\begin{equation} C_{2r,0} (A)= \frac{1}{\Delta t}\sum _{k=2}^{\infty } \frac{{\rm e}^{-\lambda }\lambda ^{k}}{k!}P_{1}^{2}P_{0}^{k-2} . \end{equation} \tag{ 14 }$$

Coincidence windows of 1, 3 and 6 ns were considered in these studies. The count rate for true and scatter events C₂(A) and for randoms C_2r(A) are calculated as

$$\begin{equation} C_{2} (A)= \frac{C_{2,0}}{1 + C_{{\rm total},0}\tau }\epsilon _{2} \times \epsilon _p^2, \end{equation} \tag{ 15 }$$

$$\begin{equation} C_{2r} (A)= \frac{C_{2r,0}}{1 + C_{{\rm total},0}\tau }\epsilon _{2r} \times \epsilon _p^2, \end{equation} \tag{ 16 }$$

where τ is the instrumental dead time, _p is the count rate correction factor for the pile-up (equation 12), and ₂ and _2r are the probabilities of a triggered event to pass the event selection criteria. C_{total, 0} is the total trigger rate including random coincidences. The ratios ₂ and _2r depend on the combined energy resolution and energy window threshold, as well as on the event reconstruction strategy used. They are calculated for each data set (simulated true plus scatter data set and random data set) as the number of events which have combined energy within the 450–600 keV energy window and define a LOR which passes through the phantom, divided by the number of triggered events. The random set was generated by combining single unrelated events in pairs. The first-stage trigger probabilities for zero, one, two photons detection, and second-stage event selection efficiencies are given in table 3. The trigger probability of detecting one or two photons is 60% for both the mouse and rat phantoms. The probability of detecting two photons depositing more than 180 keV is 22%, significantly higher for the mouse phantom than 13% found for the rat phantom due to the smaller amount of scattering produced by the mouse phantom. The amount of scattering is related to the size of the phantom.

Table 3. Trigger probabilities for zero, one, two photons (P₀, P₁, P₂), and probabilities (₂, _2r) of a triggered event to pass the event selection for non-random and random events.

		Mouse-like phantom	Rat-like phantom
Scenario	[P0] (%)	41.5	40.8
	[P1] (%)	36.3	46.6
	[P2] (%)	22.2	12.6
Efficiency	2 (%)	43.3	33.5
	2r (%)	3.39	6.26

Once the two final detection rates were calculated, a rate dependence could be applied to the output of the Geant4 simulation. This was done by scaling the simulated double (true and scatter) and random events (pair-wise combinations of single events) to obtain the total detection rate C₂(A) + C_2r(A). This scaling approach allowed us to use a single large set of simulation data to compute the behavior of the detector and its performance at various resolution limits and activities without the need to re-simulate under different detector parameters.

In order to preserve the high-resolution spatial information contained in the data produced by the LXePET scanner, we reconstructed the point source data and the micro-Derenzo phantom with a list-mode reconstruction algorithm. The main advantages of list-mode data reconstruction over rebinned data reconstruction are preservation of the maximum sampling frequency, and faster reconstruction for low-statistics frames. Data reconstructed with histogram-mode methods are compressed in the axial and radial directions to reduce the sinogram size and to accelerate the reconstruction with a consequent loss of axial and transaxial resolution (Arman et al 2005). This effect is particularly evident moving away from the axial axis in the transaxial plane. List-mode methods reconstruct the data event by event without the need of binning the data into space and time intervals thereby preventing information losses. The information preserving characteristic of list-mode reconstruction methods is particularly useful for high spatial and temporal resolution PET systems (Reader et al 2004). The computational time of histogram-mode reconstruction methods depends on the number of LOR in the sinogram, whereas reconstruction time of list-mode methods depends only on number of events recorded. List-mode methods are therefore preferred for high-resolution scanners where the number of LOR can be much higher than the number of recorded events (Arman et al 2005). List-mode image reconstruction methods are also favorable in time-of-flight PET (Pratx et al 2011), motion corrected PET (Lamare et al 2007) and dynamic and gated PET (Rahmim et al 2005). We used a 3D list-mode image reconstruction algorithm for PET based on the maximum likelihood expectation maximization (MLEM) approach (Shepp and Vardi 1982). As in Barrett et al (1997) and Parra and Barrett (1998) each detected LOR was considered as a unique projection bin with the number of counts in each projection bin g_i equal to 1. Using notations fⁿ_j and f^{n + 1}_j for the intensity vectors in voxel j for step n and the next n + 1 iteration estimates, the iteration step for the list-mode MLEM algorithm is equal to

$$\begin{equation} f_{j}^{n+1} = \frac{f_{j}^{n}}{{s_{j}}}\sum _{i}{p_{ji}}\frac{1}{\sum \limits _{k}{p_{ik}f_{k}^{n}}} , \end{equation} \tag{ 17 }$$

where p_ij is the value of the system matrix describing the probability that a given emission event i originates from a certain voxel j and s_j is the sensitivity value for voxel j. The list-mode MLEM used on-the-fly ray-driven forward and back projection with bilinear interpolation (Rahmim et al 2004). We used 20 MLEM iterations for the ²²Na point sources and 100 MLEM iterations for the Derenzo phantom. The voxel size was 0.15 × 0.15 × 0.15 mm³ and the image size was 360 × 360 × 360 voxels. The reconstruction time for point sources (5.5 million LORs on average) was less than 3 h on an Intel Xeon 2.00 GHz CPU (single core). The reconstruction speed of the list-mode MLEM algorithm can be further improved by using the ordered subsets approach (Hudson and Larkin 1994) and parallel processing.

The system performance was evaluated based on the NEMA standards (NEMA 2008). The only deviation from the NEMA protocol was the use of a list-mode MLEM reconstruction method instead of FBP reconstruction algorithm for the spatial resolution studies. As explained in section 2.6, we used the list-mode MLEM method in order to preserve the high-resolution spatial information of the scanner.

Sensitivity. The sensitivity of the system was determined with a Na²² point source embedded in a 1 cm³ acrylic cube. The source was stepped axially over the axial length of the scanner. The data were rebinned using the SSRB.

Scatter fraction and count rate performance. The scatter fraction and count rate performance were obtained with a mouse- and a rat-like phantoms. The mouse-like phantom was a 25 mm diameter and 70 mm length polyethylene cylinder with a 3.2 mm diameter hole drilled at a radial distance of 10 mm. A simulated 3.2 mm diameter and 60 mm long rod was filled with water and ¹⁸F. The rat-like phantom was a 50 mm diameter and 150 mm length polyethylene cylinder with a 3.2 mm diameter hole drilled at a radial distance of 17.5 mm. A simulated 3.2 mm diameter and 140 mm long rod was filled with water and the ¹⁸F. The data were rebinned using the SSRB.

Spatial resolution. The spatial resolution was obtained with the Na²² point source used for the sensitivity studies. The source was placed at two axial positions 0 and 12.5 mm and five radial positions 0, 5, 10, 15 and 25 mm. The simulated data were reconstructed with the list-mode MLEM iterative method.

Image quality: The image quality was studied with a micro-Derenzo phantom made from acrylic measuring 40 mm in diameter and 35 mm in length. Arrayed throughout the phantom were cylindrical rods of length 30 mm and diameters 1.6, 1.4, 1.2, 1.0, 0.8 and 0.6 mm. The rods were offset radially by 7 mm from the phantom center and filled with water and ¹⁸F. The rod-to-rod separation was set to twice the rod diameter. The simulated data were reconstructed with the list-mode MLEM iterative method.

The absolute sensitivity was calculated following the NEMA standard. A simulated ²²Na point source was used for this study. The point source was stepped axially through the scanner at 0.5 mm steps over an axial length of 100 mm. One million ²²Na decays were simulated at each step. The total absolute sensitivity for mouse applications was calculated by summing the sensitivity for the sinograms which encompass the central 7 cm. Since the axial extent of the scanner was less than the length of the rat phantom, we calculated the total absolute sensitivity for rat applications summing all the slices, as described in the NEMA standard. The absolute sensitivity at the center of the FOV (CFOV) for an energy window of [450 ,600] keV was 12.6 %. The sensitivity profile for all axial steps can be seen in figure 3. The total absolute sensitivity for mouse and rat applications were 9.4% and 7.2%. The total system sensitivity was 7.2%. For comparison, typical values of the absolute sensitivity at CFOV range from 3.4% for the micro-PET FOCUS-220 with 7.6 cm axial FOV, 250–750 keV energy window and 10 ns time window (Tai et al 2005) to 9.3% for the inveon system with 12.7 cm axial FOV, 250–625 keV energy window and 3.4 ns time window (Bao et al 2009).

Zoom In Zoom Out Reset image size

The scatter fraction and noise equivalent count rate (NECR) studies were carried out using the rat- and mouse-like phantoms following the NEMA protocol. For each phantom 50 million ¹⁸F decays were simulated. The list-mode simulated true plus scatter data set was arranged in sinograms (radial bin size 0.3 mm) and oblique slices were combined into 2D projections using the SSRB method with a 1 cm slice thickness. For each sinogram, all pixels located farther than 8 mm from the edges of the phantom were set to zero. The profile of each projection angle was shifted so that the maximum value was aligned with the central pixel of the sinogram. All the angular projections were then summed to generate a sum projection. All counts outside the central 14 mm band were assumed to be scatter counts. To evaluate the scatter inside the 14 mm central band, we used a linear interpolation. For each slice i, the number of scatter counts C_{scatt, i} was given by the total scatter counts in the sinogram (outside and inside the central 14 mm band) divided by the number of pairs in the data set. The total event count C_{TOT, i} is the sum of the pixels in the projections divided by the number of pairs in the data set. The scatter fraction is given by

$$\begin{equation} {\rm SF} =\sum _{i=1}^{{\rm NSlices}} C_{{\rm scatt},i}/C_{{\rm TOT},i}. \end{equation} \tag{ 18 }$$

The mouse (rat) scatter fraction was 12.1% (20.8%), of which 4.9% (10.5%) was due to scatter only and 7.2% (10.3%) was due to ambiguities in the Compton reconstruction algorithm. A future paper will deal with reducing the ambiguities. An example of Compton ambiguity involves multi-interaction events where one or both photons interact in only two locations and deposit the same amount of energy. To calculate the percentage of the scatter fraction due to Compton ambiguities, we selected only true events in the simulation data set.

The random set was arranged in sinograms (radial bin size 0.3 mm) and oblique slices were combined into 2D projections using the SSRB method with a 1 cm slice thickness. The number of random counts C_{random, i} for each slice was the total counts in the random coincidence sinogram within 8 mm from the edges of the phantom divided by the number of pairs in the random set.

The noise equivalent rate for each slice was calculated as follows, where C₂(A) and C_2r(A) are the rates previously calculated:

$$\begin{equation} {\rm NECR}_i(A) = \sum _{i=1}^{{\rm NSlices}} \frac{((C_{{\rm TOT},i}- C_{{\rm scatt},i})\times {C_2(A)})^2}{ C_{{\rm TOT},i} \times {C_2(A)} + C_{{\rm random},i} \times {C_{2r}(A)} } \end{equation} \tag{ 19 }$$

The NECR curves for mouse and rat phantoms are plotted in figure 4 for 1, 3 and 6 ns coincidence windows. The scatter fractions (SF), peak true counting rate (R_{t, peak}), peak noise equivalent count rate (R_{NEC, peak}), activity at which R_{t, peak} is reached, and activity at which R_{NEC, peak} is reached can be found in table 4 for mouse and rat phantoms and the three coincidence windows with an energy window of [450 600] keV. Figure 5 shows true, scatter, random, total counts and NECR as a function of activity for the mouse-like phantom with coincidence window 6 ns and dead time 0.2 μs. The simulated results show a similar usable range of activity compared with commercial micro-PET systems (1670 kcps at 130 MBq for a mouse phantom, a 350–625 keV energy window and 3.4 ns timing window—Inveon (Bao et al 2009)).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Spatial resolution was determined using the ²²Na point source with diameter 0.25 mm embedded in a 1 cm³ acrylic cube. A total of 50 million ²²Na decays were simulated and an energy window of [450,  600] keV was used. It was assumed that the source activity would be low enough that random coincidences could be ignored. The source was placed at two axial positions: 0 and 12.5 mm. Five radial positions were used for each axial position: 0, 5, 10, 15 and 25 mm. The data were reconstructed with the list-mode MLEM iterative method (voxel size 0.15 × 0.15 × 0.15 mm³, 20 iterations). The point spread functions were formed by summing one-dimensional profiles parallel to the direction of measurement and within two FWHM of the orthogonal directions. The FWHM and FWTM values were calculated through linear interpolation between adjacent pixels at one-half and one-tenth of the peak value in each direction. The point spread function for a point source at the CFOV is shown in figure 6.

Zoom In Zoom Out Reset image size

Radial, tangential and axial resolutions, reported as FWHM and FWTM, are given in figures 7–8. At the CFOV radial, tangential and axial FWHM resolutions of 0.6, 0.6 and 0.8 mm, respectively, were found. At 25 mm radial and 12.5 mm axial offset, radial, tangential and axial FWHM resolutions were 0.7, 0.7 and 0.8 mm, respectively. The results show a uniform resolution ⩽0.8 mm (FWHM) throughout the FOV in radial, tangential and axial directions. At the CFOV, the 2D FBP gave the same results of the MLEM algorithm. For comparison, typical values of spatial resolution for conventional micro-PET systems are 1.3, 1.3 and 1.5 mm (micro-PET FOCUS-220 (Tai et al 2005)). Also, the deterioration of the radial resolution towards the periphery of the FOV, which is common for crystal-based preclinical PET systems due to lack of DOI information, is absent in the LXePET.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Figure 9 shows a trans-axial slice (thickness 24 mm) of the micro-Derenzo phantom with cylindrical rods of length 30 mm and diameters 1.6, 1.4, 1.2, 1.0, 0.8 and 0.6 mm reconstructed with the list-mode MLEM method (100 iterations). The voxel size was 0.15 × 0.15 × 0.15 mm³. No attenuation or scatter corrections were applied. The source activity was low enough that random coincidences could be ignored. Rods of diameter 0.6–1.6 mm are visible.

Zoom In Zoom Out Reset image size