Assessing the Effect of Directional Bremsstrahlung Splitting on the Output Spectra and Parameters Using BEAMnrc Monte Carlo Simulation Package

Abstract

Introduction:

EGSnrc software package is one of the computational packages for Monte Carlo simulation in radiation therapy and has several subset codes. Directional bremsstrahlung splitting (DBS) is a technique that applies braking radiations in interactions in this software. This study aimed to evaluate the effect of this technique on the simulation time, uncertainty, particle number of phase-space data, and photon beam spectrum resulting from a medical linear accelerator (LINAC).

Materials and methods:

The gantry of the accelerator, including the materials and geometries of different parts, was simulated using the BEAMnrc code (a subset code in the EGSnrc package). The phase-space data were recorded in different parts of the LINAC. The DBS values (1, 10, 100, and 1000) were changed, and their effects were evaluated on the simulation parameters and output spectra.

Results:

Increasing the DBS value from 1 to 1000 resulted in an increase in the simulation time from 1.778 to 11.310 hours, and increasing the number of particles in the phase-space plane (5 590 732-180 328 382). When the DBS had been picked up from 1 to 100, the simulation uncertainty decreased by about 1.29%. In addition, the DBS increment value from 100 to 1000 leads to an increase in uncertainty and simulation time of about 0.71% and 315%, respectively.

Conclusion:

Although using the DBS technique reduces the simulation time or uncertainty, increasing the DBS from a specific value, equal to 100 in our study, increases simulation uncertainties and times. Therefore, we propose considering a specific DBS value as we obtained for the Monte Carlo simulation of photon beams produced by linear accelerators.

Introduction

One of the standard methods for treating malignant tumors is radiotherapy with megavoltage photon beams using medical linear accelerators (LINAC). The radiotherapy treatment plans must be designed, and the dose distribution must be calculated on the patient’s body before the actual treatment. ¹ Computer software packages are usually used for designing the radiotherapy treatment plans, and dose calculations in a patient model, which is usually obtained from a CT scan performed on the patient (CT simulation). ² The dosimetry calculations in the treatment planning software (TPS) use correction or convolution-based algorithms. They have dose calculation uncertainties and errors, especially in the inhomogeneity regions or boundaries between different types of tissues like air and soft tissue. These errors are relatively small (about 2%-4%) for convolution-based algorithms; however, they can reach up to 20%. ³

Monte Carlo simulation is an accurate and efficient method for estimating radiation transport inside the accelerator and patient body.^4,5 This method results in a more precise estimate of dose distribution in patients’ bodies. Monte Carlo is used as a powerful tool to model the transport of photons in radiotherapy. ⁶ Unlike other methods, the Monte Carlo method uses more straightforward rules, identifies the history of particles, and pursues produced secondary particles. The disadvantage of this method is the long time needed for the calculation to achieve acceptable results; however, recent advances in computer processing have made this method a proper way to calculate the dose distributions in the clinic.^4,7 Notably no TPS is using this method to obtain the dose distributions due to its longer processing time than convolution-based dose calculation algorithms. Another limitation is the lack of an accurate and public model for primary electron beam dosimetry.

The Monte Carlo method can simulate the spatial and energy spectrum in different planes at different parts of the device, like inside and outside the LINAC.^8,9 The spatial and energy data of the produced particles in a plane is called phase-space, showing the simulated particles’ spatial and phase data. Implementing Monte Carlo simulation for the LINAC requires a great deal of time, and its result is a file of megabytes size or more that its storage is impractical for each treatment plan. Using phase-space data could significantly reduce the simulation time and make it possible to use it in the clinic for treatment.

