Monte Carlo Simulation and a Review of the Physics of the Positron Annihilation Process in PET

In this paper, we investigate the physics of the positron annihilation process, which occurs in a PET imaging system. In particular, the diffusion of beta particles (positrons) within water was addressed. Beta particles are emitted isotropically from the same source point with random directions and randomly chosen energy levels. After traversing a certain distance within water, beta particles lose a certain amount of its energy. The inelastic collisions with atomic electrons are mainly responsible for the energy dissipation of charged particles, such as electrons and positrons (that have low mass). The energy loss rate due to inelastic process is estimated by using the Beta-Bloch formula. These results help in understanding how to develop and implement a computationally efficient Monte Carlo Simulation of the positron annihilation process.


Introduction
Positron Emission Tomography (PET) is a medical imaging technique used in clinical and preclinical investigations to produce a three dimensional image that can help to show how certain organs perform their activities (Rocha and Harbert, 1978 [25]); e.g. in clinical oncology to detect different types of cancerous tumor tissues (Strauss and Conti, 1991 [26]; Contractor et al., 2012 [8]).The radionuclides that are commonly used in positron PET, for studies and examinations performed on humans (i.e.patients) are characterized by: 1) Their short half-life time, such as Oxygen-15 (~2 min), Nitrogen-13 (~10 min), Carbon-11 (~20 min), and Fluorine-18 (~110 min), to be able to minimize the radiation dose to the patient.2) Their ideal decay characteristics.3) Their ability to be prepared with very high specific activities.4) Their chemical stability (Rocha and Harbert, 1978 [25]; Murray and Williams, 1958 [20]; Ott, 1981 [21]).
Monte Carlo simulation is an important and popular approach which plays a key role in the research field of PET imaging (Thompson et al., 1992 [27]; Buvat and Castiglioni, 2002 [3]; Buvat and Lazaro, 2006 [4]; Barba et al., 1995 [1]; Chow, 2012 [7]; Espana et al., 2009 [10]; Levin et al., 1995 [17]; Thomson and Kawrakow, 2011 [27]; Yamamoto, 2010 [29]), because of the high costs of the sophisticated hardware equipment used in the PET imaging system (especially the cyclotrons) in addition to the high costs of performing real-life experiments, including the costs related to the needed materials, laboratories and patient facilities as well as all the staff who should be involved in the different phases of the experimental tasks (performed on humans or phantoms).
According to the number of publications obtained from Google Scholar for the time period beginning from 1990 and ending in 2011 (between 1990 and 2011), the interest in using the Monte Carlo method for the simulation of the PET imaging system increased gradually during this time period, as illustrated in Table 1.
Monte Carlo Simulation became a very powerful tool that is used to achieve better understanding of the impact of each part of the PET imaging system and to gain detailed information about the performance of the current design of the simulated PET imaging system and how to improve the currently simulated design further.And this is the main advantage of using Monte Carlo Simulation, which makes it easy to optimize the different system parts and parameters with high flexibility.
However, the latter issue also has some disadvantages, mainly due to the high complexity of an actual PET imaging system.The consequence of this issue is to try to simplify and include only a subset of the parameters of the actual system in the simulation to achieve a result which is as close as possible to the actual system.Usually, the most important parameters, which have significant   1990 1991 1992 1993 1994 1995 1996 1997 1998 1999 2000 2001 2002 2003 2004 2005 2006 2007 2008  effects on the performance and efficiency of the resulting PET imaging system, are taken into account.This paper addresses one of the most important parameters of the PET imaging system.Currently, the positron annihilation process is still playing the most significant role in limiting the spatial resolution of the PET imaging system (Cal-Gonzalez et al., 2009 [5]; Levin and Hoffman, 1999 [18]; Jodal et al., 2012 [13]; Paganetti, 2012 [22]; Peng and Levin, 2012 [23]).This process begins with the decay process of positron energy which ends with the positron-electron annihilation process, which in turn results in producing and emitting two 511-keV gamma photons (in opposite directions) that can easily penetrate the body of the patient.Each such pair of gamma photons are detected by the sensors of the PET camera which surrounds certain part of the patient's body where the annihilation process occurred (Rocha and Harbert, 1978 [25]).Therefore, it is important to understand the physics behind these processes to be able to find solutions to this open problem within the field of PET imaging.
In this paper, the physical laws behind the annihilation process are investigated and analyzed to be able to understand how to perform a Monte Carlo simulation of this process in a computationally efficient way.

