Monte Carlo calculation of proton stopping power and ranges in water for therapeutic energies

Monte Carlo is a statistical technique for obtaining numerical solutions to physical or mathematical problems that are analytically impractical, if not impossible, to solve. For charged particle transport problems, it presents many advantages over deterministic methods since such problems require a realistic description of the problem geometry, as well as detailed tracking of every source particle. Thus, MC can be considered as a powerful alternative to the well-known Bethe-Bloche equation where an equation with various corrections is used to obtain stopping power and ranges of electrons, positrons, protons, alphas, etc. This study presents how a stochastic method such as MC can be utilized to obtain certain quantities of practical importance related to charged particle transport. Sample simulation geometries were formed for water medium where disk shaped thin detectors were employed to compute average values of absorbed dose and flux at specific distances. For each detector cell, these quantities were utilized to evaluate the values of the range and the stopping power, as well as the shape of Bragg curve, for mono-energetic point source pencil beams of protons. The results were found to be ±2% compared to the data from the NIST compilation. It is safe to conclude that this approach can be extended to determine dosimetric quantities for other media, energies and charged particle types.


Introduction
While traversing through an absorbing material, a charged particle undergoes various interactions with the orbital electrons and the nuclei of the atoms in the medium.The amount of energy lost is usually characterized by the concept of stopping power which is defined as the rate of energy loss by a charged particle per unit path length in an absorbing medium.It is denoted as -dE/dx and usually has the units of MeV cm2/g when named as mass stopping power.This concept plays an important role in radiation dosimetry and protection when determining certain properties of charged particle beams such as the range.It depends on the type and energy of the charged particle, as well as the properties of the absorbing material it penetrates, such as its density and atomic number.
Before expending its kinetic energy, a charged particle experiences many different interactions that may change its path in the form of elastic or inelastic scatterings.The kinetic energy of the particle is also decreased in the process either through transfers to electrons (collision loss) or to bremsstrahlung photons (radiative loss).The total stopping power for charged particle interactions is composed of two parts; namely collision stopping power and radiative stopping power.

−
The radiative stopping power is the result of Coulomb interaction between a charged particle and the nuclei of the absorber.Compared to collision losses with orbital electrons, this quantity is usually negligible, especially for heavy charged particles such protons, alpha particles, etc.Therefore, energy transfer from energetic charged particles to an absorber they traverse occurs mainly through Coulomb interactions of the charged particles with orbital electrons of the absorber atoms.
Theoretically, mass stopping power (in units of MeV cm2/g) of an absorber for a certain type of charged particle is given by the Bethe formula, developed in the realm of quantum mechanics and relativity.This formulation is usually refined for some density and shell corrections and is then read as where  and  are the atomic number of the absorber and the incoming charged particle, respectively;  is the atomic mass of the absorber;  : : Avagadro's number;  is the charge of an electron;  0  = is the rest energy of an electron (0.511 MeV);  is the relative velocity of the incoming charged particle;  is the mean ionization potential of the absorber [1].
Bethe's equation is very useful for calculating stopping power and ranges of charged particles in any material if one is interested only in point energies.For practical uses, there are data tables (provided by NIST) or computer codes (such as SRIM) which base their calculations on this equation [2,3].Fig. 1 shows a sample graph of mass collision stopping power of water against proton kinetic energy obtained from such data.The cumbersome nature of Bethe's equation can be avoided by use of other techniques.This study presents an alternative stochastic method to compute the stopping power and ranges of protons in water using the wellknown Monte Carlo simulations for proton energies encountered in therapeutic settings.A comparison with NIST data is also provided.