EGSnrc software package is one of the essential computational software that uses the Monte Carlo method for simulating the radiotherapy and radiology devices and calculating dose distributions in different materials with different geometries like the patient models. ¹⁰ BEAMnrc code (a subset code in the EGSnrc package) has an option for decreasing simulation uncertainties or times.^11,12 Directional bremsstrahlung splitting (DBS) is a technique that aims to reduce the simulation uncertainty, where the photons are split at the time of creation. DBS technique is the one that applies braking radiation in a collision.^13,14 If a charged particle undergoes a bremsstrahlung or annihilation event, then the DBS splits the event with the bremsstrahlung splitting number (NBRSPL). The resultant photons are given weight reduced by a factor of NBRSPL⁻¹. DBS computes the pass way of these photons to determine whether they are located in the specified region of interest. If a photon is located in the region of interest, then the photon is kept and considered low-weighted; otherwise, another uncertainty reduction method called “Russian Roulette” will be applied to the photon to decide about the survival of the photon regarding the energy threshold with a random number and NBRSPL⁻¹. If the random number is smaller, then the photon is kept and its weight is multiplied by NBRSPL and becomes a high-weighted photon. The high-weighted photon will be split again when it undergoes interactions and Russian Roulette will be applied again for photons that were not located in the region of interest. As a result, DBS will eliminate many low-weight photons inside the field and a few high-weight photons outside the field. ¹⁵

The main purpose of radiotherapy is to deliver the highest dose to the tumor as well as the minimum dose to the surrounding healthy tissues. Monte Carlo methods are widely accepted as an accurate technique to calculate and verify dose distribution for radiotherapy treatments. However, the main obstacle to using Monte Carlo is the simulation time needed to gain results below a desirable level of uncertainty. Since the DBS value is one of the critical parameters affecting the simulation time, this study aimed to investigate the impact of the DBS technique on simulation time, uncertainty, particle number of phase-space data, and its effect on the photon beam spectrum. The main purpose of the current work was to obtain an appropriate DBS value for simulating 6 MV photon beams resulting from Siemens linear accelerator, based on lower simulation time and the smaller number of particle histories.

Materials and Methods

BEAMnrc software package (National Research Council, Canada), which is a subset of electron gamma shower user code (EGSnrc) that come as a package under license to the National Research Council of Canada (nrc), ¹⁶ was used in this study to simulate transmitted electrons and photons in a linear accelerator (Primus model, Siemens, Germany). BEAMDP program was used to plot the output spectrum. ¹⁷ These programs are based on the Monte Carlo simulations and were run using a desktop computer with an Intel (R) Xenon processor of 2.2 GHz speed (8 processors). The first step was simulating the medical LINAC with full geometrical details. In the next step, the obtained phase-space data according to the scoring plane defined in the BEAMnrc code were recorded for different values of DBS, and the parameter of simulation and output spectrum were compared with each other to evaluate the effect of DBS on these parameters. The details of the simulation procedure and analysis were brought in the following paragraphs.

Monte Carlo simulation with BEAMnrc code

In the initial versions of the BEAMnrc software package, there was an option to use uniform bremsstrahlung splitting (UBS). Selective bremsstrahlung splitting (SBS) was added to the code in 1998 with an improved splitting routine and simulation efficiency. ¹⁵ In 2004, a further significant improvement in efficiency from the DBS algorithm was added to the code. ¹⁷

LINAC (Primus, Siemens, Germany) with full geometrical and material details of different parts was simulated in the BEAMnrc software package. The geometrical and material information was obtained from the manufacturer data and user guide books and then inserted into BEAMnrc software for 6 MV photon beams. Different components, named component module (CM) in EGSnrc software of the LINAC gantry used for simulations are illustrated in Figure 1, including the target, collimator, flattening filter, ionization chambers, and jaws (X and Y).

Figure 1.

Left: component module (CM) used in simulations, Right: positions of scoring planes in the head of the simulated LINAC with BEAMnrc code.

The initial number of particle histories was 10⁸ particles (history), modified with DBS implementation. Simulation settings were adjusted based on Jabbari et al study. ⁶ The EXACT algorithm was used for border crossing, and the PRESTA II algorithm was used for fast electrons. The spin effects and bremsstrahlung cross-sections were enabled based on the relevant tables. ¹⁵ Simple collision and bremsstrahlung simple angle collision were considered based on the NIST model. In addition, Compton collision and incident electron ionization were considered according to the Kawarakow model and Atomic comfort in the simulation. Simple angle photoelectric and Rayleigh scattering interactions were disabled due to their small contributions to the other interactions inside the LINAC gantry.^15,18

Cut-off energies for electron and photon were chosen 0.7 MeV (ECUT), and 0.01 MeV (PCUT), respectively. These values were selected based on the maximum range of 1 mm for the electron and photon in the water. The 700ICRU material software library in PEGS4 data was used for defining the materials of different parts in the LINAC gantry. ¹⁸

