Vol 28, No 1 (2023)
Review paper
Published online: 2023-02-13

open access

Page views 2395
Article views/downloads 409
Get Citation

Connect on Social Media

Connect on Social Media

The effect of voxelization in Monte Carlo simulation to validate Bragg peak characteristics for a pencil proton beam

Gumec M. Gungor Price1, Neslihan Sarigul2
Rep Pract Oncol Radiother 2023;28(1):102-113.

Abstract

Background: The purpose of this research was to show how the Bragg peak (BP) characteristics were affected by changing the voxel size in longitudinal and transverse directions in Monte Carlo (MC) simulations by using Geant4 and to calculate BP characteristics accurately by considering the voxel size effect for 68 MeV and 235.81 MeV.

Material and methods: Different interpolation techniques were applied to simulation data to find the closest results to the experimental data.

Results: When the x-size of the voxel was increased 2 times at low energy, the maximum dose increase in the entrance and plateau regions were 17.8% and 17%, respectively, while BP curve shifted to the shallower region, resulting in a 0.5 mm reduction in the curable tumor width (W80pd). At high energy, the maximum dose increase at the entrance and plateau regions were 0.4% and 0.6%, respectively, while it was observed that W80pd did not change. When the y-z sizes of the voxel were increased 2 times at low energy, the maximum dose reduction at the entrance and plateau regions was 3.4%, but no change was observed in W80pd. At high energy, when the y-z sizes of the voxel were increased 2.2 times, the maximum dose reduction at the entrance and plateau regions were 8.9% and 9.1%, respectively, while W80pd increased by 0.5 mm. When linear, cubic spline, and Akima interpolations were applied to the simulation data, it was found that the results closest to the experimental data were obtained for Akima interpolations for both energies.

Conclusion: it has been shown that the voxel size effect for the longitudinal direction was more effective at low energy than at high energy. However, the voxel size effect for the transverse direction was more effective for high energy.

review article

Reports of Practical Oncology and Radiotherapy

2023, Volume 28, Number 1, pages: 102–113

DOI: 10.5603/RPOR.a2023.0007

Submitted: 07.06.2022

Accepted: 06.02.2023

© 2023 Greater Poland Cancer Centre.

Published by Via Medica.

All rights reserved.

e-ISSN 2083–4640

ISSN 1507–1367

The effect of voxelization in Monte Carlo simulation to validate Bragg peak characteristics for a pencil proton beam

Gumec M. Gungor Price1Neslihan Sarigul2
1Arts-Sciences Faculty, Physics Department, Cukurova University, Saricam, Adana, Türkiye
2Institute of Nuclear Science, Hacettepe University, Ankara Türkiye

Address for correspondence: Gumec M Gungor Price, MSc, Arts-Sciences Faculty, Physics Department, Cukurova University, 01330 Saricam, Adana, Türkiye, tel: +44 7913 166895; e-mail: petekgungor@gmail.com

This article is available in open access under Creative Common Attribution-Non-Commercial-No Derivatives 4.0 International (CC BY-NC-ND 4.0) license, allowing to download articles and share them with others as long as they credit the authors and the publisher, but without permission to change them in any way or use them commercially

Abstract
Background: The purpose of this research was to show how the Bragg peak (BP) characteristics were affected by changing the voxel size in longitudinal and transverse directions in Monte Carlo (MC) simulations by using Geant4 and to calculate BP characteristics accurately by considering the voxel size effect for 68 MeV and 235.81 MeV.
Materials and methods: Different interpolation techniques were applied to simulation data to find the closest results to the experimental data.
Results: When the x-size of the voxel was increased 2 times at low energy, the maximum dose increase in the entrance and plateau regions were 17.8% and 17%, respectively, while BP curve shifted to the shallower region, resulting in a 0.5 mm reduction in the curable tumor width (W80pd). At high energy, the maximum dose increase at the entrance and plateau regions were 0.4% and 0.6%, respectively, while it was observed that W80pd did not change. When the y-z sizes of the voxel were increased 2 times at low energy, the maximum dose reduction at the entrance and plateau regions was 3.4%, but no change was observed in W80pd. At high energy, when the y-z sizes of the voxel were increased 2.2 times, the maximum dose reduction at the entrance and plateau regions were 8.9% and 9.1%, respectively, while W80pd increased by 0.5 mm. When linear, cubic spline, and Akima interpolations were applied to the simulation data, it was found that the results closest to the experimental data were obtained for Akima interpolations for both energies.
Conclusion: it has been shown that the voxel size effect for the longitudinal direction was more effective at low energy than at high energy. However, the voxel size effect for the transverse direction was more effective for high energy.
Key words: voxel size effect; proton therapy; Bragg peak; Monte Carlo simulations; Geant4
Rep Pract Oncol Radiother 2023;28(1):102–113

Introduction

During the past decade advancements in proton therapy have grown considerably, especially for cancer treatment [1, 2]. The shape of the proton depth dose distribution, with its low entrance dose and well-defined range ending with the high dose Bragg peak (BP), is without doubt critical to these advancements [3]. Protons can deliver a high dose to the target without inflicting any surrounding normal tissue [4], so this ability allows proton therapy to spare healthy tissue and have a conformal dose distribution [5–7]. Use of protons in radiation therapy was first treated in 1954 at the Lawrence Berkeley Laboratory [8]. Proton therapy centres are now appearing more widespread across the world and treatment modality is becoming very well known.