Materials and Methods
Monte Carlo is a well-established statistical technique for obtaining numerical solutions to physical or mathematical problems when physical measurements are either difficult or impossible.In charged particle transport problems, MC uses random numbers along with probability distributions to estimate the energy, position, direction and path-length of individual particles, as well as the type of physical interactions that particles experience when traversing a medium.The technique can yield quantities such as particle fluence across a surface or a cell and energy deposition in defined volumes.Thus, MC presents many advantages over deterministic methods since such problems require a realistic description of the problem geometry, as well as detailed tracking of every source particle.Thus, MC can be considered as a powerful alternative to the wellknown Bethe-Bloche equation for modeling charged particle interactions [4].
Various Monte Carlo packages are available today that are specifically developed for handling radiation transport problems.MCNPX (version 2.7.0) was used in this study to model the geometry and composition of a water cylinder.[5] This Monte Carlo code was developed at the Los Alamos National Laboratory (LANL, USA) and is available from the Oak Ridge National Laboratory (ORNL).Transport and interaction of various particles can be handled in a wide range of energies and three-dimensional arbitrary geometries.Particle interactions are treated either as continuous or discrete energies where the cross-section data are adopted from the ENDF/B libraries (Evaluated Nuclear Data Files, version B).MCNPX can facilitate a wide variety of source distributions, detector conditions and user-configurable output options, as well as providing a rich set of variance reduction techniques [6].
A water cylinder (r=30 cm; h=100 cm) was modeled in a vacuum sphere (r=150 cm).Starting from the surface of this water phantom, 1000 concentric detector disks (r=2 cm) were created to record particle data.The thicknesses of these disks were determined approximately from a fraction of proton range at the source energy.
In the MCNP model, the mono-energetic protons were emitted mono-directionally from a point source at one end of the water cylinder.Tally types F4 and +F6 were placed in every detector cell in the region of interest to record the average cell flux (in units of #/cm 2 /source particle) and the absorbed dose (in units of MeV/g/source particle) by all primary and secondary particles, respectively.
All the secondary particles that may arise from proton interactions were transported in the runs.The simulations were carried out for 1 billion proton histories that yielded better than 1% statistical errors.Overall, the precision in the results was satisfactory considering the extremely time-consuming nature of the problem.The typical CPU time was about 1000 min.

Results
MCNPX was used to compute absorbed dose and flux for protons in the energy range of 50-250 MeV at varying depths of water.For each incoming proton energy, the first detector was taken as the basis for stopping power calculations where the absorbed dose from +F6 tally (in units of //) was divided by the flux from F4 tally (in units of 1/ = /) to represent the stopping power (  = //) of these protons in water.Later, the listing of the detector cells was searched to find the depth at which the absorbed dose falls below 1% of the maximum.This depth value was taken to represent the range of the incoming proton energy.Finally, the Bragg curve was sketched for each proton energy as a graph of relative dose (%) versus absorber depth ( = /).
The calculated data of water obtained this way are given in Table 1 both for stopping power and range of 50-250 MeV protons.As expected, the stopping power of water for protons exhibits a decreasing relationship with increasing particle energy.This correlation can be approximated by a power equation of the form ( = = 0.997)  ] ^0_ (2 4 `= 206.7   Ad.e= . (3) On the other hand, the range of protons in water displays a direct proportionality with particle energy which can be summarized by this fit equation ( = = 0.999)  / = = 0.0023   @.eg (4) The results obtained in this study using MCNPX simulations were compared with NIST values and found good agreement within %4 for stopping power and %0.2 range data.Fig. 2

Conclusion
This study aims to conceptualize a Monte Carlo approach for calculating stopping power and range of charged particles.The results obtained from MCNPX simulations of 50-250 MeV protons in water medium suggest that stopping power of a mono-energetic proton beam can be calculated from the computed values of absorbed dose and cell flux.This approach is in good agreement with the tabular values provided in NIST databases based on computer codes that utilize Bethe's equation.The results are encouraging in that detailed studies that involve wider energy range, different particle types, and various materials need to be carried out to investigate the validity of this methodology.

Fig. 1 .
Fig. 1.Mass collision stopping power of liquid water for protons (adopted from NIST).

Fig. 2 .
Fig. 2. Stopping power and range of protons in water as a function of incoming proton energy

Table 1 .
is a plot of simulation results and NIST tables as a function of proton energy.Corresponding author: abozkurt@akdeniz.edu.trStopping power and range of protons in water.