If you don't remember your password, you can reset it by entering your email address and clicking the Reset Password button. You will then receive an email that contains a secure link for resetting your password
If the address matches a valid account an email will be sent to __email__ with instructions for resetting your password
Pencil Beam Redefination Algorithm (PBRA) was implemented with three initialization methods and compared to MC.
•
Phase space distribution of the electrons below the applicator was generated using MC simulation.
•
The agreement between PBRA dose calculation and MC is promising for development of an intraoperative TPS for IOERT.
Abstract
Purpose
Similar to other radiation therapy techniques, intraoperative electron radiation therapy (IOERT) can benefit from an online treatment planning system (TPS). Among all the analytical electron dose calculation algorithms, pencil beam redefinition algorithm (PBRA) has shown an acceptable accuracy in inhomogeneities. The input dataset for PBRA includes electron planar fluence, mean direction and root mean square (RMS) spread about the mean direction which had been introduced based on the conventional linear accelerator geometry in former studies. Herein, three methods for implementing PBRA for IOERT system are presented.
Methods
The initialization parameters were identified using Monte Carlo (MC) simulation of a dedicated IOERT system equipped with a cylindrical 10 cm applicator, irradiating a water phantom. Phase space distribution of electrons was recorded on a plane below the applicator. The input dataset was extracted for 2 × 2 mm2 pixels and energy bin width of 1 MeV.
Results
PBRA was implemented with three initialization methods and compared to MC. The 3D gamma analysis of the algorithm with the Formula method, which was in best agreement with MC in a simple water phantom, showed passing rates of more than 99 % for all nominal energies and it was 97.1 % for 8 MeV in the presence of protecting disk and irregular surface. Implementing PBRA on CUDA C++ resulted in 5 s run time for 8 MeV nominal energy in a water phantom.
Conclusions
The agreement between PBRA dose calculation and MC is promising for the development of an intraoperative TPS for IOERT.
Intraoperative radiation therapy (IORT) is a treatment modality in which the tumor bed is irradiated during the tumor removal surgery with a high dose of 9–25 Gy [
]. There are three techniques of delivering IORT, including electron (IOERT), X-ray (Low-kV IORT), and high dose rate brachytherapy (HDR-IORT). Among them, IOERT benefits from more uniform dose distribution, limited penetration depth, and short treatment time which results in its popularity compared to the other two techniques [
]. Today, the significance of IOERT treatment planning is overlooked, and it is only limited to calculation of required monitor units based on one-dimensional (1D) percent depth dose curves. Accordingly, oncologists implement atlases showing the dose distribution in water for different energies and pre-selected applicators [
]; however, dose distribution atlases are associated with doses of a flat surface and a single tissue material. These predicted doses can be significantly different from those of the lumpy surfaces caused by the weight of the applicator or the tumor cut outs [
Costa F, Sarmento S, Sousa O, Phys Med. Assessment of clinically relevant dose distributions in pelvic IOERT using Gafchromic EBT3 films. 2015;31:692-701.
]. In addition, in the case of tumor adjacent to organs at risk, a protecting disk is used with distinct materials compared with the soft tissue (e.g. Teflon and steel) [
]. This setup may cause inadequate dose coverage of the tumor bed and result in cancer recurrence or hot spots in nearby healthy tissues. Therefore, there is a need for a treatment planning system (TPS) to accurately calculate the dose distribution in IOERT.
Some TPSs have been developed for IORT. Except for some in-house developed TPSs, the only certified and clinically available option is Radiance (GMV, Tres Cantos 28,760 Madrid, Spain), in which the IORT procedure is simulated by displaying the virtual position of the applicator superimposed on the patient’s computed tomography (CT) or magnetic resonance image (MRI) [
].Because it is not possible to predict the required placement of the applicator and the tumor and the final geometry of the surface prior to surgery, the treatment planning that is done using pre-surgery images is not accurate. However, Radiance software was also used for online treatment planning and it was shown that using data based on intraoperative actual anatomical conditions can improve the IOERT plans [
García-Vázquez V, Marinetto E, Guerra P, Valdivieso-Casique MF, Calvo FÁ, Alvarado-Vásquez E, et al. Z Med Phys. Assessment of intraoperative 3D imaging alternatives for IOERT dose estimation. 2017;27:218-31.
Widespread use of image-guided surgeries has resulted in the availability of C-arm fluoroscopy, cone beam CT (CBCT), or sonography in most operating rooms. This makes it possible to image the surface geometry post-surgery before performing the dose calculation for a more accurate result. In recent years, suitable methods for online treatment planning have been proposed. The purpose of these methods is to correct and complement the pre-operative treatment planning with the aid of intraoperative images. For this purpose, algorithms with high accuracy in inhomogeneity, e.g. Monte Carlo (MC)-based codes, should be used. However, even fast MC-based algorithms are time-consuming in order of three [
Guerra P, González W, Lallena A, Cal-González J, Herranz E, Udias JM, et al. Monte Carlo based dose estimation in intraoperative radiotherapy. IEEE Nuclear Science Symposuim & Medical Imaging Conference: IEEE; 2010. p. 3069-72.
]. Notwithstanding, the speed of IORT TPS calculation is a crucial factor because it is performed during the operation in the presence of the patient and the surgery team.
Another online TPS is image-guided IORT (igIORT) software for Low kV-IORT [
]. It combines the intraoperative CBCT with pre-op CT and applies a hybrid MC algorithm with Radiance for dose calculation. It results in a planning time of 11–16 min for a vertebra treatment. Furthermore, the accuracy of the dose distribution suffers from metal artifacts in the CBCT images.
Radiance has been implemented using three dose calculation algorithms, including pencil beam (PB) and MC for IOERT and hybrid MC for Low-kV IORT [
]. PB algorithm is widely used to calculate broad beam electron dose distribution in clinical cases based on the Fermi-Eyges theory of multiple Coulomb scattering [
]. However, significant limitations have been observed for non-slab geometries that result in lateral and depth dose discontinuities. Pencil beam redefinition algorithm (PBRA) has been proposed to overcome the limitations of ‘‘semi-infinite’’ slab approximation in the PB model [
]. Considering these two factors is essential in IOERT due to irregular surface geometry caused by the surgery and presence of the applicator and collimation systems and the protecting disk.
Because of previous successful implementations of PBRA in arc therapy [
], its application in IOERT has been investigated. Herein, some methods to supply input parameters for the PBRA algorithm have been studied and compared. The gold standard used for dose distribution comparison was MC simulation.
Materials and methods
Theory
PBRA algorithm improves the accuracy of PB by redefining the parameters of each pencil beam at each depth. These parameters include electron planar fluence, mean direction and root mean square (RMS) spread about the mean direction of each pencil beam, which are defined in every pixel and energy bin at each depth [
]. A dedicated IOERT accelerator, LIAC (Light Intraoperative Accelerator–SordinaSpA, Italy) is studied to obtain the input parameters of the PBRA algorithm.
], the initial beam is assumed to be mono-energetic, and the fluence is independent of the position. The mean angle of each pencil beam was only dependent on the divergence from the virtual accelerator source to the studied pixel. The RMS angular spread was obtained from an in-air measurement of penumbra. Both latter parameters were considered to be independent of the electron beam energy.
] modified the developed PBRA by considering the initial beam with a poly-energetic spectrum. They assumed a uniform fluence emerging from a single virtual source and a discrete poly-energetic spectrum for the initial electron beam. In addition, the RMS spread of the beam radial distribution (σθ) was assumed independent of position based on the methodology proposed by Hogstrom et al. [
]. However, the collimation used in dedicated IOERT accelerators is mainly long cylindrical applicators without a clear SCD definition. Therefore, a new strategy for its calculation is needed. In addition, formerly implemented PB and PBRA algorithms did not take initial scattered electrons into account, while this contribution in IOERT is even higher than in conventional electron linear accelerators (linacs) with different energy and angular distributions [
]. To consider the true fluence and direction distributions in IOERT-specific accelerators, MC simulation was used. Then, the PBRA input parameters were determined as a function of position and energy.
PBRA initial parameters
The input dataset
MC simulation of LIAC 12 MeV IOERT system has been performed previously for electron beams with various energy (nominal energy of 6, 8, 10, and 12 MeV). The beam collimation systems and the simulation results have been validated by their experimental counterparts in the water phantom [
]. In this study, the simulations with the 10 cm diameter applicator, which was mentioned above, are performed. Phase space distribution (PSD) of the electron particles was recorded using the PTRAC card of MCNPX code [
]. The PSD surface was positioned under the changeable applicator and recorded the position, direction, and energy of 5,000,000 crossing source particles.
To supply the initial parameters for PBRA and to redefine them in every plane in the Z direction, electrons were grouped in 2 × 2 mm2 pixels. The energy was considered discrete with 1 MeV bins. The planar fluence was equal to the number of electrons in each pixel and energy bin. Direction of electrons should be calculated in both X_Z and Y_Z planes, separately. According to Shiu et al. [
], in the X_Z plane, the mean and RMS spread of the particle directions should be calculated based on their angle with the Z axis (φ angle in Fig. 1). Determination of this angle from the output data from the PTRAC card was done based on the following equations. The initial directions in the Y_Z plane are the transpose of the X_Z plane’s directions, according to symmetry.
(1)
where u, v, and w are cosine of the crossing particle directions, given by the PTRAC card, with X,Y, and Z axes, respectively.
Fig. 1The coordinate system used for defining the electron pencil beam position and direction.
]. Therefore, in order to obtain the mean and RMS spread about the mean directions in each pixel and energy bin, the histogram of the angles of the recorded electrons was fitted with a Gaussian function. The mean and standard deviation (SD) of the fits corresponded to the mean direction and RMS spread about the mean direction, respectively. To obtain better goodness of fit measures, 20 angle bins per each energy bin were considered, an appropriate initial guess for fitting parameters was applied to the curve fitting procedure, the recorded particles of the two semicircles about the Y axis were combined because of system’s symmetry, and the fits with less than 10 data points were ignored.
MATLAB (MathWorks, R2017a) was used for the analysis of the PSD data of the MCNP code and fitting the angular distribution curves. In addition, the mean and SD angle changes as a function of electron position and energy were obtained. Some equations were proposed for the position and energy-dependent direction parameters based on fitting the data points.
Finally, PBRA implementations with different sets of initial parameters were compared. The three sets of initial parameters include:
1)
Parameters extracted from the PSD data using the PTRAC card (PTRAC).
2)
Fitting formulas for the mean and SD angles presented in section 3.2 (Formula).
For the third set, the mean angle for each pencil beam in the X – Z plane was calculated based on the beam divergence equal to arctan (x/SSD), where SSD stands for the geometric source–surface distance in LIAC that is usually equal to 713 mm [
] on RMS of projected scattering angle. Because of the ambiguity in the determination of SCD in the equations for IOERT, as described in section 2.1, SSD was replaced. It should be noted that in all three PBRA implementations, the planar fluence was imported from the PSD surface data of PTRAC card of MC simulation with discretized energy.
For comparison of PBRA performance using the three different input datasets, a water phantom of 20 × 20 × 15 cm3 dimension with a mesh dimension of 2 × 2 × 5 mm3 was simulated in both MCNP and PBRA, and the recorded energy depositions were compared in terms of gamma index.
Implementation of PBRA algorithm
Determination of correction factors
Energy-dependent correction factor (C(E)) plays a crucial role in the accuracy of PBRA dose calculation. It accounts for the neglected physics of electron interactions (e.g., small-angle scattering approximation) [
] was used for C(E) determination, i.e. fitting a high-order polynomial to the energy dependent correction factor values. The gold standard in this study was the dose distribution in the water phantom calculated by the validated MC. The electrons were transported while the applicator was placed perpendicular to the phantom surface and the PDD curve was calculated using the energy deposition tally. The absorbed doses were recorded in 2 × 2 × 5 mm3 cells. The number of histories was enough to achieve less than 3 % uncertainty. The degree of polynomial function fitted to the correction factors versus energy was chosen in such a way to obtain the least possible dose difference between the PBRA and MC simulation.
Consideration of the backscatter of the protecting disk
], so the back scattered electrons are not taken into account. To consider the effect of the backscattered electrons of the protecting disk on the dose distribution, MCNP was used. The disk was simulated in a water phantom at a 90 % depth dose [
Baghani HR. Possibility evaluation of intraoperative imaging, treatment planning and simultaneous dosimetry during breast cancer intraoperative radiotherapy and comparing the results of preoperative and postoperative treatment planning in this treatment technique. Iran: Shahid Beheshti University; 2014.
], and energy deposition mesh tally was used to calculate the dose with and without the protecting disk at different nominal energies. Backscattering factor (BSF) was calculated in percentage in each voxel as follows:
where is the dose in the presence of the protecting disk and is the dose without the disk.
In the presence of the protecting disk, corrected dose of PBRA can be achieved by multiplying BSF by the PBRA calculated dose of the voxels on the top of the disk.
PBRA parallelization
In order to implement the algorithm faster, PBRA was written in C++, and the Formula-based initialization method was used based on the obtained results. It was parallelized using CUDA toolkit 8.0 in the Nvidia Tesla K20X graphic card.
Implementing the PBRA requires 6 nested loops. The outermost loop indicated the depth and was not parallelized, because the beam characteristics were redefined at each depth. The innermost loop iterated on the energy bins with the highest value of 13; therefore, it was not efficient to be parallelized in GPU. Another 4 parallelized loops were responsible for accounting the X and Y coordinates of the pixels in the given depth and the pixels of the previous depth. The calculations inside these loops were independent so the loops could work in parallel. The number of threads that can run simultaneously is the number of pixels in the current depth multiplied by the number of pixels in the previous depth.
In the case of 10 × 10 cm2 field divided by 2 × 2 mm2 pixels, the number of threads was 514. A unique combination of 4 numbers from 0 to 50 should be extracted from each thread number to specify the X and Y coordinates of the current and previous depth pixels. In other words, the number should be converted from base 10 to base 51, and each digit represents the required number. Fig. 2 shows the number of thread which is responsible for calculating the deposition energy of the electron which is traveling from pixel (1,49) in the previous depth to pixel (50,0) in the current depth.
Fig. 2An example of the relation between the thread number and the coordinates of the pixels.
In addition to using MCNP to supply PBRA initial parameters, MCNP code was used to validate the dose distribution calculation performed by PBRA (implemented by MATLAB). It was accomplished by calculating the 2D percent depth dose (PDD) and transverse dose profile (TDP) and 3D absolute dose distribution. A 20 × 20 × 15 cm3 cubic water phantom was used with the applicator perpendicular to the phantom surface. The energy deposition was recorded in voxels with the size of PBRA voxels (i.e., 2 × 2 × 5 mm3) with the largest dimension in the depth direction (Z axis). Comparison between the PBRA results and MC simulation was performed in terms of gamma index with criteria of 3 % dose difference (DD) and 3 mm distance to agreement (DTA) [
To show the applicability of PBRA in a clinical case, an IOERT treatment of breast was simulated in MCNP using the mathematical phantom Revised ORNL (Adult) [
Permission to use and modify the relevant file was obtained from the author of the article.
The phantom was rotated horizontally and the applicator was placed on the left breast. The breast size was changed to fit the average breast size of Iranian women [
]. The part of breast under the applicator was considered flat with two incisions in the middle and the corner of the field of view (FOV). The depth of the incisions was 5 mm and the protecting disc was placed at the depth of 1.5 cm (Fig. 9).
Baghani HR. Possibility evaluation of intraoperative imaging, treatment planning and simultaneous dosimetry during breast cancer intraoperative radiotherapy and comparing the results of preoperative and postoperative treatment planning in this treatment technique. Iran: Shahid Beheshti University; 2014.
], for the patients with 1 to 2 cm tumor bed, electrons with energy of 8 MeV are used. Due to the complex geometry, this simulation was very time consuming. Therefore, the energy cut off card of 0.55 MeV for electrons and 0.010 MeV for photons [
Righi S, Karaj E, Felici G, Di Martino F, J Appl Clin Med Phys. Dosimetric characteristics of electron beams produced by two mobile accelerators, Novac7 and Liac, for intraoperative radiation therapy through Monte Carlo simulation. 2013;14:3678.
] was applied and energy deposition mesh tally with a mesh dimension of 2 × 2 × 5 mm3 was used to calculate the energy deposition. The same voxel size, energy and geometry were used in PBRA with formula-based initialization. To consider the back scatter of the protecting disk in the PBRA, the method described in section 2.3.2 was used. Comparison between the dose distributions of PBRA and MCNP was made in terms of gamma index with criteria of 3 % DD and 3 mm DTA.
Results
Determination of direction parameters
Fig. 3. shows two typical Gaussian fits to angle distributions of initial electron beams recorded on the PSD plane in 9–10 MeV and 0–1 MeV energy bins within the pixel on the central axis. They were recorded from the 10 MeV nominal energy electron source and fitted using the curve fitting tool of MATLAB. It can be seen that although there is about 25 orders of magnitude difference in the fluence of the two energy bins, both histograms are Gaussian with a peak near zero degree. Due to the greater probability of large-angle scattering at low energies [
Svensson H, Almond P, Brahme A, Dutreix A, Leetz HK. Radiation dosimetry: electron beams with energies between 1 and 50 MeV. Report 35. International Commission on Radiation Units and Measurements (ICRU)2016. p. NP-NP.
], the peak of 0–1 MeV energy bin shows more significant deviation from zero. The curve fittings were done for every pixel and energy bin on the PSD plane. The correlation coefficient (R2) of the fits for energies higher than 3 MeV was higher than 0.9 and, for all energies, it was higher than 0.7. The lower R2 value for lower energies was not crucial due to their negligible electron fluence.
Fig. 3Angle distribution of (a) 9–10 MeV and (b) 0–1 MeV energy bins in the pixel on central axis for 10 MeV LIAC operation mode.
Fig. 4 shows the dependence of the mean angle () on energy in different pixels on the X axis. For better understanding, at each nominal energy, the curves of 5 pixels are highlighted. The red and yellow curves are associated with pixels near the edge of the applicator on opposite sides (46 mm from the center), the green curve is related to the pixel on the center of FOV, while the blue and purple curves are related to the pixels in between (22 mm from the center). The dependence of mean angle () on energy in different pixels shows a regular divergent pattern from a point shown by an arrow on each nominal energy plot. The irregular pattern before that point can be neglected because of providing less than 10 % of the total fluence on the PSD plane. In order to achieve a relationship between and energy and position, the curve fitting was done with the data points weighted by the third power of their fluence. Equation (2) shows the applied fitting formula based on the fact that the linear trend of the curves obeys the beam divergence based on the triangular relations.
(2)
where a and b are the fitting parameters and × is the X coordinate of the pixels relative to the central axis. The rightmost term in equation (2) shows the convergence of the curves. It was one of the fitting parameters, which resulted in a value close to E0/2 (E0 is the nominal energy of the accelerator). By replacing the parameter with E0/2, the 95 % confidence intervals (CI) did not change considerably and even improved in some cases. The parameter values as well as their 95 % confidence intervals (CI) for different E0 are presented in Table 1.
Fig. 4Mean initial electron angle versus energy for 6, 8, 10, and 12 MeV nominal energies. Different colors show the position of the 2 × 2 mm2 pixels on the phantom surface along the X-axis. The bold lines and the arrow is described in the text.
Eq. (2) can be approximated with a linear function of x. This provides the advantage of more straightforward implementation and less fitting parameters.
(3)
The variables can be defined as before and the fitting parameter, A, as well as goodness of fit values are presented in Table 1.
Based on Fig. 5, energy dependence of SD of the angles is exponential and approximately independent of position. Therefore, the mean and standard deviation of the SD were calculated as a function of X positions. The average was used for curve fitting, while the inverse of standard deviation was used as the weighting factor.
(4)
where a and b are fitting parameters. The fitting parameters and their 95 % CI for different E0 values are given in Table 2.
Fig. 5Standard deviation of initial electron angles versus energy for 6, 8, 10, and 12 MeV nominal energies. Different colors show the position of the 2 × 2 mm2 pixels on the phantom surface along the X-axis.
Comparison between the three initialization methods of PBRA was performed based on the procedure described in Section 2.2.2. The energy-dependent correction factors for the three methods were calculated using the same polynomial degree but based on the corresponding dataset, for each nominal energy. Gamma passing rates for different initialization methods are presented in Table 3 in the whole electron beam FOV and up to the last nonzero depth dose in PBRA. In Table 1, PTRAC, Formula, and Hogstrom indicate the three initialization methods as described in section 2.2.2. It can be seen that all the passing rates are more than 96 % except for Hogstrom initialization method for 6 MeV. Moreover, Formula-based initialization results in the most accurate dose distribution among the others. Thus, it was chosen as the initialization method for PBRA-based dose calculations performed in the subsequent sections. Most of the gamma indices more than one occurred in the field edge.
Table 3Gamma passing rates for different initialization methods.
Fig. 6, Fig. 7 show the PDD and TDP curves obtained from MC simulation and PBRA-derived dose distribution using formula-based initial parameters, respectively. The deposited dose in each voxel was assigned to the middle point of each cell in depth direction. Both PDD and TDP curves were normalized to their maximum values on the central axis. In all MCNP simulations, the maximum uncertainty was 2.8 %. In all PDD curves except in the 10 MeV curve, obtained gamma indices were less than 1 up to the last non-zero dose point, and up to two last non-zero points for the 10 MeV curve. All TDP curves passed the gamma index criteria within the FOV of 10 cm.
Fig. 6Comparison of PDD curves obtained by PBRA and MC in the water phantom for 6, 8, 10, and 12 MeV nominal energies.
In addition to one dimensional analysis, the 3D dose distribution was calculated using both MC and PBRA methods. Gamma analysis of the 3D dose distribution was performed with a mesh dimension of 2 × 2 × 5 mm3. The obtained gamma values are reported in the whole electron beam FOV and up to the last nonzero depth dose in PBRA (Fig. 8). The voxels outside the FOV were not taken into account where MC error was not acceptable. The 3D gamma passing rates for 6, 8, 10, and 12 MeV were 99.3 %, 99.9 %, 99.7 %, and 99.9 %, respectively.
Fig. 83D absolute dose distribution and gamma indices for 6, 8, 10, and 12 MeV nominal energies.
{Han, 2006 #28}, in (a) Y-Z plane at X = 0, (b) X-Z plane at Y = 0, and (c) X-Y plane at Z = 3 mm from the applicator surface. PBRA and MCNP isodoses of 95 %, 100 %, and 105 % in the Y-Z plane are shown. Dot dash lines indicate PBRA dose distribution and solid lines are related to MCNP. The red color in the figures represents the air in the applicator and incisions, grey shows the breast tissue, green and orange colors show the protecting disk layers from Teflon and steel, cherry red shows the chest and vertebrae, blue demonstrates the lung, and pink shows the heart.
An IOERT treatment of breast was simulated in Revised ORNL Phantom for 8 MeV nominal energy. The 95 %, 100 %, and 105 % isodoses for PBRA and MCNP in the YZ plane are shown in Fig. 9. The dose is normalized to the dose of a point in 1 cm from the applicator surface and 3 cm from the edge of the field. The shift relative to the central axis was due to presence of a scar on the central axis. Under the protecting disk, the dose of PBRA was zero while it was not the case for MCNP. However, the statistical error of the simulation was large and unacceptable. Also, there was an unacceptably large error in MCNP-derived dose distribution in the incisions made under the applicator, too. Excluding the aforementioned regions, the maximum error was 1.4 %. The gamma passing rate in the acceptable voxels between the applicator surface and the protecting disk was 97.1 %.
Discussion
Mean and SD of deflection angles
Herein, the applicability of PBRA for dose calculation in IOERT has been examined. In order to provide the input dataset for PBRA calculation, MC simulation was used. Fig. 4 shows that the mean angle of the electron’s trajectory is proportional to energy with a position dependent divergence after a point about half of the nominal energy (shown by an arrow). The reason for this proportionality is that a larger deflection angle of electrons is likely to be found with the increase in distance of the pencil beam pixel from the central axis. As can be seen in Fig. 4 for different nominal energies, the average angle in the central axis pixel is about zero (green highlighted curves). For the region after the arrow marker, the highest mean angle was found for the farthest pixel (half FOV diameter, 5 cm from the central axis). Based on the typical IOERT SSD of 71.3 cm, the beam divergence in this pixel is equal to degrees, which is in accordance with the data presented in Fig. 4.
In contrast with mere position-dependent approach used by Shiu et al. [
], it was found that the mean angle is a function of electron energy.
According to Fig. 5, the RMS spread about the mean direction of electrons demonstrates an exponential relationship versus energy while an inverse dependence was assumed by Boyd et al. for conventional linac [
Svensson H, Almond P, Brahme A, Dutreix A, Leetz HK. Radiation dosimetry: electron beams with energies between 1 and 50 MeV. Report 35. International Commission on Radiation Units and Measurements (ICRU)2016. p. NP-NP.
]. Therefore, the standard deviation of the Gaussian fits to the deflection angles increases.
Different initialization methods
Gamma indices indicate that PBRA initialization using the Formula method is more consistent with the MC calculation. The use of Formula initialization method outperforms the counterparts in the high dose gradient regions rather than the middle low gradient regions. Therefore, in general, the use of fitting formulas to Gaussian direction distribution can be considered a good choice for direction parameter initialization in PBRA, even better than storing PSD data of an IOERT accelerator. Using energy and position – dependent mean and RMS spread about the mean formulas increase the simplicity and decrease the calculation time as well as required memory. It is worth noting that the coefficients of the presented equations will vary for different accelerators and applicator shapes and sizes, but the overall relationship and their rationale will not be changed.
PBRA validation
Gamma passing rates presented in Fig. 8 and Table 3 show a good agreement between the dose distributions of PBRA and MC simulation. Fig. 7 illustrates a high discrepancy in the penumbral region of the dose profiles. This is a result of modeling the scattering process in the form of a Gaussian distribution in the Fermi-Eyges model [
]. However, the scattering is in fact completely unstructured and irregular. Due to more scattering in the corner of the FOV, the error introduced by this approximation is maximized. Another source of discrepancy is neglecting the bremsstrahlung radiation and delta electrons and assuming local deposition of electron energy which lead to over-estimation of the calculated dose. PBRA deficiency in dose estimation in the penumbral regions can be seen in the work presented by Boyd et al. [
] as well. A distance as high as 1.5 cm between the calculated and measured isodose lines at the maximum depth dose of 2 cm for 9 MeV energy can be observed in their presented figures. The performance of the PBRA in inhomogeneity which was tested in a typical incision was acceptable excluding the air within the incision and the tissue under the protecting disk.
PBRA was parallelized in CUDA C++. It took 5 s to calculate the dose distribution based on PBRA in a 10 × 10 × 8 cm3 water phantom with a grid size of 2 × 2 × 5 mm3. The time consumed in the presence of inhomogeneities will not be changed considerably because neither the number of threads working simultaneously, nor the number of variables which should be calculated in a thread will be changed. Calculation time depends on the field size and the depth at which the electron fluence becomes zero (e.g. after the depth of shielding disk or the energy dependent electron range).
Conclusion
In this paper, PBRA was implemented for IOERT treatment planning. The algorithm passed the gamma criteria except near the field edges. The initial parameters required to feed into this algorithm were extracted using three methods and they were compared. It was found that PBRA is capable of dose calculation with an acceptable accuracy and easy implementation using the proposed fitting formulas for its initialization.
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
References
Gunderson L.L.
Willett C.G.
Calvo F.A.
Harrison L.B.
Intraoperative irradiation: techniques and results.
Costa F, Sarmento S, Sousa O, Phys Med. Assessment of clinically relevant dose distributions in pelvic IOERT using Gafchromic EBT3 films. 2015;31:692-701.
García-Vázquez V, Marinetto E, Guerra P, Valdivieso-Casique MF, Calvo FÁ, Alvarado-Vásquez E, et al. Z Med Phys. Assessment of intraoperative 3D imaging alternatives for IOERT dose estimation. 2017;27:218-31.
Guerra P, González W, Lallena A, Cal-González J, Herranz E, Udias JM, et al. Monte Carlo based dose estimation in intraoperative radiotherapy. IEEE Nuclear Science Symposuim & Medical Imaging Conference: IEEE; 2010. p. 3069-72.
Baghani HR. Possibility evaluation of intraoperative imaging, treatment planning and simultaneous dosimetry during breast cancer intraoperative radiotherapy and comparing the results of preoperative and postoperative treatment planning in this treatment technique. Iran: Shahid Beheshti University; 2014.
Righi S, Karaj E, Felici G, Di Martino F, J Appl Clin Med Phys. Dosimetric characteristics of electron beams produced by two mobile accelerators, Novac7 and Liac, for intraoperative radiation therapy through Monte Carlo simulation. 2013;14:3678.
Svensson H, Almond P, Brahme A, Dutreix A, Leetz HK. Radiation dosimetry: electron beams with energies between 1 and 50 MeV. Report 35. International Commission on Radiation Units and Measurements (ICRU)2016. p. NP-NP.
02076-2&id=gr1.jpg){kind=link}
02076-2&id=gr2.jpg){kind=link}
02076-2&id=gr3.jpg){kind=link}
02076-2&id=gr4.jpg){kind=link}
02076-2&id=gr5.jpg){kind=link}
02076-2&id=gr6.jpg){kind=link}
02076-2&id=gr7.jpg){kind=link}
02076-2&id=gr8.jpg){kind=link}
02076-2&id=gr9.jpg){kind=link}