Efficient simulation of voxelized phantom in GATE with
embedded SimSET multiple photon history generator
(Running heads: Efficient simulation of voxelized phantom in GATE)
Hsin-Hon Lin
a, Keh-Shih Chuang
a*, Yi-Hsing Lin
b, Yu-Ching Ni
c, Jay Wu
d, Meei-Ling Jan
ca
Department of Biomedical Engineering & Environmental Sciences
National Tsing-Hua University, Taiwan
b
Health Physics section
Kuosheng nuclear power plant, Taiwan Power Company
c
Health Physics Division, Institute of Nuclear Energy Research
Atomic Energy Council, Taiwan
d
Department of Biomedical Imaging and Radiological Science
China Medical University, Taiwan
*Corresponding author:
Keh-Shih Chuang, PhD
Professor
Department of Biomedical Engineering & Environmental Sciences
National Tsing-Hua University
101, Sec. 2, Kuang-Fu Rd., Hsinchu 30013, Taiwan
Tel.: +886-3-574-2681; fax: +886-3-571-8649.
Abstract
GEANT4 Application for Tomographic Emission (GATE), is a powerful Monte Carlo simulator that combines the advantages of the general-purpose GEANT4 simulation code and the specific software tool implementations dedicated to emission tomography. However, the detail physical modelling of GEANT4 is highly computationally demanding, especially when tracking particles through voxelized phantoms. To circumvent the relatively slow simulation of voxelized phantom in GATE, one can use another efficient Monte Carlo code to simulate photon interactions and transport inside a voxelized phantom. Simulation system for emission tomography (SimSET), a dedicated Monte Carlo code for PET/SPECT systems, is well-known for its efficiency in simulation of voxel-based objects. An efficient Monte Carlo workflow integrating GATE and SimSET for simulating pinhole SPECT has been proposed to improve the voxelized phantom simulation. Although the workflow achieved a desirable speedup, it sacrifices the ability to simulate decaying radioactive sources such as nonpure positron emitters or multiple emission isotopes with complex decay scheme and lacks the modeling of time-dependent process due to the inherent limitation of SimSET photon history generator (PHG). Moreover, a large volume of disk storage is needed to store the huge temporal photon history file produced by SimSET that needs to be transported to GATE. In this study, we developed a multiple photon emission history generator (MPHG) based on SimSET/PHG to support most of medically important positron emitters and incorporated the new generator codes inside GATE to improve the simulation efficiency of voxelized phantom in GATE while eliminating the need for the temporal photon history file. The validation of this new code based on a MicroPET R4 system was conducted for 124I and 18F with mouse-like and rat-like phantoms. Comparing GATE/MPHG with
GATE/GEANT4 indicates that there is a slight difference in energy spectra for energy below 50 keV due to the lack of simulation of x-ray from 124I decay in the new code. The spatial resolution, scatter
fraction, count rate performance are in good agreement between the two codes. For the case studies of
18F-NaF (124I-IAZG) using MOBY phantom with 111 mm3 voxel sizes, the results show that
GATE/MPHG can achieve acceleration factors of approximately 3.1 (4.5), 6.5 (10.7) and 9.5 (31.0) as compared with GATE with regular navigation method, compress voxel method and parameterized tracking technique, respectively. In conclusion, the implementation of MPHG in GATE allows for improvement of the efficiency of voxelized phantom simulations and is suitable for studying clinical and preclinical imaging.
1. Introduction
Monte Carlo (MC) simulation in emission tomography is an essential tool to hasten the development of PET and SPECT imaging systems and the corresponding image reconstruction and correction methods, as well as the optimization of data acquisition and processing protocols (Buvat and Lazaro, 2006; Zaidi, 1999). Among these MC codes in emission tomography (Buvat and Castiglioni, 2002), GEANT4 Application for Tomographic Emission (GATE) (Jan et al., 2004), is a powerful Monte Carlo simulator that combines the advantages of the general-purpose GEANT4 simulation code (Agostinelli et al., 2003; Allison et al., 2006) and the specific software tool implementations dedicated to emission tomography (Pietrzyk et al., 2012). The platform can model time-dependent process such as decay kinetics, dead time and movement, while benefitting from the same versatility and support as that of the general-purpose simulation codes.
However, the detail physical modelling of GEANT4 is very computation-demanding especially in tracking particles through voxelized phantoms. The underlying GEANT4 behind GATE offers three methods to volume arrangement in voxelized phantom for tracking particle: G4VPVParameterised, G4VNestedParameterisation and G4PhantomParameterisation classes (Schümann et al., 2012). In GATE, the general particle tracking algorithm (called parameterized tracking) is developed based on G4VPVParameterised class. As the method regard each voxel in a voxelized phantom as a large of cuboidal volume set, it needs to spend much time on the determining the next voxel from the total volume set and updating the paths of particle at each voxel boundary encountered, resulting in the low efficiency for simulation. To tackle the inefficient simulation in voxelized phantom, two approaches using compressed voxels techniques (Taschereau and Chatziioannou, 2008) and improved voxel navigation (called regular navigation algorithm) (Rehfeld et al., 2009) had been proposed in GATE. The compressed voxels technique optimize the original parameterized tracking method by combining adjacent voxels with identical physical properties into larger voxels. This method can reduce the memory and central processing unit requirements for high resolution objects. Since GEANT4.9.1, a
new navigation algorithm based on G4PhantomParameterisation classes (Arce et al., 2008), is developed for the tracking of particles in voxelized volumes. The regular navigation algorithm optimizes the voxel navigation by only searching the neighboring voxel without a large memory overhead and skipping the steps when the material of next voxel is the same. The improved navigation is introduced in GATE and can speed up to about 2~7 (dependent on voxel sizes) compared to the compressed voxels technique.
To circumvent the relatively slow simulation of voxelized phantoms in GATE, one can also use other efficient Monte Carlo codes to simulate photon interactions and transport inside a voxelized phantom. Simulation system for emission tomography (SimSET) (Harrison et al., 1993), a dedicated Monte Carlo code for PET/SPECT systems, is well-known for its efficiency in the simulation of voxel-based objects (Barret et al., 2005). It uses variance reduction techniques such as forced detection (Haynor et al., 1990) and important sampling techniques (Haynor et al., 1991) to improve efficiency. An efficient Monte Carlo workflow integrating GATE and SimSET for simulating pinhole SPECT was proposed (Chen et al., 2008; Mok et al., 2010) and significant computational speedup of about 10 folds was achieved by the workflow. Although the workflow achieved desirable speedup, it sacrificed the ability of simulating the decaying radioactive sources such as nonpure positron emitters or multiple emission isotopes with complex decay scheme (Tang et al., 2009; Zhu and El Fakhri, 2009) and lacked the modeling of time-dependent process due to the inherent limitation of SimSET photon history generator (PHG). Moreover, high volumes of disk storage capacity for SimSET-GATE workflow were often needed to store the huge temporal photon history file produced from SimSET, even though the history file can be useful to efficiently study or optimize different detector configuration for saving the simulation time within the phantom. To overcome these problems, we developed a multiple photon emission history generator (MPHG) based on SimSET/PHG to support most of medically important positron emitters and incorporated the new generator codes inside GATE to improve the simulation efficiency of voxelized phantom in GATE.
The Monte Carlo modeling of MPHG will be described in the first part of this article, followed by the technical aspects of embedding MPHG in GATE. In the second part, we report the validation and efficiency of this new code based on different voxelized phantoms using various positron emitters as compared to GATE. Finally, we discuss some of the advantages and limitations of the new code.
2. Materials and methods
2.1. Multiple photon emission history generator
The basic architecture of the MPHG code remains identical to SimSET/PHG. Original PHG can be divided into three parts: (1) radioactive isotope simulation, (2) physical modeling and (3) photon tracking in voxelized phantom. In MPHG, major changes are made in part (1) to extend the capability of simulation of decaying radioactive sources with complex decay scheme, while preserving the efficient way of dedicated physical modeling and photon tracking.
(i) Decay scheme and isotope definition
A new input parameter file for isotope definition is needed for the MPHG codes. The file contains the half-life, atomic number and a decay scheme of each isotope involved. For each particle in the decay scheme, its type (positron or photon), energy, and branching ratio from high-state to low-state transition are defined. The decay scheme is modeled as a finite state-machine (Laedermann and Décombaz, 2000), as shown in figure 1. The model consists of a set of states and a transition function between those states and actions. Computation begins in the start state with an input state symbol of “0”. It continues changing to a new state depending on the corresponding transition function until reaching the final state. The transition function determines the next state sampled from a uniform distribution, which is apportioned pro rata among the subsequent states. The actions associated with state transition will generate the photon or positron, and then, the following cascade or annihilated photons are tracked through the attenuation distribution using the native photon tracking algorithm of
SimSET. It’s worth noting that once the action is electron capture, it has no particle generated but descends to a new state.
(ii) Photon cross section table
The photon cross section table in SimSET was based on PETSIM’s parametric model (Picard et al., 1992) incorporated with the Evaluated Photon Data Library’94 version (EPDL94) (Kaplan et al., 1997). The data files were trimmed to fit the main needs of the nuclear medicine, spanning energies from 1 keV to 1 MeV. However, some of cascade gamma ray’s energies were higher than 1 MeV for positron emitters with complex decay scheme. In order to accommodate the transport of cascade gammas with high energy, we built the new photon cross section table for energies ranging from 1 eV to 4 MeV based on the most recent EPDL97 version (Cullen et al., 1997). The EPDL database developed by the Lawrence Livermore National Laboratory (LLNL) includes cross section data for photoelectric absorption, coherent and incoherent scatterings, and pair production for energies ranging from 1 eV to 100 GeV. It has been recommended to serve as a standard and supersede earlier cross section libraries used in Monte Carlo codes to simulate medical imaging system (Zaidi, 2000). In addition, for most nonpure positron emitters, the energy of emitted gamma is lower than 4MeV and the cross section of pair production under the energy level is considerably small so that probability of pair production is negligible. Comparisons of attenuation coefficients (along with individual contributions from photoelectric absorption, Compton scattering and coherent scattering) in water as simulated by SimSET and the EPDL97 are plotted in figure 2.
(iii) Positron range
Positron range in SimSET (Harrison et al., 1999) is distributed according to an empirical model given by Palmer and Brownell (Palmer and Brownell, 1992). This model assumes that the equilibrium particle density resulting from sampling of beta-decay energy spectra can be represented by a three-dimensional Gaussian distribution centered at the origin. The SimSET package has included several
positron range tables for common nuclear medicine isotopes (18F, 11C, 82Rb and 68Ge). To extend the
variety of isotopes, especially in the non-pure isotopes of high atomic number, we adopted an analytical model for the positron emission energy density N(E). For a positron from an isotope of the atomic number Z with endpoint energy Emax, the probability density N(E) is given by
N
(
E)
=(
Emax−E)
2WpF(
Z ,W)
(2.3)where E is the emission energy in keV, W=1+E/511, p (
=
√
W
2−1
) is the momentum and F(Z,W) is the Fermi function.SimSET neglects the relativistic correction factor and employs a non-relativistic approximation form for the Fermi function, which is valid for conventional positron emitters of low Z, is given by:
F( Z ,W )=
2πη
1−exp (−2 πη)
(2.4)where
η=−
α ZW
p
, α (=1
137
) is the fine structure constant. However, the effect of the Coulombfield may significantly distort the beta spectrum for isotopes of high atomic number. In this work, MPHG used a more accurate approximation form with Coulomb correction factor (Venkataramaiah et
al., 1985):
F( Z,W )=
2 πη
1−exp(−2πη)
[
W
2
(
1+4γ
2)
−1
]
s (2.5)where
γ=αZ
andS=
(
1−α
2Z
2)
1/2−1
. Figure 3 shows the normalized positron spectra corresponding to 18F and 124I, with and without relativistic correction. It can be seen the differencerelativistic correction made is bigger for 124I due to its higher atomic number.
(iv) Temporal information and decay generation
Since the native photon history generator in SimSET didn’t support the temporal simulation, it tracked the given simulated decays and scaled the particle with a weighting factor calculated by the given simulated decay numbers, activity distribution and length of scan time. Therefore, it cannot describe
some physical mechanisms in terms of time-dependent phenomena. In addition, the decay generation in SimSET/PHG is voxel-by-voxel in sequence using a pre-calculated decay maps, it may conflict the randomized decay generation in GATE and affect the simulation results of time-dependent process, especially in coincidence detections of PET. Hence, additional options of randomized decay generation is implemented inside MPHG for connecting the MPHG and GATE seamlessly. To enable the temporal information in MPHG while keeping the synchronization with the time clock of GATE, we imported the decay time generated from GATE into MPHG when decay occurs, and the subsequent time of photon flight is calculated based on the travel length of emitted photon inside the object (Harrison et al., 2004).
2.2. Implementation of GATE/MPHG
In order to embed the MPHG in GATE, we developed the GateSourceVoxelSimSETInserter program inherited form GateSourceVoxellized class in the GATE code. The GateSourceVoxelSimSETInserter reads the SimSET parameters to initialize the MPHG engine. Once initialized, the routine GateSourceVoxellized::GeneratePrimaries (G4Event* event) will evoke MPHG to generate a cascaded decay process. The generated particles of each decay transport within the voxelized phantom as governed by SimSET. Once the tracked photons leave outside the object, the photons are projected and recorded onto a bounding cylinder (referred to as target cylinder in SimSET) and the cascade process of the event is finished in MPHG. Synchronously, GATE will produce the propagating photons at the same positions on the bounding cylinder that inherits all the photon information of the event from MPHG. Then, GATE tracks the propagated photons through the collimator and detector level according to the user’s instructions. The output file supports the LMF, ASCII, ECAT and Root formats according to GATE data output streams. The implementation flowchart of GATE/MPHG is illustrated in figure 4. At the time of this development, the MPHG is based on the version of SimSET/PHG 2.9.1, and GATE is 6.1.
2.3. Positron emitters
In this work, two representative isotopes of pure and nonpure positron emitters were simulated. 18F is
the most commonly used in PET for its attractive properties, one of which is that it can be easily substituted into biomolecules. 18F is almost pure positron emitter with 97% positron emission and 3%
electron capture. 124I radionuclide had been reported to be useful for imaging applications in oncology
as a surrogate for 131I prior to radioimmunotherapy. 124I has a complex decay scheme including 24
electron capture decay lines and 5 beta plus decay lines, resulting 81 cascade gamma emissions (Bhat, 1992). As a good compromise between computing efficiency and accuracy, we modeled the emitted particles if their absolute emission probabilities were more than 0.01 %. A total of 66 cascade gamma emissions were modeled, and the resulting gamma abundance of 124I was 99.53% in our simulations.
2.4. Simulation models of MicroPET R4 scanners and physics
The scanner geometry of MicroPET R4 system (Knoess et al., 2003; Lartizien et al., 2007) consisted of 4 rings with 96 detector modules, each coupled to an 8×8 array of 2.1×2.1×10 mm3 LSO crystals
resulting in 32 rings with 192 elements. The 96 modules were arranged in four rings of 24 modules each, with a diameter of 7.4 cm. The crystal pitches were 2.423 mm and 2.426 mm in the axial and transverse directions, respectively. The external lead shields were modeled with the opening of 120 mm at the front and 134.5 mm at the rear. The simulated energy resolution for each crystal was randomly sampled from a uniform distribution over the interval [17% 35%] and a 0.91 quantum efficiency was applied. A back compartment was also simulated to account for the back scattering induced by the light collection system. To compare the object simulation between GATE/MPHG and GATE/GEANT4 directly, the dead time models at different levels of the electronic read-out and transfer bandwidth were neglected in this study.
The underlying GEANT4 engine behind GATE is a general purpose Monte Carlo simulation package (Agostinelli et al., 2003). Three physics models in GEANT4 (Geant4.9.3.p02) are available for electromagnetic processes. The standard process model for photon interactions is employed in this study for its efficiency, while ensuring sufficient accuracy in PET simulation. To avoid infrared divergence, some electromagnetic processes (ex: gammas and electrons) require production thresholds below which no secondary particles will be generated. This threshold is defined as a range cutoff, which is internally converted to energy for each individual material. The range cutoff for all GATE system is set to 1 mm by default. In addition, GATE codes offers three options of source type for PET: ion source, particle source and back-to-back source. The ion source is used since it could afford the description of complex decay schemes with the most realistic way of simulating a radionuclide among these options. To validate the accuracy of the proposed MPHG model, the simulation results of energy spectra, spatial resolution, scatter and cascade gamma fractions, and count rate performance were compared with GEANT4 model.
2.5.1 Energy Spectra Energy spectra were acquired for the MicroPET R4 system using a 1 mm
diameter 124I line source in the center of a cylindrical water phantom (5 cm in diameter and 10 cm in
height). The acquisition of energy window and coincidence time window setting were set 1-1000 keV and 6ns, respectively. The energy spectra of detected photons with energy between 1 keV and 1000 keV were extracted from both MC simulation data for comparisons.
2.5.2 Spatial resolution The spatial resolutions for 124I and 18F were measured using an ideal point
source (10 MBq) embedded in a 10 mm acrylic cube phantom. Simulations were performed with the point source located at the radial distance of 0, 10 and 25 mm from the center. Data were acquired with an energy window of 350–750 keV and a coincidence time window of 6 ns. The simulated list mode data at each location were histrogramed into sinogram data sets with span 3 and maximum ring difference 31. The 3D sinograms were reconstructed by FORE + 2D FBP using a ramp filter with
cutoff at the Nyquist frequency and a zoom factor of 5 to achieve a reconstructed voxel size of 0.162 × 0.162 × 1.21 mm3. Data were not corrected for normalization, scatter, attenuation, and source
dimension. The response function was formed by summing all 1-dimensional profiles that were parallel to the radial, tangential, and axial directions. According to the NEMA procedure (NEMA, 2008), a parabolic fit of the peak point and its 2 nearest neighboring points were used to determine the maximum value of the response function. Linear interpolation between adjacent pixels was used to determine the position of 50% and 10% of the maximum intensity for the measurement of full width at half maximum (FWHM) and full width at tenth maximum (FWTM).
2.5.3 Scatter and cascade gamma fractions A 50 mm diameter, 150-mm long polyethylene (HDPE)
phantom (density 0.96 g/cm3) with a hole at 17.5 mm distance from the axis and a 25-mm diameter,
70-mm long HDPE phantom with a hole at 10 mm distance from the axis were used representing objects with size of a rat (referred to as the “rat-like phantom”) and a mouse (referred to as the “mouse-like phantom”), respectively. The line source with diameter of 1.0 mm was inserted in the drilling of each phantom. Both 18F and 124I were simulated on both rat-like and mouse-like phantoms.
The total activity for each simulation was set as 10 μCi, which was low enough to ignore the random coincidence (below 0.5% of the true event). Data were acquired in three energy window settings of 350-650 keV, 350-600 keV, and 425-650 keV. The simulated list-mode data were binned into 2D sinograms by single slice rebinning SSRB method (Daube-Witherspoon and Muehllehner, 1987) and sinogram of true, scatter, gamma and random components were stored separately. The gamma component was specific to nonpure positron emitters and was defined as the double coincidence involving at least one cascaded gamma originating from the same decay. For each simulation, all 63 slices in the sinograms were summed first and the total counts of the summed sinograms for each component were scored, separately. The gamma fraction was calculated as the ratio of the gamma counts to the total counts. However, in order to cross validate the real impact of scatter coincidence at the same phantom between 18F and 124I, the scatter fraction was calculated by taking the ratio of
scattered coincidences to the sum of scatter and true coincidence events. It’s worth noting that to separate the coincidence events involving cascaded gamma and annihilated gammas, the analysis module of GATE code was modified to tag the single events originating from annihilated or cascaded gamma ray photons (De Beenhouwer et al., 2009).
2.5.4 Count rate performance A count rate experiment was performed using the two aforementioned
phantoms (rat-like and mouse-like phantom) to validate the GATE/MPHG against GATE/GEANT4. These two phantoms with line source filled with 18F and 124I were simulated for activities varying from
0.2 to 8 mCi. Data were acquired with an energy window of 350–650 keV and two settings of coincidence time window (6 and 10 ns). The simulated list-mode data were rebinned into sinogram of true, scatter, gamma and random components separately. The count rates of prompts (total coincidence) and random coincidences were calculated as a function of the activity concentration.
2.6 Efficacy of GATE/MPHG model
To demonstrate the efficacy of GATE/MPHG in the efficient simulation of voxelized phantoms with reasonable accuracy, realistic simulations of biological distributions in MocroPET was performed. Current GATE platform supports three options to simulate the voxelized phantom, including general parameterized tracking (‘parameterized’), compressed voxel method (‘compressed’), and regular navigation method (‘regular’) for imaging applications. Another nested parameterized method based on G4VNestedParameterisation class was used for radiotherapy applications in GATE at the time of this development. A new version of GATE (version 7.0) had been released on 12 May 2014 which has unified the navigation engine for radiation therapy and imaging applications (OpenGATE collaboration, 2014).
The Mouse whole-body MOBY voxelized phantom (Segars et al., 2004) was used in comparing the efficacy among GATE/MPHG and the above three methods offered in GATE/GEANT4. All the
parameter settings for each method were under the same system configuration and acquisition protocols.
It’s worth noting that the tracking of electron in GEANT4 enabled the calculation of radiation dose and produced more accurate results with the cost of its efficiency. In contrast, SimSET is fast but cannot be used for dose estimation (Chuang et al., 2014). To be comparable with GATE/MPHG, which ignores the simulation of electrons, the electron range cutoff in phantom region of GATE/GEANT4 was set to 500 mm. This range cutoff setting is sufficiently large to suppress the production of electron in phantom. The settings of range cutoff outside the phantom still remain 1 mm. All simulations were performed on four 64 bit Fedora Linux 15 computers (AMD Phenom™ II X6 CPU @2.8 GHz, 8GB). Simulations of each method was assigned to one core in the computers and the CPU time was recorded.
The input attenuation distribution of MOBY phantom was segmented into six materials including air, lung, spine bone, rib bone, skull and soft tissue. Two simulations of different radiopharmaceuticals distribution using 18F-NaF and 124I-IAZG were presented in this study. The first simulation protocol
consisted in a static 5 min acquisition of the MOBY phantom corresponding to a radiotracer distribution of 10 µCi of 18F-NaF in the bones. The activity distribution in bone region was assumed
uniform. In the second simulation protocols, the radiotracer distribution of the MOBY phantom emulated a 124I-labeled iodo-azomycin-galactoside (124I-IAZG) uptake at 3 hours post injection in a
tumor-bearing mouse (Zanzonico et al., 2004). The spherical-shaped tumor with 2.5 mm radius was inserted in the right hind leg. The activity concentration ratios of the 450 µCi extended source distribution were 51:9:10:8:18:11:24:46:255:8 for tumor, blood pool, heart, lung, liver, spleen, kidney, stomach, intestine and background, respectively. Various image sets with voxel sizes of 0.125×0.125×0.125, 0.25×0.25×0.25, 0.5×0.5×0.5 and 1×1×1 mm3 (corresponding to voxel
dimensions of 256×256×900, 128×128×450, 64×64×225, and 32×32×113) were generated from MOBY software to investigate the influence of the voxel number on the simulation time.
3. Results
3.1 Validation of GATE/MPHG with GATE/GEANT4 model
3.1.1 Energy Spectra The energy spectra of detected photons for 124I are plotted in figure 6. Due to
the limited energy resolution of the simulated PET detector, there is no apparent energy peak at 603 keV, where high abundance of cascaded gamma with energy of 603 keV can be detected. The results shows good agreement between GATE/GEANT4 and GATE/MPHG models except for energies below 50 keV. The discrepancy at energies range below 50 keV is mainly due to the simulation of x-ray emissions for 124I in GATE/GEANT4, but not in GATE/MPHG.
3.1.2 Spatial resolution The FWHM’s and FWTM’s of the point source profiles at different radial
offsets from GATE/MPHG and GATE/GEANT4 simulations are listed in table 1. Good agreements in both FWHM and the FWTM values in all three directions between GATE/MPHG and GATE/GEANT4 can be found. Although the discrepancy between the two models for 124I is larger
than that for 18F, most errors in FWHM between the two models remain below 5%. Since the transport
of photons in MicroPET scanner for both models are based on the original GATE/GEANT4 codes, the cause of this difference on spatial resolution is mainly due to the different positron model employed between MPHG and GEANT4. The modeling of positron range in MPHG is based on the analytical modeling proposed by Palmer and Brownell (1992), while GATE is based on the Monte Carlo simulations using cross-sections for the annihilation process initially described in Heitler (1954). Recent reports have indicated the difference in positron range simulation between the two models is greater for isotopes with higher positron range (Lehnert et al., 2011).
3.1.3 Scatter fraction and cascade gamma fraction Scatter and/or gamma fractions of the line source
inside HDPE phantoms for 18F and 124I with various energy window settings are listed in table 2 and
table 3. For both isotopes, there is excellent agreement between the scatter fractions for the two phantoms at different energy window settings with GATE/MPHG and GATE/GEANT4. Also, the
gamma fractions of 124I on both phantoms agree well between GATE/MPHG and GATE/GEANT4.
The agreement also indicates that the neglect of pair production and non-essential cascade gammas at the energy window of typical PET in the simulation is practical and justified.
3.1.4 Count rate performance Figure 7 shows simulated prompt and random coincidence count rates
of 18F on the MicroPET R4 system for the mouse-like and rat-like phantoms as a function of the total
activity at 6 and 10 ns coincidence time windows. The figure shows that prompt and random count rate curve with various activities between GATE/MPHG and GATE/GEANT4 agree well. Figure 8 shows similar results for 124I. There is a minor systematic underestimation of count rates from
GATE/MPHG compared with those from GATE/GEANT4. The slight discrepancy might be due to the fact that the positron range is only modeled in the object in MPHG i.e., positrons leaving from the object are not propagated into collimator and detector of GATE. The detection of the escaping positron or their associated annihilated gammas increases additional detector signal and coincidence. This effect will be amplified when the positron rage of the isotope is large and the object size is small such as in the case of mouse-like phantom with 124I as shown in figure 8 (a)
and 8 (c).
3.2 Efficacy of GATE/MPHG model The summed sinograms of 18F-NaF and 124I-IAZG for the four
voxelized tracking methods were shown in figure 9 and figure 10, respectively. Figure 11 shows the profiles summed over all projection angles of the sinograms for the four methods. The profiles of the sinogram show good agreement among MPHG, compressed voxel, regular navigation, and the parameterized method. There is a larger average relative error (~8%) in the profiles of figure 11(b) for
124I-IAZG between MPHG and other tracking methods. This slightly larger discrepancy of MPHG
found in the profile of regions outside the phantom is due to the negligence of coincidence caused by the escaping positrons as mentioned in the previous section. The execution CPU times
as a function of voxel numbers in phantoms of 18F-NaF and 124I-IAZG were compared in figure 12(a)
and 12(b), respectively. The results indicated that the simulation time using parameterized and compressed voxel methods was increased more quickly with voxel grid size than the MPHG and regular navigation methods. These trends conform to those presented in Rehfeld et al (2009), although the scanners studied are different. The ratios of improvements in efficiency of GATE/MPHG as compared to the regular navigation method for 18F-NaF (124I-IAZG) under different voxel grid sizes
almost stay constant and 3.1 (4.5) in average. By contrast, the speedup ratios using GATE/MPHG in comparison to the compressed voxels and parameterized tracking method for 18F-NaF (124I-IAZG)
cases improve from 6.5 to 21.7 (10.7 to 30.0) and 9.5to 35.9 (31.0to 170.3) when the matrix sizes increase from 32×32×113 to 256×256×900. The improvement factors of GATE/MPHG for 124I are greater than that for 18F. It may be attributed to the efficient simulation of decay scheme in
MPHG in comparison to the time-consuming simulation process of radioactive source in GATE.
4. Discussion and Conclusion
The work presented here demonstrates a newly developed MPHG program based on SimSET/PHG to model the simplified decay process of radionuclide and the subsequent photon transport in object. Although no neutrinos and de-excitation process (Auger effect or x-ray emission) are simulated, the program can meet the most nuclear medicine simulation needs for a wide variety of radionuclides while preserving the efficiency of the simulator.
We compared our simulations with the well-validated GATE/GEANT4 simulation package for imaging 18F and 124I isotopes with MicroPET R4. There are in good agreement in terms of the energy
models. The main discrepancy of 124I energy spectra was observed in the energy range below 50 keV
due to the lack of simulation in de-excitation process. Fortunately, it didn’t affect the coincidence acquisitions within the PET energy window. A comparison of spatial resolution reveals that the GATE/MPHG positron range is comparable withGATE/GEANT4 positron rage for 18F, but is slightly
lower than GATE/GEANT4 positron range for 124I. Despite the underestimation of positron range for 124I in MPHG, the relative error of averaged spatial resolution in terms of FWHM and FWTM between
the two models remains below 5%.
MPHG uses an analytical model to calculate the positron range instead of directly simulating the transport of positron. Only photons leaving the object will be transferred from MPHG to GATE; no positrons (or other particles) will be transferred. As a result, the photons from positrons that would escape from the object and interact with the tomograph will not be simulated by GATE/MPHG, resulting in the slight underestimation of count rate compared with GATE/GEANT4. However, it is reasonable to neglect the escaping positron because they will be blocked by the outer casing of PET in real-world scenarios.
GATE is a highly object-oriented software; therefore MPHG can be readily embedded inside GATE codes without changing the original architecture of GATE. The new MPHG engine offers an alternative option in the simulation of voxel phantoms. In comparison to the regular navigation method, the speedup ratios using MPHG can achieve approximately 2.2 (3.5) in spatial resolution studies of 18F (124I), where only point source in air was simulated, and approximately 3.1 (4.4) in
scatter and gamma fraction studies of 18F (124I), respectively. Using realistic MOBY phantoms with 18F-NaF and 124I-IAZG radiopharmaceuticals distribution as input, the efficiency gain for a complete
MicroPET simulation is about 3-fold for 18F-NaF and 4-fold for 124I-IAZG as compared to the regular
navigation method, the fastest method currently available in GATE. The acceleration of GATE/MPHG is even more significant when compared with compressed voxel and parameterized tracking methods, especially for 124I. The computation speedup makes GATE/MPHG method very useful for tasks
requiring large object simulations, such as realistic simulations of biological distributions, the validation of reconstruction algorithms, and the studies of cardiac or respiratorymotions in dynamic PET.
When compared to other tracking methods for voxel phantoms, the main benefits of this method are its efficiency in photon tracking and the simplification of decay processes while still retaining the most imaging characteristics of the radionuclide. Although the underlying SimSET/PHG behind MPHG can only simulate the photon tracking in voxelized phantom as well as the modeling of positron range. However, in nuclear medicine, most imaging instruments, such as PET, SPECT and gamma camera are employed to capture the gamma ray. Hence, GATE/MPHG can meet most needs for simulations of imaging applications in nuclear medicine. On the other hand, when charge particle simulations are required, e.g. for dosimetry applications in radionuclide therapy (Chuang et al., 2014), other tracking methods with careful verification of the transportation step size should be considered, such as compressed voxel or parameterized tracking (Taschereau and Chatziioannou, 2008).
Although we only demonstrated the feasibility of the proposed method with MicroPET imaging of
18F and 124I, the method is versatile and can also be applied to other radioisotopes for diagnostic
nuclear medicine by inputting a user-defined decay scheme. We preformed preliminary tests, not detailed above, of an 111In line source inside HDPE phantoms (5 cm radius and 15 cm length) with the
SPECT detector built in SPECT benchmark of GATE. The MPHG method achieved an approximately 5.3-fold speedup when compared with the regular tracking method. In addition, the detection probability for the cascaded gammas with higher energy or lower abundance within the energy window is usually small. Hence, it’s possible to further optimize the simulated cascaded gamma lines to accelerate the simulation by only modeling the major cascaded gammas.
Despite the lack of time modeling in native mode of SimSET/PHG, SimSET can order and stamp the list mode file with their detection times through SimSET’s post-processing utilities (Harrison et
al., 2005). This function required huge storages of the list mode files for reassigning the decay and
our developed program, as one of our goals was to avoid list mode files. By contrast, the most distinctive feature of GATE apart from other MC codes is the modeling of time information throughout the simulation. But this information also makes it difficult to apply the variance reduction techniques for time-dependent phenomena such as time-of-flight, detector dead time, random coincidence rate and time-based coincidence detection. As a result, most variance reduction techniques were developed for SPECT (De Beenhouwer et al., 2008; Descourt et al., 2010), while few variance reduction techniques have been introduced for PET in GATE. Therefore it is important to improve the efficiency for voxelized phantom in PET, which we have done here using MPHG. On the other hand, as SimSET/PHG is equipped with variance reduction techniques, further work could investigate applying these variance reduction techniques in GATE/MPHG.
We have established a new MC simulation tool to improve the simulations of voxelized phantom in GATE by embedding our MPHG program, based on SimSET/PHG. The development should be of interest to the user communities of both SimSET and GATE for studying clinical and preclinical imaging. We also plan to release the source code to the communities as open source software in the future.
Acknowledgements
The authors would like to thank Cynthia Chuang for reading this manuscript and to Robert Harrison, University of Washington Medical Center, for help with the SimSET package. They also thank Dr. Uwe Pietrzyk and Michaela Gaens, Institute of Neurosciences and Medicine (INM), Research Center Juelich, Germany, for their helps with ion source management in GATE. Our thanks to the National Science Council of Taiwan, for financially supporting this research under Contract No. NSC102-NU-E-007-007-NU and NSC101-2211-E-007-012-MY3.
References
250-303
Allison J et al 2006 Geant4 developments and applications IEEE Trans. Nucl. Sci. 53 270-8
Arce P, Apostolakis J and Cosmo G 2008 A technique for optimised navigation in regular geometries.
IEEE Nuclear Science Symp. Conf. Record pp 857-9
Barret O, Carpenter T A, Clark J C, Ansorge R E and Fryer T D 2005 Monte Carlo simulation and scatter correction of the GE Advance PET scanner with SimSET and Geant4 Phys. Med. Biol.
50 4823-40
Bhat M 1992 Evaluated nuclear structure data file (ENSDF). In: Nuclear Data for Science and
Technology: Springer pp 817-21
Buvat I and Castiglioni I 2002 Monte Carlo simulations in SPET and PET QJ Nucl. Med 46 48-61 Buvat I and Lazaro D 2006 Monte Carlo simulations in emission tomography and GATE: An
overview Nucl. Instrum. Methods Phys. Res. A 569 323-9
Chen C L, Wang Y, Lee J J S and Tsui B M W 2008 Integration of SimSET photon history generator in GATE for efficient Monte Carlo simulations of pinhole SPECT Med. Phys. 35 3278-84 Chuang K S, Lu J C, Lin H H, Dong S L, Yang H J, Shih C T, Lin C H, Yao W J, Ni Y C, Jan M L
and Chang S J 2014 Improvements on a patient-specific dose estimation system in nuclear medicine examination Radiat. Prot. Dosimetry 158 1-7
Cullen D E, Hubbell J H and Kissel L 1997 EPDL97: The evaluated photon data library,’97 version
UCRL-50400 6 1997
Daube-Witherspoon M E and Muehllehner G 1987 Treatment of axial data in three-dimensional PET
J. Nucl. Med. 28 1717–24
De Beenhouwer J, Staelens S, Vandenberghe S and Lemahieu I 2008 Acceleration of GATE SPECT simulations Med. Phys. 35 1476-85
De Beenhouwer J, Staelens S, Vandenberghe S, Verhaeghe J, Van Holen R, Rault E and Lemahieu I 2009 Physics process level discrimination of detections for GATE: Assessment of
contamination in SPECT and spurious activity in PET Med. Phys. 36 1053-60
Descourt P, Carlier T, Du Y, Song X, Buvat I, Frey E, Bardies M, Tsui B and Visvikis D 2010 Implementation of angular response function modeling in SPECT simulations with GATE
Phys. Med. Biol.55 N253-66
Harrison R L, Alessio A M, Kinahan P E and Lewellen T K 2004 Signal to noise ratio in simulations of time-of-flight positron emission tomography. IEEE Nuclear Science Symp. Conf. Record pp 4080-3
Harrison R L, Gillispie S B, Alessio A M, Kinahan P E and Lewellen T K 2005 The effects of object size, attenuation, scatter, and random coincidences on signal to noise ratio in simulations of time-of-flight positron emission tomography. IEEE Nuclear Science Symp. Conf. Record pp 1900-4
Harrison R L, Kaplan M S, Vannoy S D and Lewellen T K 1999 Positron range and coincidence non-collinearity in SimSET. IEEE Nuclear Science Symp. Conf. Record pp 1265-8
Preliminary experience with the photon history generator module of a public-domain simulation system for emission tomography. IEEE Nuclear Science Symp. Conf. Record pp 1154-8
Haynor D R, Harrison R L, Lewellen T K, Bice A N, Anson C P, Gillispie S B, Miyaoka R S, Pollard K R and Zhu J B 1990 Improving the efficiency of emission tomography simulations using variance reduction techniques IEEE Trans. Nucl. Sci. 37 749-53
Haynor D R, Harrison R L and Lewellen T K 1991 The use of importance sampling techniques to improve the efficiency of photon tracking in emission tomography simulations Med. Phys. 18 990-1001
Heitler W 1954 The quantum theory of radiation (London: Oxford University Press)
Jan S et al 2004 GATE: a simulation toolkit for PET and SPECTPhys. Med. Biol. 49 4543-61
Kaplan M S, Harrison R L and Vannoy S D 1997 Coherent scatter implementation for SimSET. IEEE
Nuclear Science Symp. Conf. Record pp 1303-7
Knoess C, Siegel S, Smith A, Newport D, Richerzhagen N, Winkeler A, Jacobs A, Goble R N, Graf R and Wienhard K 2003 Performance evaluation of the microPET R4 PET scanner for rodents
Eur. J. Nucl. Med. Mol. Imaging 30 737-47
Laedermann J-P and Décombaz M 2000 Simulation of nuclear decay Appl. Radiat. Isot. 52 419-25 Lartizien C, Kuntner C, Goertzen A, Evans A and Reilhac A 2007 Validation of PET-SORTEO Monte
Carlo simulations for the geometries of the MicroPET R4 and Focus 220 PET scanners Phys.
Med. Biol. 52 4845-62
Lehnert W, Gregoire M-C, Reilhac A and Meikle S R 2011 Analytical positron range modelling in heterogeneous media for PET Monte Carlo simulation Phys. Med. Biol. 56 3313-35 Mok G S, Du Y, Wang Y, Frey E C and Tsui B M 2010 Development and Validation of a Monte
Carlo Simulation Tool for Multi-Pinhole SPECT Mol. Imaging Biol. 12 295-304
NEMA 2008 Performance measurements for small animal positron emission tomographs NEMA Standards Publication NU 4–2008, Technical Report (Rosslyn, VA: National Electrical Manufacturers Association)
OpenGATE collaboration 2014 http://www.opengatecollaboration.org/
Palmer M R and Brownell G L 1992 Annihilation density distribution calculations for medically important positron emitters IEEE Trans. Med. Imaging 11 373-8
Picard Y, Thompson C J and Marrett S 1992 Improving the precision and accuracy of Monte Carlo simulation in positron emission tomography IEEE Trans. Nucl. Sci. 39 1111-6
Pietrzyk U, Zakhnini A, Axer M, Sauerzapf S, Benoit D and Gaens M 2012 EduGATE–basic examples for educative purpose using the GATE simulation platform Z Med Phys.23 65-70
Rehfeld N S, Stute S, Apostolakis J, Soret M and Buvat I 2009 Introducing improved voxel navigation and fictitious interaction tracking in GATE for enhanced efficiency Phys. Med. Biol. 54 2163-78
Schümann J, Paganetti H, Shin J, Faddegon B and Perl J 2012 Efficient voxel navigation for proton therapy dose calculation in TOPAS and Geant4 Phys. Med. Biol. 57 3281-93
Segars W P, Tsui B M, Frey E C, Johnson G A and Berr S S 2004 Development of a 4-D digital mouse phantom for molecular imaging research Mol. Imaging Biol. 6 149-59
Tang J, Rahmim A, Lautamäki R, Lodge M A, Bengel F M and Tsui B M 2009 Optimization of Rb-82 PET acquisition and reconstruction protocols for myocardial perfusion defect detection Phys.
Med. Biol. 54 3161-71
Taschereau R and Chatziioannou A F 2008 Compressed voxels for high-resolution phantom simulations in GATE Mol. Imaging Biol. 10 40-7
Venkataramaiah P, Gopala K, Basavaraju A, Suryanarayana S and Sanjeeviah H 1985 A simple relation for the Fermi function J. Phys. G: Nucl. Phys. 11 359-364
Zaidi H 1999 Relevance of accurate Monte Carlo modeling in nuclear medical imaging Med. Phys. 26 574-608
Zaidi H 2000 Comparative evaluation of photon cross-section libraries for materials of interest in PET Monte Carlo simulations IEEE Trans. Nucl. Sci. 47 2722-35
Zanzonico P, O’Donoghue J, Chapman J D, Schneider R, Cai S, Larson S, Wen B, Chen Y, Finn R and Ruan S 2004 Iodine-124-labeled iodo-azomycin-galactoside imaging of tumor hypoxia in mice with serial microPET scanning Eur. J. Nucl. Med. Mol. Imaging 31 117-28
Zhu X and El Fakhri G 2009 Monte Carlo modeling of cascade gamma rays in 86Y PET imaging:
Table 1. Comparisons of FWHMs and FWTMs from GATE/GEANT4 and GATE/MPHG simulation results for different isotopes.
Isotope s Radial offset (mm) Radial resolution (FWHM/FWTM) Tangential resolution (FWHM/FWTM) Axial resolution (FWHM/FWTM)
GEANT4 MPHG GEANT4 MPHG GEANT4 MPHG
18F 0 1.854/4.020 1.851/4.157 1.852/4.048 1.865/4.207 1.941/4.582 1.914/4.606 10 2.135/4.352 2.138/4.461 2.211/4.427 2.209/4.573 3.029/8.123 2.991/8.196 25 3.177/6.086 3.190/6.163 2.411/5.280 2.406/5.475 3.489/9.994 3.437/9.978 124I 0 3.563/8.600 3.335/8.498 3.574/8.692 3.308/8.685 4.120/9.155 3.876/9.008 10 3.731/8.669 3.547/8.841 3.777/8.902 3.541/8.781 4.941/10.46 4.704/10.41 25 4.581/9.391 4.444/9.301 4.172/9.324 3.991/9.243 4.752/12.48 4.688/12.37
Table 2. Comparison of scatter fractions from GATE/GEANT4 and GATE/MPHG simulation results for 18F.
Phantoms Energy window (keV) Scatter fraction (%)
GEANT4 MPHG 350-600 10.16 10.09 Mouse-like 350-650 9.79 9.71 425-650 5.79 5.81 350-600 23.46 23.09 Rat-like 350-650 22.66 22.35 425-650 13.20 13.34
Table 3. Comparisons of scatter and gamma fractions from GATE/GEANT4 and GATE/MPHG simulation results for 124I
Phantoms Energy window (keV) Scatter fraction (%) Gamma fraction (%)
GEANT4 MPHG GEANT4 MPHG 350-600 10.04 9.96 19.24 19.26 Mouse-like 350-650 9.65 9.59 25.29 24.88 425-650 5.74 5.87 23.51 23.02 350-600 23.98 23.16 30.58 30.45 Rat-like 350-650 23.12 22.42 37.80 37.12 425-650 13.56 13.20 35.32 34.92
Figure 1. Full cascaded lines of decay scheme with 3 layers and the corresponding finite state machine.
Figure 2. Cross sections of photon interactions in water as calculated from EPDL97 (lines) and SimSET (symbols). The curves represent the total linear attenuation coefficient, photoelectric absorption (PE), Compton scattering (CS), and
Rayleigh scattering
(RS).Figure 3. Simulated positron spectra of 18F and 124I with (MPHG) and without (SimSET) relativistic
(a) (b) (c)
Figure 5. (a) Coronal views of the attenuation map, (b) sagittal views of the activity map for 18F-NaF
and (c) coronal views of the activity map 124I-IAZG including a spherical lesion in right hind leg,
Figure 6. Comparisons of energy spectra of detected photons of 124I line source inside a cylindrical
(a) (b)
(c) (d)
Figure 7. Comparisons of the prompt and random coincidence count rate of 18F from GATE/MPHG
and GATE/GEANT4. Upper row: 6 ns coincidence time window for (a) mouse-like and (b) rat-like phantom. Lower row: 10 ns coincidence time window for (c) mouse-like and (d) rat-like phantom.
(a) (b)
(c) (d)
Figure 8. Comparisons of the prompt and random coincidence count rate of 124I from GATE/MPHG
and GATE/GEANT4. Upper row: 6 ns coincidence time window for (a) mouse-like and (b) rat-like phantom. Lower row: 10 ns coincidence time window for (c) mouse-like and (d) rat-like phantom.
(a) (b) (c) (d)
Figure 9. Summed sinograms of the full 3D sinogram simulated with (a) MPHG, (b) Regular, (c) Compressed, (d) Parameterized methods for 18F-NaF acquisitions of MOBY phantom with voxel sizes
(a) (b) (c) (d)
Figure 10. Summed sinograms of the full 3D sinogram simulated with (a) MPHG, (b) Regular, (c) Compressed, (d) Parameterized methods for 124I-IAZG acquisitions of MOBY phantom with voxel
(a) (b)
Figure 11. Comparisons of the profiles in the sinogram summed over the projection angles (figure 9 and 10) for (a) 18F-NaF acquisitions and (b) 124I-IAZG acquisitions of MOBY phantom, respectively.
(a) (b)
Figure 12. Total CPU time needed for (a) 5 min simulations of 10 µCi 18F-NaF and (b) 10 s
simulations of 450 µCi 124I-IAZG for the four different tracking methods. The matrix dimensions were
32×32×113, 64×64×225, 128×128×450, and 256×256×900 (corresponding to voxel sizes of 1×1×1 mm3, 0.5×0.5×0.5 mm3, 0.25×0.25×0.25 mm3, and 0.125×0.125×0.125 mm3), respectively.