Different proton therapy centres substantially vary from one to another, which affects the beam line design and the primary part of the depth dose characteristics of the proton beam [9]. Determination of characteristics of all radiation beams is vital. In proton therapy, the quality of a proton beam is mostly expressed in terms of its physical properties; the deposition of the major proportion of the dose within a few millimetres in a well-defined depth (BP), followed by a sharp dose fall-off and a subsequent negligible dose deposition thereafter, but one of the biggest factors affecting BP dosimetric characteristics is a geometrical collection efficiency [10], which can be defined as the partial volume effect associated with the process of voxelization [11]. BP is defined by its characteristics. Percentage dose ranges, such as proximal depth of the 50% dose (R50p) and distal depth of the 50% dose (R50d), proximal depth of the 80% dose (R80p), distal depth of the 80% dose (R80d), distal depth of the 90% dose (R90d), depth of the peak (Rpeak), distal depth of the 20% dose (R20d), a width of 80–20% distal fall-off (W80–20dd), full with at half maximum (FWHM) of BP, width between proximal and distal depth of the 80% dose (W80pd) and peak to entrance dose ratio (a division of maximum dose of the BP by entrance dose), are parameters for depth dose characteristics [2].

Monte Carlo (MC) dose calculation is considered to be a gold standard for simulating doses delivered in radiation therapy. MC simulations have increasingly been used for dose calculations in proton therapy due to inherent accuracy. Several MC transport codes have been published for proton radiotherapy [12, 13]. MC simulations, based on the particles using theoretical models and experimental data section, take into account the physics of particle interaction particle by particle. In photon therapy, the difference between dose distributions obtained with pencil beam (PB) algorithms and MC generated dose distributions can be significant for certain treatment areas and beam configurations. This may be even more significant in proton therapy because of steep dose gradients, which are often near critical structures. MC simulation has a much larger clinical impact in proton therapy than in photon therapy [14]. The main reason for the increased interest in proton therapy is the physical characteristics of the depth dose curve with a dose peak (BP) at a well-defined depth in tissue.

Different MC codes, such as GEANT4 [15] (Geometry and Tracking), FLUCA [16], PENELOPE [17], and MCNPX [18] have been used by several authors for proton therapy. GEANT4 is a MC particle physics simulation toolkit, implemented in object-oriented C++. It is very important to find an optimum value of ionization potential energy and to define voxel dimensions for Geant4 MC toolkit at energies of interest in particles. The proton range of energetic ions depends on the mean ionization potential (I-value) of the stopping medium which has a strong influence on the BP spatial position [19]. Uncertainties in the mean ionization potential energy are presumably significant. International Committee for Radiological Units (ICRU) Report 49 [20] recommends mean ionization potential energy, or I-value of water as 75 ± 3 eV. Another study has reported 80 ± 2 eV [21]. The stopping power tables by Janni [22] are based on a value of 81.8 eV while a more recent experimental study resulted in a value of 78.4 ± 1.0 eV [23]. Andreo [24] has demonstrated that a variation of I-values between 67.2 eV and 80 eV can have a substantial impact on the proton beam range in water. The value I = 75 ± 3 eV for water from ICRU Reports 37 [25] and 49 [20] has been accepted as the norm for many years. Then, a lower value I = 67.2 eV has been assumed in ICRU Report 73 [26] for water. Increased I values for water have also been suggested by Kramer (I = 77 eV) [27], Faddegon (I = 80–83 eV) [28], Dingfelder [29] (I = 81.8 eV), and Emfietzoglou et al. [30] (I = 80–85 eV) for different proton energies. Researchers have also suggested different I of water in literature [24, 31, 32]. The I value of water is a subject of growing interest and values between 67.2 eV and 85 eV were reported [21–24, 33, 34]. It was concluded that the range uncertainty due to uncertainties in I values is in the order of ± 1.5–2.0% [24].

The other important parameter is voxel size which, if not taken into consideration, can give erroneous results, like statistical errors, as well as prolong the computational time for commissioning. The effect of voxel size was studied for photon therapy [35–37] but to our knowledge, this effect has not been extensively studied for proton therapy. It has been reported in the literature that the voxel size effect on dose distribution in photon therapy varies between 12% and 17% [36]. The voxel size effect has been utilized to establish a relationship with both efficiency [10] and step size [38] by using different ion chambers. When calculating dose distribution, the choice of voxel size is a compromise between computational time and accuracy. Larger voxels affect the computational time and smaller voxels, should be the same size as that of the CT data about 1mm, affect the accuracy of depth dose calculation due to reduced spatial resolution [35]. Hence, simulations used for evaluation should be performed with excellent resolution to avoid any bias. While the voxel size was chosen as 0.5 × 1 × 1 mm3, 1 × 1 × 1 mm3 in some studies [39, 40], an effective radius of 4.1 cm and 6. 0 cm was used in another study [10] for energies between 90 and 226 MeV proton therapy. In the literature, possible lengths in the transverse direction of the voxel dimension for an ideal proton PB are defined [41]. The transverse length can vary between 0.85 cm and 1.35 cm at low proton energies such as 68 MeV, and between 0.75 cm and 11.47 cm at high proton energies such as 235.38 MeV [41].

The aim of this study was to investigate how BP characteristics were affected by changing the voxel size in the longitudinal and transverse directions and to calculate BP characteristics accurately by considering the voxel size effect for 68 MeV and 235.81 MeV. By applying different interpolation techniques to the simulation results, it was also aimed to investigate which of these techniques would be the most suitable for simulation data to evaluate BP characteristics.

Materials and methods

Experimental set-up

The experimental data, which was supplied by the manufacturer of VARIAN for ProBeam system (Munich, Germany), was used to choose a suitable I value of water for the proton energies used in MC simulations. Depth-dose measurements were taken for 68 MeV and 235.81 MeV proton energies at the surface of the water tank by the manufacturer. The measurement was performed with PTW 34070 Bragg Peak Chamber. Measured depth–doses were normalized to the maximum doses.

