A semi‐empirical model for the therapeutic range shift estimation caused by inhomogeneities in proton beam therapy

The purpose of this study was to devise a simple semi‐empirical model to estimate the range shift in clinical practices with high‐Z inhomogeneity in proton beam. A semi‐empirical model utilizing the logarithmic dependence on Z in stopping power from Bohr's classical approach has been developed to calculate the range shift due to the presence of inhomogeneity. Range shift from metallic plates of atomic number Z of various thicknesses were measured in water using a parallel plate ionization chamber and calculated with the FLUKA Monte Carlo code. The proton range shifts for bone and polymethyl methacrylate (PMMA) were estimated using the semi‐empirical model and compared with Monte Carlo calculation. The semi‐empirical equation to determine range shift and water equivalent thickness is presented. The model predicts a shift of the therapeutic range to within 2.5% accuracy for initial proton energies of 50 to 250 MeV and atomic numbers from 3.3 (effective Z for water) to 82. This equation is independent of beam energy, and thus provides range shift from high‐Z materials without the knowledge of proton energy. The proposed method of calculating the therapeutic range shift accurately requires only knowledge of the effective or actual atomic number of the inhomogeneity and the thickness of the inhomogeneity along the beam direction. The model generalizes the range shift calculation for any material based on its effective atomic number, and permits reliable prediction of the range shift for material combinations where no data is currently available. The proposed model can be readily implemented in routine clinical practice for proton range shift estimation and quality assurance on the treatment planning. PACS numbers: 87.53.‐j, 87.55.‐x, 87.53.Bn, 87.55.D‐, 87.55.Qr, 87.55.K‐

The effect of high-Z materials in proton therapy has been studied by various investigators. Urie et al. (2) presented a qualitative approach and documented the presence of dose perturbation caused by realistic inhomogeneities. Recently, Herrmann et al. (3) and Cheung et al. (4) investigated the effect of Ni-Ti and ceramic carbon-coated fiducial markers on charged particle therapy. The significant dose perturbation caused by such inhomogeneities has raised concerns over the usage of such fiducial markers close to the particle range. (5) Impact on proton transport due to the thickness of an inhomogeneity is commonly estimated by its water equivalent thickness (WET) (6) or water equivalent distance (WED). (7) Various approaches have been reported to predict the WET in particle beams. (8)(9)(10)(11)(12) A method to compensate for the inhomogeneities has also been suggested. (2) However, Gottschalk (8) commented on and cautioned concerning the use of analytical method proposed by Zhang and Newhauser (11) due to differences in mean excitation energy that could produce error in the computation of ranges. A detailed description of the excitation energy and computation of range has also been proposed by Bichsel. (13) Nichiporov et al. (5) recently presented experimental measurements of the range shift for a number of inhomogeneities at the depths corresponding to the middle of a spread out Bragg peak (SOBP). The usage of these methods for WET calculations requires a significant amount of tabulated data from the PSTAR stopping power database or the ICRU Report 49, (14,15) as well as a knowledge of the initial energy of proton beams.
Clinical practice typically operates using the concept of the required (modulated) proton range, instead of using the actual energy of the proton beam. Therefore, in order to be useful in a routine clinical practice, the technique to calculate range shift should not require knowledge of the proton beam energy. For this purpose, a semi-empirical model for the estimation of the proton range shift in the presence of inhomogeneity is developed in this study.