The simulation parameters were entered according to the photon energy of 6 MV. We used Source # 19 with full width at half maximum (FWHM) = 0.35 cm and primary electron energy of 6.5 MeV. Actually, we applied different primary electron energies with different beam FWHM values. ⁶ For each set of primary electron energy and FWHM values, the simulated beam was irradiated to a cubic water phantom with the dimensions of 30 × 30 × 30 cm³, and with a voxel size of 3 × 3 × 3 mm³. The depth dose was calculated and compared with the measurements at the same geometry (ie, source to surface distance = 100 cm and field size = 10 × 10 cm²). The highest compliance rate among different primary electron energies and FWHM value sets with the measurements was chosen for simulation as an optimum. The optimal parameters were primary electron energy of 6.5 MeV and FWHM = 0.35 cm. Figure 2 shows the calculated depth dose simulated by optimal parameters versus the measurement data. Experimental measurements were performed using a Semiflex 3D chamber (sensitive volume of 0.07 cm³, PTW-Freiburg, Germany) and BEAMSCAN water phantom (PTW-Freiburg, Germany) with a radiation field size of 10 × 10 cm². The Kolmogorov-Smirnov statistical test was used to compare the experimental and simulation depth dose data. Higher P-values showed higher compliance between the data. Furthermore, mean differences between the experimental and simulation depth dose data were calculated in percentage. Lower differences represent higher compliance between the data. The mean percentage differences must be lower than 2% to accept the simulation parameters.

Figure 2.

Comparing the depth dose percentage curve resulted from Monte Carlo simulation and experimental data.

Evaluating the effect of changing DBS value on the simulation parameters and outputs

BEAMDP program was used to draw the spectrum for 6 MV energy with different DBS values (0, 1, 10, 100, and 1000) based on the phase-space plane’s data. Furthermore, the simulation time, particle numbers in the phase-space plane, and total reported simulation uncertainty were recorded and compared for different DBS values. The phase-space data at 3 scoring planes were recorded, and their output parameters were compared. Three scoring planes were defined below the target (plane1), after the flattening filter (plane 2), and after the jaws (below the x jaws, plane 3). These planes are shown in Figure 1.

Results

The energy spectra for different phase-space planes (scoring planes) with a DBS value of 0 are illustrated in Figure 3, and for the values of 1, 10, 100, and 1000 are shown in Figure 4. Increasing the DBS value resulted in increasing the intensity of produced particles in the energy spectrum, and lower uncertainty and errors in calculations in all scoring planes. The energy fluence distributions of the DBS value 1 and 10 had a maximum difference of 9.2% and 4.6%, respectively, with the energy spectrum of DBS = 1000. However, the energy spectrum resulting from the DBS value of 100 had only a maximum difference of 1.8% with the energy spectrum of DBS = 1000.

Figure 3.

Energy spectrums obtained with DBS = 0 from different scoring planes: (a) energy spectrum obtained from scoring plane 1, (b) energy spectrum obtained from scoring plane 2, and (c) energy spectrum obtained from scoring plane 3.

Figure 4.

Energy fluence distribution obtained with different DBS values (1, 10, 100, and 1000) from different scoring planes: (a) energy fluence distribution obtained from scoring plane 1, (b) energy fluence distribution obtained from scoring plane 2, and (c): energy fluence distribution obtained from scoring plane 3.

Table 1 represents the simulation time and the particle number in the phase-space plane below the LINAC gantry (plane 3) for different DBS values. Furthermore, the average calculated uncertainty values are presented in Table 2 from other scoring planes.

Table 1.

Simulation time (including CPU time and elapsed time) and number of particles in phase-space obtained in the scoring plane 3.

NDBS	CPU time (h)	Elapsed time (h)	Particle in phase-space
1	1.778	1.79	5 590 732
10	1.834	1.838	7 168 149
100	2.727	2.734	22 898 045
1000	11.310	11.513	180 328 382

Table 2.

Calculated simulation uncertainties in different scoring planes with different DBS values.

Average error%	1DBS%	10DBS%	100DBS%	1000DBS%	NO DBS%
Scoring plane1	2.99	2.90	2.91	3.86	5.95
Scoring plane2	6.48	5.88	5.89	6.70	9.65
Scoring plane3	8.96	7.30	5.77	6.13	16.59