MC simulation geometry

In the following study, MC simulations were performed using Geant4 v10.01.p02 [42]. G4EmStandardPhysics_option3 standard electromagnetic package [15] was used. In this study, energy type was selected as mono [43, 44] and 68 MeV and 235.81 MeV energies were used to model a proton PB to evaluate the impact of voxelization on the characteristics of BP curve. Initial beam energy spread σE was accepted as zero [45]. Nozzle wasn’t modelled, proton beam source was placed at the entrance surface of the water tank. For production thresholds, a range cut of 0.1 mm was used for photons, electrons, and positrons. There was no range cut on protons [46]. Protons were tracked to the end of their range.

The transportation of particles in Geant4 was done in steps. One step was defined by the distance between two voxel boundaries if there was no discrete interaction in the present voxel [47]. Very small steps need to be used in order to ensure accurate simulation [48]. Range cut and maximum allowed step size should be equal to or lower than the voxel size, around 1 mm for clinical applications [49]. Maximum step size was set to 1 mm. ρwater = 0.998 g/cm3 was set in the source file according to ICRU 2016 [50].

The beam was modelled using several parameters, which were defined by beam source position and structure, beam direction, angular distribution, and energy spectra. The parameters related to the visualization of the beam were defined in the macro file. Beam direction was set from –x direction to +x direction. Protons in a beamlet incident on a patient or a phantom are very nearly monoenergetic and are distributed laterally, essentially as a narrow Gaussian function of position relative to the beamlet’s central axis [44]. Therefore, beam energy was selected as mono-energetic and it was defined in Gaussian profile. The incident surface was in the y–z plane and the beam was traveling along the x axis without any angular dispersion. The proton source was set at d = –20 cm depth on the x-axis in the y-z plane. The beam shape was set to circle. The monoenergetic proton source used is entirely a point source, therefore the beam radius was defined as r = 0 mm in the macro file. The beam thickness σr = 3 mm, which was found as a result of proton fluence measurements taken in air, was added to the macro file. 2 × 107 protons were traced to obtain less than 1% statistical accuracy. The homogenous water phantom was cubic with dimensions of 40 × 40 × 40 cm3. Simulation results were deposited in a voxelized (small cubic volumes) water phantom. Dose from all primary and secondary particles was considered. The coordinates of the modelled water phantom shown in Figure 1 were expressed in a Cartesian coordinate system.

154885.png
Figure 1. Voxelized water phantom with an incident proton beam in the yz plane. The beam was travelling along the x-axis denoted as “Depth” in the plots shown later. Colours: protons in blue, gammas in green, neutrons in yellow, electrons in red
Selection of voxel dimension

Firstly, the suitable I value of water was selected for the high and low proton energies to be used in the study. For 235.81 MeV proton energy 76 eV, 77 eV and 78 eV and for 68 MeV proton energy 85 eV, 86 eV and 87 eV mean ionization potential energies were tested. BP characteristics were measured for different voxel sizes by selecting the most appropriate I according to the results of the graphs obtained.

The water phantom was divided into voxels. In the setup geometry, x side of the voxel was varied between 0.05 cm and 0.2 cm. Voxel sizes at transverse direction (y–z size of the voxel) were varied between 1 cm to 7.4 cm to understand the contribution of Multiple Coulomb scattering to the Bragg Peak. Profiles were obtained with 0.5 mm resolution for 68 MeV and 1 mm resolution for 235.81 MeV in the beam direction (longitudinal direction). In MC simulations, plenty of voxel variations were tested and the selected ones were listed in Table 1. These voxel dimensions were chosen to test the voxel effect at 86 eV and 77 eV of I values for 68 MeV and 235.81 MeV, respectively.

Table 1. Voxel dimensions used in the study for low and high proton energies

Voxel dimensions [cm3]

Proton energy/Ionization potential energy

E = 68 MeV

I = 86 eV

E = 235.81 MeV

I = 77 eV

0.05 × 1 × 1

v

0.1 × 1 × 1

v

0.1 × 2 × 2

v

0.1 × 3.4 × 3.4

v

0.2 × 3.4 × 3.4

v

0.2 × 7.4 × 7.4

v

The characteristics of Bragg peak

In this study eleven dose range parameters were used to characterize the shape of a BP (as shown in Figure 2): Rpeak, R90d, R80p, R80d, R50p, R50d, R20d, FWHM, W80–20dd, W80pd, and peak-to-entrance dose ratio was used for comparison between the experimental and simulation data. The maximum doses in the plateau region were also compared. The plateau region is defined as the region from the entrance to the proximal depth of the 36% dose point. Characteristic parameters and regions of a BP curve are shown in Figure 2.

154897.png
Figure 2. Beam characteristics of Bragg peak (BP)
Statistical methods

All data processing, analysis, and visualization were performed using OriginPro 2017 SR2 [51]. In this study, the comparison between dose values has been done using a percentage dose normalization equation (1)

154869.png (1)

Di represents evaluated doses in measurement, Dref represents the reference doses in simulation and is the dose difference between the simulation and measurement. Doses are normalised to the maximum dose. To obtain BP characteristics, three different interpolation techniques were also used: linear (L), cubic spline (CS), and Akima (A) interpolations. All comparisons were assessed according to ICRU Report 78 [43].

Results and Discussion