A. A model for range shift calculations
We start by assuming the geometry in Fig. 1, where Fig. 1(a) represents the depth-dose characteristics of proton beam in a homogenous water phantom where R w represents the range (90%) in water. Consider a slab of thickness t M (in cm) of a material, M, placed in the front of a water tank ( Fig. 1 (b)). A unidirectional proton beam with energy in the range 50 to 250 MeV is incident normally to the surface M. The amount of energy deposited in the slab corresponds to the amount of energy deposited in a slab of water equivalent thickness (WET), t W , is defined as: where ρ M is the mass density of the material, ρ W is the mass density of water, Z is the atomic number of the material M, and M W S is the mass stopping power ratio of material M to water W.
Implicitly, the mass stopping ratio depends on Z and E.
Equation (1) for the water equivalent thickness corresponds to the definition given by Zhang et al. (11,12) Figure 1(b) represents the geometry used in WET definition.
A clinical example of the geometry in Fig. 1(b) is the placement of fiducial markers on the patient's surface. A more commonly encountered situation is that of an inhomogeneity inside a patient, such as surgical clips, spinal prosthesis, metallic breast implants, dental fillings, and hip prosthesis, as illustrated by the schematic diagram in Fig. 1(c).
The proton range R in in Fig. 1(c) is defined as the range of protons in an inhomogeneous system consisting of water embedded with a layer of material M. Medin and Andreo (16) showed that the mass stopping power ratio of air to water, W air S , changes slowly as a function of depth up to the proximal slope of the Bragg peak in a homogeneous medium where the Bragg peak starts forming. The analysis of the ICRU Report 49 data (14,15) shows that the mass stopping power ratio of a material to water has a weak dependence on proton energy above 50 MeV for low-Z (for example, about 2% for Z = 13) and above 80 MeV for high-Z (about 5% for Z = 82) materials. Thus, in our semi-empirical model, the slab M is placed before the Bragg peak in water for a given energy, and the thickness of the slab t M satisfies the condition t M < RM. The observed range shift, Δx, can be defined as , where R W is the range of protons in water, and R r W is the reduced range of protons in water due to the presence of material M (see Fig. 1). The difference, , is equal to the water equivalent thickness, WET, of the material (Eq. (1)). The observed proton's range shift, and consequently the therapeutic range shift R 90 , in a water phantom embedded with a layer t M of a material M with atomic number Z is then given by: The function, WET(t M ,Z) has a monotonic dependence on atomic number Z that was fitted by a function as shown in Eq. (3a): on atomic number Z of the material.
For a particular case where ρ W = 1 g/cm 2 , Eq. 3(a) can be simplified further as: Final equations in the Results and Discussion section below will also be written in a simple form omitting the ρ W .

B. A fitting procedure
According to Bohr's classical approach for a universal form of the stopping power for ions, the stopping number L Bohr expresses logarithmic functional dependence on target atomic number. (17) A semi-empirical model with a logarithmic dependence on atomic number was used to determine the stopping powers for ions with energies of 0.1-1.0 MeV/u in elemental targets. (18) In the present work, we used the functional form A-B*Ln(Z) to represent α(Z).
The tabulated stopping power data from NIST (15) were used for fitting α(Z), with a least squares fitting routine from the Golden Software Grapher (v. 7).

c. Generalization of semi-empirical model
Equation (3) generalizes range shift for any material based on the effective atomic number. The calculations presented in Results and Discussion below uses the effective atomic number, Z eff , for protons. Please note that Z eff depends on the type of radiation and interaction probability.
A number of studies have discussed the definition of effective atomic number (19)(20)(21)(22)(23)(24)(25)(26) for photons and protons over a wide range of energies. The conventional method to determine Z eff of a material for photon transport is based on the photoelectric coefficient per electron or total photon energy absorption cross section per electron in composition, (19,22) or total photon interaction cross section per atom (21,23,24) for each element in the material. The estimated value of the Z eff varies with the photon energy due to the weighting of the photoelectric interaction process and depends on the calculation technique used. The weighted average of the number of electrons per atom gives an estimation of the average atomic number of a material consisting of multiple elements. (22) On the other hand, charged particles undergo multiple interactions during their transport through an absorbing medium. Collisions with the atomic electrons dominate in the region of the therapeutic energies. The procedure of the Z eff calculation for charged particles involves the utilization of the total stopping power. (20,24) The values of the Z eff thus calculated are different from those obtained for the conventional megavoltage photons. For instance, the classical value for Z eff for water for photons given by Johns and Cunningham (19) is 7.4, (19,22) while it is 3.3 for protons as determined Prasad et al. (20) The values for Z eff calculated by Prasad et al. is used in the Results and Discussion section.

d. Validation of the semi-empirical model
The experimental setup used to validate the model corresponds to the geometry presented in Fig. 1(b). The depth-dose data were collected for the various high-Z materials: Aluminum (Al), Titanium (Ti), Copper (Cu), Tin (Sn), and Lead (Pb) using a parallel plate ion chamber. Slabs of different high-Z materials of various thicknesses were attached to the upstream face of the water tank during the depth-dose scans. A set of measurements was taken to measure WET and range shift.
The values of the therapeutic range R 90 (defined as the 90% depth in water for a given proton energy) were determined from the depth-dose scans. The water equivalent thickness WET(t M ,Z) was derived from these data. The measured data on WET(t M ,Z) was compared to the results obtained by Eq. (3) with the use of its final form that is given by Eq. (4) in the Results and Discussion section.
The results of range shift from measurements and from the semi-empirical model were compared with the results of simulation using the general-purpose particle transport code FLUKA (27,28) version 2008.3c.0. The capability of FLUKA for proton therapy calculations has been validated in a number of studies. (29)(30)(31)(32) The model of the beam used in simulation of the range shift experiments with pure materials was verified by independent measurements on depth-dose distribution in a water equivalent phantom.
One of the FLUKA defaults, PRECISIO, was used to customize the physical model used in the simulation. The initial proton transport was simulated with a cutoff energy at 100 keV. USRBIN cards were used for scoring the particle fluence and dose with the parameter DOSE. The geometry of the irradiated phantom was defined with combinatory geometry package, which is a part of the FLUKA code. Fifty million (5*10 7 ) initial protons were used in the simulation.  size impinged normally on the surface of the phantom. The range shift is estimated as the shift of the 90% isodose line.