Radiative Energy Dissipation Effects on Positron Trajectory Calculation
Two main types of interactions exist.The first one is inelastic collisions with atomic electrons of the material.The other one is elastic scattering from nuclei in without any radiation or excitation.For the range of energies of beta (β) particles of our interest (between 10 4 -10 7 eV), almost all the betas are deflected because of elastic scattering from nuclei (Evan, 1955 [11]; Meyerhof, 1967 [19]).Additionally, there exist less important interactions which are called Emission of Cherenkov Radiation, Nuclear Reactions and Bremsstrahlung.In practice, we can ignore the Bremsstrahlung because the resulting energies are out of the range of interest.The Cherenkov Radiation is emitted only if the speed v of a charged particle, such as an electron or a positron which traverses a dielectric medium, is greater than the speed of light c.In our study the speed of the β-particle is usually not greater than c, therefore the contribution of the Cherenkov Radiation is negligible.
There are two types of elastic atomic collisions: soft and hard collisions.In a soft collision corresponding to a not-so-close encounter ( ) b a  , which is responsible for only the excitation of an atom, the incoming particle's trajectory almost never deviate from its original track.The energy transferred by an excitation is dissipated in the emission of low-energy radiation of frequencies of 10 18 -10 21 Hz as well as in atomic and molecular vibrations.While in a hard collision, the incoming particle's trajectory is changed considerably.When this kind of process takes place, there occurs a large energy transfer to an atomic electron which is called a recoiling ionization electron.If the transferred energy is enough to produce its own secondary ionization, this recoiling electron is called a knock-on ora delta electron (δ-electron).