In this study, the effect of voxel size on BP characteristics was examined. To calculate BP characteristics, the value at the x % dose point indicated by Rx must be known. If these points are not calculated in the simulation data set, the missing data is completed by applying the interpolation technique. L, CS, and A techniques were tested for low and high proton energies. Simulation data sets with values of 0.1 × 1 × 1 cm3 and I = 86 eV for 68 MeV and 0.2 × 3.4 × 3.4 cm3 and I = 77 eV for 235.81 MeV were selected. The difference between interpolations was examined for W80pd, which indicates the treatable tumor thickness and it was also an important parameter in clinical practice because it was known that the 80% dose range (R80d) does not change with varying initial energy spread [52]. For 68 MeV, the maximum difference between L and A, L and CS and A and CS were 5.9%, 7.7% and 2.6%, respectively. For 235.81 MeV in the same region, the difference between A and CS was less than 0.5% while the maximum difference between L-A and L-CS was 1.2%. When the experimental data and the interpolated simulation data were compared for 68 MeV for W80pd, there was a 0.1 mm difference between L and experimental data and a 0.2 mm difference between CS and experimental data while there was no difference between A and experimental data. When experimental data was compared with the interpolated simulation data for 235.81 MeV for W80pd, there was a 0.6 mm difference between L and experimental data, 0.8 mm difference between CS and experimental data and A and experimental data. When peak to entrance dose ratios was compared between experimental data and the interpolated simulation data for 235.81 MeV, 2.6% and 1.3% differences were seen between L and experimental data and CS and experimental data, respectively, while there was no difference between A and experimental data. Therefore, the A interpolation technique was used in all mathematical calculations in this study.

Selection of ionization energy

In this study, the characteristics of BP curves were obtained as a result of changing the voxel sizes used in dose deposit examined. To do this, at the first step, different mean ionization potential energies were tested for 2 therapeutic proton energies 68 MeV and 235.81 MeV by using 0.1 × 1 × 1 cm3 and 0.2 × 3.4 × 3.4 cm3 voxel dimensions, respectively. As can be seen in Figure 3, BP curves were obtained for low and high proton energies using different I values and simulation data were compared with the measurement. The I values were changed to obtain a relevant sensitivity curve for BP position.

154912.png
Figure 3. Bragg peak (BP) curves obtained using different ionisation energies for 68 MeV (A), 235.81 MeV (B)

In Figure 3A, the peak to entrance dose ratio differences between the experimental data and the simulations were 4.0%, 2% and 0% and maximum dose differences at the plateau region were 4.1%, 2.5%, and 1.6% for I of 85, 86 and 87eV, respectively. Rpeak (BP position) differences were 0.3 mm, 0.2 mm and 0.1 mm for 85, 86 and 87 eV, respectively. While there was a 0.1 mm and 1.2 mm difference in W80pd for 85 and 87 eV, respectively, there was no difference in W80pd for 86 eV. The differences in the peak to entrance dose ratio, the maximum dose at the plateau region, and Rpeak value seemed small when compared with the experimental data for I = 87 eV. Since treatable tumor thickness was an important parameter in clinical practice, I value was chosen as 86 eV for 68 MeV proton energy.

In Figure 3B, there was not any difference for the peak to entrance dose ratio between the experimental data and the simulations for I of 76 eV, 77 eV, and 78 eV. Maximum dose differences at the plateau region were 3.3%, 3.2% and 2.9% and Rpeak (BP position) differences were 0.2 mm, 0.1 mm and 0.6 mm for I of 76 eV, 77 eV and 78 eV, respectively. While there was a 0.8 mm difference in W80pd for 76 eV, a 0.6 mm difference was calculated for both 77 eV and 78 eV. When evaluating the calculation results for 235.81 MeV proton energy, the closest I value to the experimental data was chosen as 77 eV.

In the present study, water was used as a stopping medium. One of the parameters required for a precise calculation of the BP position was the I value of the stopping medium. In this study, it was observed that as I increases, the BP curve shifts to the right (deeper region) due to the relationship between range and stopping power and the result was consistent with the literature [53].

Effect of voxel size

In this section, the effect of voxel size on BP curves was tested using Geant4 MC software. First by keeping the voxel sizes constant in the y–z directions (in transverse direction), and selecting the voxel sets with different x sizes (in the same direction as the proton beam), the effect of voxel size in the longitudinal direction was examined for low and high proton energies (Fig. 4). Then by keeping the x length of the voxel constant (in longitudinal direction) and increasing the lengths in the y–z directions, the effect of voxel size in the transverse direction was examined for both energies (Fig. 5).

154923.png
Figure 4. Bragg peak (BP) curves from Monte Carlo (MC) simulation were obtained using different x dimensions for 68 MeV (A), 235.81 MeV (B)
154938.png
Figure 5. Bragg peak (BP) curves from Monte Carlo (MC) simulation were obtained using different yz sizes for 68 MeV (A), 235.81 MeV (B)

To investigate x size effect on the BP curve for 68 MeV, 0.05 × 1 × 1 cm3 and 0.1 × 1 × 1 cm3 voxel sizes were used. By keeping the y–z lengths of the voxel constant and increasing the x length of the voxel from 0.05 cm to 0.1 cm, it was observed that the dose at the entrance increased by 17.8% and the maximum dose difference at the plateau region was 17%. BP shifted to the left so range characteristics were observed in the shallower region. When it was examined in terms of characteristics, Rpeak value shifted 0.5 mm to the left. The difference for R80d was 0.2 mm and 1.4 mm for R50p. W80pd, FWHM and W80–20dd increased 0.5 mm, 1.2 mm and 0.1 mm, respectively.