According to the results, increasing the DBS values led to an increase in the simulation times and the number of particles in the phase-space plane. Additionally, increasing the DBS value resulted in a reduction in uncertainty values; however, increasing the DBS value from 100 to 1000 increased uncertainty. It can be concluded that increasing the DBS from a particular value equal to 100 in our study will cause higher simulation uncertainties.

Discussion

Monte Carlo is the most accurate method for dosimetry calculations for patients; however, simulation time is usually prolonged in Monte Carlo radiotherapy simulation, even with new fast-processing computers. Therefore, time is one of the essential parameters limiting the Monte Carlo simulation application in clinical practice.^19-21

The variance of the statistical calculations in the Monte Carlo simulations is an index related to total simulation uncertainty. To decrease the variance and consequently the simulation uncertainty to acceptable levels, the Monte Carlo simulations usually must execute many times; hence, the simulation time will be very long even with fast processing computers. Therefore, variance reduction techniques are usually applied to decrease simulation iteration (number of history) and time. In photon beam simulation, these techniques involve “splitting” bremsstrahlung interactions so that each bremsstrahlung event produces NBRSPL photons, each having weight NBRSPL⁻¹. Bremsstrahlung splitting can significantly decrease uncertainty. The amount of decrease due to bremsstrahlung splitting in uncertainty is more significant than the increase in CPU time/history required by this variance reduction technique. Therefore, the overall efficiency in photon beam simulation can be increased. It is notable that the efficiency of the DBS technique was approved in previous investigations.^16,17

Our results demonstrated that by changing the number of DBS from 0 to 1000, the simulation time (elapsed time), CPU time, and the number of particles in the phase-space file were increased about 7.1, 7, and 36.9 times, respectively, as we expected. The DBS technique splits the collisions into separate collisions and produces photons regarding their value. Therefore, increasing the DBS value results in a higher number of particles. Consequently, these particles must be tracked until their cut-off energies, and this will increase the simulation time. Increasing the DBS value from 0 to 100 increased the simulation uncertainty by about 1.29% on average for the assessed scoring planes. However, increasing the DBS value from 100 to 1000 increases the uncertainty by about 0.71%, and it also increases the simulation time by 315%. Therefore, it is evident that this change in the DBS value is not appropriate for simulation. We did not investigate the DBS values broadly due to the limited time for simulation. For the selected DBS values, it was found that a value equal to 100 is more appropriate for the simulation of the 6 MV photon beam of a linear accelerator.

Previous studies such as Mohammed et al ¹¹ referred to the role of DBS in yielding and reducing of variance in the BEAMnrc code. In another study, it was reported that DBS could increase the photon fluency efficiency for a simulated 6 MV photon beam (with a field size of 10 × 10 cm²) over 8 times higher in comparison with the optimized SBS technique, and over 20 times compared to the UBS technique for the same DBS values. ¹⁷ They also reported that total dose efficiency in a central-axis depth-dose curve improves by 6.4 over SBS at all depths in the phantom. In agreement with the previous investigations, our result showed the efficiency of the DBS technique in reducing the simulation uncertainty.

The results obtained from the current study could be helpful for researchers who want to simulate photon fields. They can use our method for finding the appropriate DBS value, or our initial proposed DBS value (100) for the simulations. There are several limitations in the current study, including the limited number of DBS values due to the limited time for conducting this investigation; therefore, as a suggestion for future studies, these values could be extended to find the appropriate DBS value for simulating photon beams with different field sizes. Furthermore, the dose distribution in the water or patient mimicking phantoms could be calculated to evaluate the effect of changing the DBS value on the precision of the dose calculations.

Conclusion

The effect of the DBS technique on simulation output parameters and the spectrum was investigated for the linear accelerator photon beam (6 MV resulted from a Siemens Primus) to obtain an appropriate DBS value regarding the lower simulation time and the number of particle histories. Generally, using the DBS technique reduces the simulation time or uncertainty. However, increasing the DBS from a particular value (100 in our study) increases simulation uncertainties and times. Therefore, we propose considering a specific DBS value as we obtained in the current research for the Monte Carlo simulation of LINAC photon beams.