Energy Dissipation by Inelastic Collision
The inelastic collisions with atomic electrons are mainly responsible for the energy dissipation of charged particles of low mass, such as electrons and positrons.As a traversing particle moves from one point to another in matter, its energy will be decreased through different interactions with the atomic electrons that it encounters.The energy transferred from the particle to the atomic electron will cause ionization and excitation.The rate, at which this energy loss occurs, depends on several parameters, such as the type and energy of these particles and the atomic-number density of the medium through which the particle is traversing.However, at a high primary energy level, excitation is much more probable to occur than ionization.Ionization and excitation processes can occur in all atomic shells, but mostly in the inner ones.If an inner-shell electron is ejected, then filling this vacancy will be accompanied by an emission of an electron from this atom.
Usually, at most 50% of the kinetic energy of the primary particle can be transferred to a free atomic electron.The average energy lost per distance unit traversed by a β-particle is called the stopping power.This energy loss rate due to inelastic processes, denoted by d d E x (in-keV/µm), is estimated according to the Bethe-Bloch formula for betas (Leo, 1987 [16]; Knoll, 1989 [14]) as follows: ( ) with  where 0 r is the classical electron radius (2.82 × 10 −13 cm), ( ) , N is the atomic-number density for the medium, Z is the atomic number, and I is the average excitation potantial of the medium in eV ( according to Leo (1987) [16].The average energy loss S is valid for both electrons and positrons.It can be noticed from the formulae that A and B change slowly while S increases rapidly as the particle slows down.

Hard Collisions and Knock-On (δ) Electrons Production
While a β -particle with kinetic energy T E passes through a medium of atomic-number density e N , hard electron-positron collisions occur leading to a δ-electron ray of energy E δ .The number of these δ-electrons created per unit distance is [24]: T W E is the incoming particle energy and is given by: ( ) where h (Planck's constant) = c (speed of light) = 1.This formulais acceptable for T E E δ  .δ-electronsare produced at energies higher than the cut-off energy level of 50 keV.

Multiple Coulomb Elastic Scattering from Nuclei and Moliere's Theory
With respect to the multiple Moliere's theory, the probability of having an electron of momentum p and velocity v, scattered at angle θ and the angular interval dθ after traversing a thickness t within a material of atomic number Z and atomic-number density N, can be written as: ( ) ( ) where y is a dummy variable, 0 J is the zeroth-order Bessel function, where c χ characterizes the minumum single scattering angle which is also called the angle parameter and is given by: ( ) The screening angle a χ ′ describes the scattering by a single parameter: χ is given by: where λ′ is the electron DeBroglie wavelength, 0 a is the Bohr radius, h p λ′ = is Planck's constant and For instance, in a waterequivalent tissue at 25˚C, a β-particle of 80 keV energy can traverse through approximately 0.14 mm within the tissue where nearly 30 collisions occur and B ≈ 5; Hence, Moliere's theory can be used here, but not in the case of a 50 keV β-particle with a traversing distance of 0.06 mm in water, where around 20 collisions might occur (Levin and Hoffman, 1999 [18]).For all angles, ( ) f θ can be expanded as a power se- ries in 1 B − by variable-change to ( ) in which where 1 2 .u B y = If the scattering angle is much bigger than 0 χ , the distribution function ( ) f θ turns into the Rutherford single scattering law; i.e., ( ) ( ) Hence, the ratio of the total probability to the Rutherford scattering probability becomes: where (0) f , (1)  f and (2)  f are given as follows: , and 16 ln 0.4 In practice, it is enough to include the terms of ( ) in the calculations to get an accurate result.

Initial-Energy Spectrum of Emitted Positron
This energy spectrum is required to be able to determine the β-particles trajectories for a certain isotope.Fermi revealed the expressions for the β-energy spectra as well as for the half-life of β-decay (Fermi, 1934 [12]; Evans, 1955 [11]; Konopinski 1966 [15]; Wu and Moskowski, 1966 [28]; Cengiz and Almaz, 2004 [6]).For instance, the end point or maximum energy for 18 F is 635 KeV while its mean energy is 250 KeV and its is 110 min.All isotopes are allowed or super-allowed β-decays and their energy spectra distributions are given by: ( ) ( ) where ( ) , F Z E is the Fermi function of the energy E of the emitted positron.Z is the atomic number of the beta decay daughter.N(E) is the number of decays at an energy between E and E + dE, g is a coupling constant, E and E max are the total and the maximum (end point) βenergies, respectively (in units of mc 2 ), p is the momentum of β-particle (in units of mc) (Wu and Moskowski, 1966 [28]; Daniel, 1968 [9]).

Monte Carlo Simulation
The Monte Carlo simulation begins with simulating the medium where the positrons are emitted with energies randomly chosen according to the distribution described previously for the selected isotope.
To achieve a sufficient simulation accuracy of the positron annihilation process, the following three types of processes (which affect the energies and directions of the emitted positrons) should be taken into account: 1) Energy loss by inelastic collisions (i.e. by ionization and excitation processes), as described in subsection (II.B).2) Hard collisions and δ-electrons production, as described in subsection (II.C). 3) Elastic scattering by nuclei, as discussed in subsection (II.D).
Figure 1 shows a 2D-projectionof the simulation results of 100 tracksfor 18 F positronsgenerated from a point source in water.These results correspond to those obtained by Levin and Hoffman (1999) [18].
The simulation algorithm was implemented in Mat lab ® .The needed input parameters were selecting the isotope and the medium in addition to defining the 3D-distribution of this isotope.The algorithm proposed and implemented in this work performs the computations of the positron energy loos in small steps (taking into account the processes discussed previously).Furthermore, our algorithm has the new feature or functionality that it generates, at the beginning of each step, a little region or neighborhood around the current position of the positron.This feature leads to the advantage that it becomes easy to simulate what happens within that region with high spatial resolution.

Conclusions
A user friendly implementation of the simulation algorithm for the positron annihilation process was carried out.The obtained results were satisfactory and of sufficient accuracy.From Figure 1, it becomes very clear how the positron annihilation process degrades the spatial resolution of the PET imaging system, because each point source or voxel of the imaged target object is converted into a much larger 3D-spot or clouds.This leads to a considerably blurred image, because the voxels (each of which have a certain isotope concentration) of the imaged target object are positioned beside each other, while what the PET camera sensors detect are 3D-spots or clouds (each of which is centered at one of these voxels) that overlap each other extensively.
Hence, the next future step for this research work is to use this efficient simulation tool to study the positron annihilation process in order to understand and find out how to reduce the positron range to be able to enhance the spatial resolution of the PET imaging system.
It is of course also possible to study the other processes and parameters that affect the spatial resolution of the PET imaging system, such as the effect of the area of the sensing head of the detector, and the effect of the non-collinearity of the annihilation photons pair, which means that the angle between these photons is not always 180˚.

χ
is the critical scattering angle, below which, deviations from the Rutherford scattering law occur.0 Moliere's theory is only valid for 0 Ω20 > where the constant value B (defined as: b B lnB = − ) > 4.5.

Figure 1 .
Figure 1.A 2D-projection of the simulation results of 100 18 F positron tracks emitted from a point source centered in a little water volume.