For 235.81 MeV, 0.1 × 3.4 × 3.4 cm3 and 0.2 × 3.4 × 3.4 cm3 voxel sets were used to analyze x size effect on the BP curve by keeping the length of the y-z size of the voxels constant and increasing the x size from 0.1 cm to 0.2 cm. BP curves showed similar behaviour to that observed at low energy, but the differences between BP curves were not as large and sharp as they were at low energy. It was observed that the dose at the entrance increased by 0.4% and the maximum dose difference at the plateau region was 0.6%. BP shifted to the left so range characteristics were observed in the shallower region. When it was examined in terms of characteristics, Rpeak value shifted 0.6 mm to the left. The difference was 0.5 mm for R80d and 0.6 mm for R50p. W80pd, FWHM, and W80–20dd values remained same.

To investigate the y-z size effect on the BP curve for 68 MeV, 0.1 × 1 × 1 cm3 and 0.1 × 2 × 2 cm3 voxel sizes were used (Fig. 5A). When the x length of the voxel was kept constant and the lengths in the y-z directions were increased 2 times, a dose reduction was observed at the entrance and plateau region. The maximum dose difference was 3.4% at both the entrance and plateau region. It was observed that BP curve range characteristics shifted to the shallow region only less than 0.1 mm. FWHM increased 0.1 mm, while W80pd and W80-20dd remained the same.

To investigate y–z size effect on the BP curve for 235.81 MeV, 0.2 × 3.4 × 3.4 cm3 and 0.2 × 7.4 × 7.4 cm3 voxel sets were used (Fig. 5B). When the x length of the voxel was kept constant and the length in the y-z sizes were increased 2.2 times, a dose reduction was observed at the entrance and plateau region. At the entrance, the dose difference was 8.9% and the maximum dose difference was 9.1% at the plateau region. Rpeak, R80d, and R50p shifted less than 1.8 mm to the right. W80pd and FWHM increased 0.5 mm and 0.2 mm, respectively, while W80-20dd remained the same.

To perform the geometrical collection efficiency in the simulations, the voxel size should be large enough to ensure the entire beam laterally broadened by scattered and secondary contributions. The decrease in the entrance and plateau was due to primary protons undergoing nuclear interactions whereas the sharp decrease at BP is mainly due to the stopping power of primary protons. In this study, it has been seen that by increasing the y-z lengths of the voxel, BP curve was broadened by 0.5 mm for high proton energy. This will reflect the contribution from scattering and secondary particle production on the BP curve when the proton enters the phantom.

Analytical dose calculation algorithms typically project the range based on the water equivalent depth in the patient neglecting the position of inhomogeneities relative to the Bragg peak depth [54–56]. Furthermore, such algorithms are less sensitive to complex geometries and density variations, e.g. at bone-soft tissue interfaces. Consequently, analytical algorithms are not able to correctly predict the effect of range degradation caused by multiple Coulomb scattering [57–60]. In these cases, MC techniques might reduce the range uncertainty by several mm [61]. Hence, MC technique’s advantages were used in this study. When the x size of the voxel was increased, an increase in the entrance and plateau dose was observed due to the inclusion of energy deposition from non-elastic nuclear interactions, because secondaries deposit their energy further upstream and BP range characteristics have shifted to the left (shallow) region. Therefore, these interactions raise the entrance dose. While a slight increase was observed in each of W80pd, FWHM, and W80–20dd at low energy, no increase was observed in W80pd, FWHM, and W80–20dd at high energy. This could be due to increasing the voxel size and would reduce the relative uncertainty but may introduce errors due to reduced spatial resolution.

Comparison of BP characteristics with chosen voxel sizes

After investigating the effect of voxel size in longitudinal and transverse directions, simulations that were the closest to the experimental results were selected. BP characteristics obtained in the simulation results using 0.1 ×1 × 1 cm3 for 68 MeV and 0.2 × 3.4 × 3.4 cm3 for 235.81 MeV voxel sets were compared with the experimental ones and the results were given in Table 2.

Table 2. The difference between simulation results and the experimental data of BP characteristics

E = 68 MeV

E = 235.81 MeV

Simulation

Experimental

Difference

Simulation

Experimental

Difference

Peak to Ent.

4.9

4.8

1%

3.9

3.9

0.0

Dose ratio

R50p [mm]

35.2

35.1

0.1

319.4

317.2

2.2

R80p [mm]

37.4

37.6

0.2

335.8

335.0

0.8

Rpeak [mm]

38.2

38.4

0.2

340.8

340.7

0.1

R90d [mm]

38.7

38.9

0.2

343.2

343.2

0.0

R80d [mm]

38.9

39.1

0.2

344.2

344.2

0.0

R50d [mm]

39.4

39.5

0.1

346.3

346.6

0.3

R20d [mm]

39.7

40.0

0.3

348.9

349.2

0.3

W80pd [mm]

1.5

1.5

0.0

8.4

9.2

0.8

FWHM [mm]

4.2

4.4

0.2

26.9

29.4

2.5

W8020dd [mm]

0.8

0.9

0.1

4.7

5.0

0.3

As seen in Table 2, when the simulation results were compared with the experimental data for low and high energies, the difference was found to be acceptable according to the protocol [43]. In the protocol, the dose differences between experiment and simulation should be less than 2.3% and range differences should be less than 1 mm. In this study, while the difference in peak to entrance ratio was 1% for 68 MeV, there was no difference for 235.81 MeV.

For a monoenergetic proton beam, the 80% fall-off position coincides with the mean projected range of a proton, i.e. the range at which 50% of the protons have stopped. Furthermore, the 80% fall-off position is independent of the beam’s energy spread [61] and it was also demonstrated in the literature that initial energy spreads had a minimum effect on R80d so R80d can be used as a good index to represent the initial incident mean beam energy [9]. When differences in R80d (clinical range) were calculated, it was seen that while the difference was 0.2 mm for 68 MeV, there was no difference for 235.81 MeV. In proton therapy applications, another important parameter was W80pd. For this parameter, the difference was 0.8 mm for high energy while there was no difference for low energy. In this study, the only parameter that was not in good agreement with the protocol was the FWHM which was obtained for high energy. The treatment nozzle was not modelled in this study; therefore, 0–1% of the initial energy was not added into the simulations as an initial energy spread. The slight shift of FWHM was influenced by the change of R50p and R50d; in other words, the width between R50p and R50d depends on the magnitude of the added initial energy spread which varies between 0–1% of the mean initial energy [45, 61].