A. A semi-empirical equation for range shift
The function α(Z) calculated using tabulated stopping power data (15) is plotted in Fig. 3 (symbols) for various initial energies of proton beams. Since the dependence on energy is weak at E > 50 MeV for low-Z and E > 80 MeV for high-Z, a single set of fitted parameters may be used over the energy interval up to 250 MeV of the initial proton energies. The fitting parameters A and B in the functional form A-B*Ln(Z) were determined to be 1.192 and 0.158, respectively. The fitting function is presented as a solid line in Fig. 3.
A residual value is defined as the difference between the fitted value and the actual measured data values for a given Z or E value. The residual sum of squares (SSE) is the sum of the squares of all the residual values and is a measure of the quality of the fit in a least squares fitting method. The SSEs over the atomic number range from water to lead are 0.00265, 0.000205, 0.000337, 0.000839, and 0.00136 for initial proton energies of 50, 100, 150, 200, and 250 MeV, respectively. On the other hand, the variation of SSE, cumulative over the energy range from 50 to 250 MeV, with the atomic number is from 5.6*10 -5 for water to 0.0012 for lead suggesting that the model is extremely accurate.
Using Eq. (3), the final form of the empirical equation for calculation for the water equivalent thickness WET(t M ,Z) is given as: for various initial energies of protons. The symbols are the calculated values using the tabulated data on mass stopping power. (14) The curve is the fitted function for α(Z).
The observed range shift Δx in the presence of an inhomogeneity in an irradiated volume is then calculated according to Eq. (2) and (4) as: To examine the effect of heterogeneities on WET(t M ,Z ) and the corresponding range shift, the values obtained from Eq. (4) are compared with those from Eq. (1) using the tabulated data on mass stopping power. (15) The accuracy of the analytical model is shown in   (19) ) and bone (mass density 1.85 gm/cm 3 , Z eff = 4.9 (19) ) are used. Proton beam is defined by range of R 90 = 15.84 cm in water.
where the relative percent difference δ between two calculations of WET is presented. For the initial proton energy of 100 MeV, the difference between the two equations is about 1% for low-and medium-Z values (Al, Ti, Cu, Sn). It is about 1% for Tungsten (W), and 1.8% for Lead (Pb). The value of δ increases with energy to 2.2% for 180 MeV, and reaches 2.5% at 200 MeV. However δ is still significantly below 1% for low-Z materials (Al) at 200 MeV.
The WET values between the two equations in the low-energy region from 50 to 100 MeV agree to about 2% for low-Z materials and about 2.3 % for high-Z materials from 70 MeV to 100 MeV. Considering the compounds, the WET values for bone agrees with an accuracy of ~1% for all energies; however, PMMA has an error of about 3% in the considered interval of energies from 50 MeV to 250 MeV. Generally, the fitting parameters predict the WET and, consequently, the shift of the therapeutic range R 90 and 90% depth to within ± 3% accuracy for initial proton energies of 50 to 200 MeV, and atomic numbers from 3.3 (effective atomic number Z eff for water (20) ) to 82.

IV. concLuSIonS
A semi-empirical model for range shift as shown in Eqs. (4) and (5) provides a simple method to accurately estimate the range shift. The proposed method requires only knowledge on the effective or actual atomic number of the inhomogeneity, the physical density, and the thickness of the inhomogeneity along the beam direction. The results obtained from the semi-empirical model are in good agreement with measurements and Monte Carlo simulation. Relative percent difference δ between the semi-empirical fitting by Eq. (4) and calculations using the data from PSTAR (14) in Eq. (1).
The model generalizes the range shift calculation on any material based on its effective atomic number, not just the materials listed in ICRU Report 49, and permits reliable prediction of the range shift for material combinations where no data are currently available. The proposed model can be readily implemented in routine clinical practice in order to check the accuracy of a treatment plan.