Conclusion

In this paper, we investigated how BP characteristics were affected by changing the voxel size in the longitudinal and transverse directions for 68 MeV and 235.81 MeV. The mean ionization potential energies that have been consistent with experimental results were determined. It has also been shown that the interpolation effect was lower at high energies, while it was higher at low energies. Additionally, BP characteristic parameters obtained from simulation for chosen voxel size (0.1 × 1 × 1 cm3 for 68 MeV and 0.2 × 3.4 × 3.4 cm3 for 235 MeV) were compared with experimental results.

It has been shown that the voxel size effect for the longitudinal direction was more effective at low energy than at high energy. However, the voxel size effect for the transverse direction was more effective for high energy. When the voxel size in the beam direction was increased at both low and high energy, the doses at the entrance and plateau regions increased and BP curves shifted towards the shallow region, but when the voxel size was increased in the transverse direction, it was observed that the dose at the entrance and plateau regions decreased in both high and low energy. The difference between BP’s obtained as a result of increasing the voxel size at low energy was quite higher compared to the one at high energy. When BP characteristic parameters obtained from simulation for a chosen voxel size were compared with experimental results, differences were in good agreement with the ICRU 78 protocol except for FWHM value. The result of this study may be helpful in improving dose calculation accuracy for high and low energy pencil proton beams by choosing suitable voxelization.

Author statement

G.G.P.: conceptualization, investigation, software, validation, formal analysis, writing original Draft. N.S.: conceptualization, formal analysis, writing original draft.

Competing interest

The authors declare that they have no known competing financial interests or personal relationships that could have influenced the work reported in this paper.

Acknowledgement

The authors would like to thank Cukurova University Rectorate for supporting this work through the PhD project FEF-2015-3886.

References

  1. Schippers JM, Lomax AJ. Emerging technologies in proton therapy. Acta Oncol. 2011; 50(6): 838–850, doi: 10.3109/0284186X.2011.582513, indexed in Pubmed: 21767183.
  2. Newhauser WD, Zhang R. The physics of proton therapy. Phys Med Biol. 2015; 60(8): R155–R209, doi: 10.1088/0031-9155/60/8/R155, indexed in Pubmed: 25803097.
  3. Pedroni E, Bacher R, Blattmann H, et al. The 200-MeV proton therapy project at the Paul Scherrer Institute: conceptual design and practical realization. Med Phys. 1995; 22(1): 37–53, doi: 10.1118/1.597522, indexed in Pubmed: 7715569.
  4. Kimstrand P, Traneus E, Ahnesjö A, et al. A beam source model for scanned proton beams. Phys Med Biol. 2007; 52(11): 3151–3168, doi: 10.1088/0031-9155/52/11/015, indexed in Pubmed: 17505095.
  5. Liu H, Chang JY. Proton therapy in clinical practice. Chin J Cancer. 2011; 30(5): 315–326, doi: 10.5732/cjc.010.10529, indexed in Pubmed: 21527064.
  6. Grevillot L, Bertrand D, Dessy F, et al. A Monte Carlo pencil beam scanning model for proton treatment plan simulation using GATE/GEANT4. Phys Med Biol. 2011; 56(16): 5203–5219, doi: 10.1088/0031-9155/56/16/008, indexed in Pubmed: 21791731.
  7. Weber DC, Wang H, Cozzi L, et al. RapidArc, intensity modulated photon and proton techniques for recurrent prostate cancer in previously irradiated patients: a treatment planning comparison study. Radiat Oncol. 2009; 4: 34, doi: 10.1186/1748-717X-4-34, indexed in Pubmed: 19740429.
  8. Tobias CA, Lawrence JH, Born JL. Pituitary Irradiation with High-Energy Proton Beams A Preliminary Report. Cancer Research 1958.
  9. Cai SY, Chao TC, Lin MJ, et al. Depth dose characteristics of proton beams within therapeutic energy range using the particle therapy simulation framework (PTSim) Monte Carlo technique. Biomed J. 2015; 38(5): 408–413, doi: 10.4103/2319-4170.167076, indexed in Pubmed: 26459793.
  10. Bäumer C, Koska B, Lambert J, et al. Evaluation of detectors for acquisition of pristine depth-dose curves in pencil beam scanning. J Appl Clin Med Phys. 2015; 16(6): 151–163, doi: 10.1120/jacmp.v16i6.5577, indexed in Pubmed: 26699567.
  11. Koketsu J, Kumada H, Takada K, et al. 3D-printable lung phantom for distal falloff verification of proton Bragg peak. J Appl Clin Med Phys. 2019; 20(9): 86–94, doi: 10.1002/acm2.12706, indexed in Pubmed: 31538716.
  12. Hérault J, Iborra N, Serrano B, et al. Spread-out Bragg peak and monitor units calculation with the Monte Carlo code MCNPX. Med Phys. 2007; 34(2): 680–688, doi: 10.1118/1.2431473, indexed in Pubmed: 17388186.
  13. Koch NC, Newhauser WD. Development and verification of an analytical algorithm to predict absorbed dose distributions in ocular proton therapy using Monte Carlo simulations. Phys Med Biol. 2010; 55(3): 833–853, doi: 10.1088/0031-9155/55/3/019, indexed in Pubmed: 20071765.
  14. Paganetti H. Monte Carlo simulations will change the way we treat patients with proton beams today. Br J Radiol. 2014; 87(1040): 20140293, doi: 10.1259/bjr.20140293, indexed in Pubmed: 24896200.
  15. Agostinelli S, Allison J, Amako K, et al. Geant4 — a simulation toolkit. Nuclear Instruments and Methods in Physics Research Section A: Accelerators, Spectrometers, Detectors and Associated Equipment. 2003; 506(3): 250–303, doi: 10.1016/s0168-9002(03)01368-8.
  16. Ferrari A, Sala PR, Fasso A, et al. FLUKA: A Multi-Particle Transport Code. 2005, doi: 10.2172/877507.
  17. Sempau J, Acosta E, Baro J, et al. An algorithm for Monte Carlo simulation of coupled electron-photon transport. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms. 1997; 132(3): 377–390, doi: 10.1016/s0168-583x(97)00414-x.
  18. Pelowitz DB. MCNPX User’s Manual. Cell. 2011(April): 1–76.
  19. Perales A, Cortés-Giraldo MA, Schardt D, et al. PO-0791: Determination of water mean ionization potential for Geant4 simulations of therapeutical ion beams. Radiother Oncol. 2017; 123: S420–S421, doi: 10.1016/s0167-8140(17)31228-8.
  20. Berger MJ, Inokuti M, Andersen HH, et al. Stopping Power and Ranges for Proton and Alpha particles-Report 49. J Int Commission Radiat Units Measur. 1993.
  21. Bischel H, Hiraoka T. Energy loss of 70 MeV protons in elements. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms. 1992; 66(3): 345–351, doi: 10.1016/0168-583x(92)95995-4.
  22. Janni J. Proton range-energy tables, 1 keV-10GeV. Atomic Data and Nuclear Data Tables. 1982: 27.
  23. Kumazaki Y, Akagi T, Yanou T, et al. Determination of the mean excitation energy of water from proton beam ranges. Radiat Measurements. 2007; 42(10): 1683–1691, doi: 10.1016/j.radmeas.2007.10.019.
  24. Andreo P. On the clinical spatial resolution achievable with protons and heavier charged particle radiotherapy beams. Phys Med Biol. 2009; 54(11): N205–N215, doi: 10.1088/0031-9155/54/11/N01, indexed in Pubmed: 19436099.
  25. International Commission on Radiation Units and Measurement. ICRU Report 37: Stopping Powers for Electrons and Positrons. ICRU 1984.
  26. Stopping of ions heavier than helium. J ICRU. 2005; 5(1): 0, doi: 10.1093/jicru/ndi001, indexed in Pubmed: 24170851.
  27. Krämer M, Jäkel O, Haberer T, et al. Treatment planning for heavy-ion radiotherapy: physical beam model and dose optimization. Phys Med Biol. 2000; 45(11): 3299–3317, doi: 10.1088/0031-9155/45/11/313, indexed in Pubmed: 11098905.
  28. Faddegon BA, Shin J, Castenada CM, et al. Experimental depth dose curves of a 67.5 MeV proton beam for benchmarking and validation of Monte Carlo simulation. Med Phys. 2015; 42(7): 4199–4210, doi: 10.1118/1.4922501, indexed in Pubmed: 26133619.
  29. Dingfelder M, Hantke D, Inokuti M, et al. Electron inelastic-scattering cross sections in liquid water. Radiat Phys Chem. 1998; 53(1): 1–18, doi: 10.1016/s0969-806x(97)00317-4.
  30. Emfietzoglou D, Pathak A, Papamichael G, et al. A study on the electronic stopping of protons in soft biological matter. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms. 2006; 242(1-2): 55–60, doi: 10.1016/j.nimb.2005.08.001.
  31. Paul H, Geithner O, Jäkel O. The Influence of Stopping Powers upon Dosimetry for Radiation Therapy with Energetic Ions. Adv Quantum Chem. 2007: 289–306, doi: 10.1016/s0065-3276(06)52013-1.
  32. Emfietzoglou D, Garcia-Molina R, Kyriakou I, et al. A dielectric response study of the electronic stopping power of liquid water for energetic protons and a new I-value for water. Phys Med Biol. 2009; 54(11): 3451–3472, doi: 10.1088/0031-9155/54/11/012, indexed in Pubmed: 19436107.
  33. Soltani-Nabipour J, Sardari D, Cata-Danil GH. Sensitivity of the Bragg peak curve to the average ionization potential of the stopping medium. Rom J Phys. 2009; 54(3–4): 321–330.
  34. Kondev FG, Dracoulis GD, Kibédi T. Atomic Data and Nuclear Data Tables. 2015: 103–104.
  35. Mora G, Pawlicki T, Maio A, et al. Effect of Voxel Size on Monte Carlo Dose Calculations for Radiotherapy Treatment Planning. Advanced Monte Carlo for Radiation Physics, Particle Transport Simulation and Applications. 2001: 549–554, doi: 10.1007/978-3-642-18211-2_87.
  36. Yani S, Dirgayussa I, Rhani M, et al. The effect of voxel size on dose distribution in Varian Clinac iX 6 MV photon beam using Monte Carlo simulation. AIP Conference Proceedings. 2015, doi: 10.1063/1.4930646.
  37. Li D, Chen B, Ran WYu, et al. Selection of voxel size and photon number in voxel-based Monte Carlo method: criteria and applications. J Biomed Opt. 2015; 20(9): 095014, doi: 10.1117/1.JBO.20.9.095014, indexed in Pubmed: 26417866.
  38. Kueng R, Frei D, Volken W, et al. Adaptive step size algorithm to increase efficiency of proton macro Monte Carlo dose calculation. Radiat Oncol. 2019; 14(1): 165, doi: 10.1186/s13014-019-1362-5, indexed in Pubmed: 31500647.
  39. Newpower M, Schuemann J, Mohan R, et al. Comparing 2 Monte Carlo Systems in Use for Proton Therapy Research. Int J Part Ther. 2019; 6(1): 18–27, doi: 10.14338/IJPT-18-00043.1, indexed in Pubmed: 31773045.
  40. Park SH, Jung WG, Suh TS, et al. Variation of Bragg Curve Characteristic Induced by Changing the Position of Inhomogeneous Material: Geant4 Simulation Study. J Korean Phys Soc. 2011; 58(2): 187–197, doi: 10.3938/jkps.58.187.
  41. Gottschalk B, Cascio EW, Daartz J, et al. On the nuclear halo of a proton pencil beam stopping in water. Phys Med Biol. 2015; 60(14): 5627–5654, doi: 10.1088/0031-9155/60/14/5627, indexed in Pubmed: 26146956.
  42. Carrier JF, Archambault L, Beaulieu L, et al. Validation of GEANT4, an object-oriented Monte Carlo toolkit, for simulations in medical physics. Med Phys. 2004; 31(3): 484–492, doi: 10.1118/1.1644532, indexed in Pubmed: 15070244.
  43. DeLuca PM, Wambersie A, Whitmore G. Prescribing, recording, and reporting proton-beam therapy. J ICRU. 2007, doi: 10.1093/jicru/ndm021.
  44. Mohan R, Grosshans D. Proton therapy - Present and future. Adv Drug Deliv Rev. 2017; 109: 26–44, doi: 10.1016/j.addr.2016.11.006, indexed in Pubmed: 27919760.
  45. Paganetti H, Jiang H, Parodi K, et al. Clinical implementation of full Monte Carlo dose calculation in proton beam therapy. Phys Med Biol. 2008; 53(17): 4825–4853, doi: 10.1088/0031-9155/53/17/023, indexed in Pubmed: 18701772.
  46. Guan F. Application of Dynamic Monte Carlo Technique in Proton Beam Radiotherapy using Geant4 Simulation Toolkit. 2012.
  47. Fippel M, Soukup M. A Monte Carlo dose calculation algorithm for proton therapy. Med Phys. 2004; 31(8): 2263–2273, doi: 10.1118/1.1769631, indexed in Pubmed: 15377093.
  48. I. VK. A Monte Carlo Package for the UCLH Proton Therapy Centre.
  49. Grevillot L, Frisson T, Zahra N, et al. Optimization of GEANT4 settings for Proton Pencil Beam Scanning simulations using GATE. Nuclear Instruments and Methods in Physics Research Section B: Beam Interactions with Materials and Atoms. 2010; 268(20): 3295–3305, doi: 10.1016/j.nimb.2010.07.011.
  50. Prescribing, Recording, and Reporting Proton-Beam Therapy. J ICRU. 2019; 7(2): 1–8, doi: 10.1093/jicru_ndm021.
  51. OriginLab Corp. https://www.originlab.com/.
  52. Bortfeld T. An analytical approximation of the Bragg curve for therapeutic proton beams. Med Phys. 1997; 24(12): 2024–2033, doi: 10.1118/1.598116, indexed in Pubmed: 9434986.
  53. Turner JE. Atoms, Radiation, and Radiation Protection. 3rd ed. Wiley‐VCH Verlag GmbH 2007 .
  54. Petti PL. Differential-pencil-beam dose calculations for charged particles. Med Phys. 1992; 19(1): 137–149, doi: 10.1118/1.596887, indexed in Pubmed: 1320182.
  55. Petti PL, Petti PL. Evaluation of a pencil-beam dose calculation technique for charged particle radiotherapy. Int J Radiat Oncol Biol Phys. 1996; 35(5): 1049–1057, doi: 10.1016/0360-3016(96)00233-7, indexed in Pubmed: 8751415.
  56. Urie M, Goitein M, Wagner M. Compensating for heterogeneities in proton radiation therapy. Phys Med Biol. 1984; 29(5): 553–566, doi: 10.1088/0031-9155/29/5/008, indexed in Pubmed: 6330772.
  57. Sawakuchi GO, Titt U, Mirkovic D, et al. Density heterogeneities and the influence of multiple Coulomb and nuclear scatterings on the Bragg peak distal edge of proton therapy beams. Phys Med Biol. 2008; 53(17): 4605–4619, doi: 10.1088/0031-9155/53/17/010, indexed in Pubmed: 18678928.
  58. Urie M, Goitein M, Holley WR, et al. Degradation of the Bragg peak due to inhomogeneities. Phys Med Biol. 1986; 31(1): 1–15, doi: 10.1088/0031-9155/31/1/001, indexed in Pubmed: 3952143.
  59. Goitein M. The measurement of tissue heterodensity to guide charged particle radiotherapy. Int J Radiat Oncol Biol Phys. 1977; 3: 27–33, doi: 10.1016/0360-3016(77)90223-1, indexed in Pubmed: 96060.
  60. Goitein M, Sisterson JM. The influence of thick inhomogeneities on charged particle beams. Radiat Res. 1978; 74(2): 217–230, indexed in Pubmed: 96480.
  61. Paganetti H. Range uncertainties in proton therapy and the role of Monte Carlo simulations. Phys Med Biol. 2012; 57(11): R99–117, doi: 10.1088/0031-9155/57/11/R99, indexed in Pubmed: 22571913.