Introduction

In many modern radiation therapy treatment planning systems, model-based algorithms, such assuperposition convolution (S/C) algorithms, are reliable computational dose distribution methods with a reasonable speed [ 1 - 3 ]. In the S/C algorithm, the dose distribution is a result of convolving the ‘total energy released per unit mass’ (TERMA) by ‘energy deposition kernel (EDK)’ in all points of the dose calculation volume [ 2 , 4 ].

$D=T\left(\mathrm{E,r}\right)⨂K\left(\mathrm{r-r}\right)$ (1)

Where T(E,r ́) is TERMA, K(r-r ́) is EDK, E is photon beam energy, r is interaction site and r ́is dose receiving point. EDK is an invaluable tool for qualitative aspects of dose distributions and shows the energy distribution pattern around the single interaction site of the primary photon. It can be determined using empirical and analytical methods, but the Monte Carlo is the only practical approach. EDK plays a key role in the precision of the model-based dose calculation method, thus accuracy and any applied assumptions or approximations in EDK calculation will seriously affect dose calculation end results accuracy [ 5 - 7 ].

In 1984, Ahnesjo derived EDK theoretically by using deconvolution in measured narrow beam dose distributions [ 8 ]. TR Mackie et al. alternatively calculated monoenergetic EDK using Monte Carlo simulation [ 9 ]. Subsequently, TR Mackie et al. produced EDK using electron gamma shower (EGS) in water with multiple scattering for relatively high energy ranges as a database for convolution dose calculation in radiation therapy [ 6 ].

In 2000, after developing and extending the original EGS code, NRC’s EGS general-purpose software toolkit was released [ 10 ]. Three years later, Mainegra-Hing E et al. adapted and introduced EDKnrc as one of the NRC’s derived spherical geometry user code which can calculate energy deposition kernels for photons or electrons forced to interact at the center of a spherical phantom in the entire specific voxelized phantom [ 11 ].

EDKnrc is used for EDK calculation in the homogenous sphere in energy ranges between a few keV up to 55 MeV and considers different contributions from each photon scattering order (eg. primary, first scatter, second scatter, multiple scattering, and radiative transfer) [ 12 ]. However, in the patient body containing various inhomogeneous media, adjacent voxels contain various types of inhomogeneities. Therefore, generating EDK using EDKnrc for all voxels made of different medium is not practical.

To adapt water EDK to entire non-water media, O’Connor’s theorem and Mohan et al. density scaling method are usually used to scale EDK by the relative electron density between the interaction and dose receiving points, voxel-by-voxel [ 13 - 15 ]. Nevertheless, some studies in the literature have highlighted the low accuracy of this method near the low-density media interfaces such as the lung in S/C dose calculation [ 15 , 16 ]. Regarding this fact, M. K. Woo et al. investigated the validity of the density scaling method for 2 and 6 MV photon beams in primary electron transport for air and cork and their effects on dose calculation accuracy. They compared their results with kernels obtained by EGS code. The final results showed a large discrepancy between these two methods, especially for air and 6 MV beam, but no solutions were provided for this problem [ 15 ]. Jessie Y. Huang performed a comprehensive research on various refinements for EDKs such as a partial study on the density scaling method [ 2 ].

As mentioned, the proper scaling of polyenergetic kernels, especially in heterogeneous environments, is one of the most important links in the chain of model-based algorithms. However, to the best of our knowledge, published studies on kernel scaling accuracy have only been touched upon and except for a few general theoretical relationships; there has been no well-documented analysis of the mentioned scaling formula. Thus, this study aims to investigate EDK scaling formula accuracy in the presence of lung and bone inhomogeneity in comparison with EDKnrc-generated data to ensure its generalizability to the heterogeneous environment. In addition, based on the nature of the kernels, a number of correction factors were proposed and implemented on the theoretical scaling formula to create more consistency between these data.

Material and Methods

This part of the current theoretical-practical study follows some steps, including generation of EDKs using EDKnrc user code, generation of inhomogeneous EDKs using density scaling method, density scaling method correction.

Generation of EDKs using EDKnrc user code

Water EDK

EDKnrc user code of the EGSnrc package was used to generate EDKs according to ‘Mackie’s method [ 7 ]. In this code, photons are forced to collide at the center of a 60/ρ cm radius theoretical spherical phantom, in which ρ describes the mass density of the spherical phantom. This phantom is voxelized into 1152 voxels with 48 cones, separated by a polar angle of 3.75° and 24 radial shells. Consequently, the 3.75-degree angle cone is solid, and each of its voxels is in the form of a truncated solid cone while in other angles, these solid shapes turn into the shell geometry (Figure 1).

Figure 1. The voxels of energy deposition kernel (EDK) phantom: a) One voxel in a certain angle and radius (by permission of Mackie, T.R, 1988), b) Six voxels in 3.75 degree for 0 to 4 cm radius.

At first step, the 6 MV water EDK (as total and primary stored energy data in each voxel) was calculated, in which data value was maximum in the central voxels, and gradually decreased by increasing the radius. To find the maximum radius with a non-negligible amount of stored energy, its angular and radial distributions were plotted and analyzed.

Inhomogeneous EDKs

Different areas of the human body are composed of different heterogeneities. Thus, to simulate these various heterogeneities in the path of the photon beam, six different special inhomogeneous 6 MV EDKs were produced similar to water EDK. Four of six EDKs consisted of some spherical shells of bone and lung in the certain radii rather than water. For the rest, the interaction point medium (sphere of 0.05 cm radius) was changed into lung or bone in addition to the features of the previous EDK (Table 1). For all simulation 300 × 10⁶ histories were run to achieve acceptable uncertainties in the small voxels near the point of interaction. Photon and electron cutoff energies (PCUT and ECUT) were set to 0.01 and 0.521 MeV, respectively [ 11 ].

Considering all the above details, the total and primary EDKs data representing deposited energy per particle in each scoring voxel were calculated. To make reasonable comparisons between calculated and EDKnrc-generated EDKs, all the 1152 voxels data for each EDK were divided by each voxel volume and imported into EXCEL software (Version 2018). The volume for each voxel of phantom was calculated analytically for a series of truncated spherical cones, according to the following equations:

$\mathrm{volume}\left(\mathrm{i,j}\right)=\frac{2}{3}\pi \left[\right({\mathrm{cos\theta }}_{\mathrm{i-1}}-{\mathrm{cos\theta }}_i\left)\right(r_j^3-r_{\mathrm{j-1}}^3\left)\right]$ (2)

θ_i=0°,3.75°,11.25°,…,180°

r_j=0,0.002,0.01,…,60cm

In which θ_i and θ_i-1 are semi-vertex inner and outer angles and r_j and r_j+1 are small and big radii of truncated conical shell voxels, respectively and volume (i,j) refer to these voxel volumes. Normalization was then performed using dedicated written MATLAB code, version 2013.

Generation of inhomogeneous EDKs using density scaling method

Six series of water EDK data were considered that each of them had heterogeneities in special radii similar to EDK-generated ones in Table 1. To convert these water EDKs data to inhomogeneous media kernels, they were scaled according to density scaling method. The density scaling method devised by Mohan et al. is based on the assumption of O’Conner’s theorem, i.e. the energy propagation is linear between interaction (S point) and interesting points (r point), related equations are as follows [ 17 , 18 ]:

$h_{\mathrm{E,inhom}}\left(\overrightarrow{s},\overrightarrow{r}\right)={\rho }_{\mathrm{rel}}\left(\overrightarrow{s}\right){\rho }_{\mathrm{rel}}^2\left(\overrightarrow{s},\overrightarrow{r}\right)h_{\mathrm{E,}{H}_{2}O}\left({\rho }_{\mathrm{rel}}\right(\overrightarrow{s},\overrightarrow{r}).(\overrightarrow{r}-\overrightarrow{s}\left)\right)$ (3)

${\rho }_{\mathrm{rel}}\left(\overrightarrow{s},\overrightarrow{r}\right)={|\overrightarrow{r}-\overrightarrow{s}|}^{-1}{\int }_{\overrightarrow{s}}^{\overrightarrow{r}}{\rho }_{\mathrm{rel}}\left(\overrightarrow{s}\right)ds,{\rho }_{\mathrm{rel}}=\frac{\rho \left(\overrightarrow{s}\right)}{\rho H_2}$ (4)

Where $h_{\mathrm{E,inhom}}\left(\overrightarrow{s},\overrightarrow{r}\right)$ is inhomogeneous media kernel, $h_{\mathrm{E,}{H}_{2}O}\left({\rho }_{\mathrm{rel}}\right(\overrightarrow{s},\overrightarrow{r}).(\overrightarrow{r}-\overrightarrow{s}\left)\right)$ is water EDK, and ${\rho }_{\mathrm{rel}}\left(\overrightarrow{s}\right)$ is relative electron density of interaction site. ${\rho }_{\mathrm{rel}}\left(\overrightarrow{s},\overrightarrow{r}\right)$ is provided by equation 4. Setting the EDKnrc user code-generated data as a reference, the relative difference between these values and density scaled EDK was calculated and plotted versus radii.

Density scaling method correction

Equation 3 considers the second power of average relative electron density of the path $\left({\rho }_{\mathrm{rel}}^2\right(\overrightarrow{s},\overrightarrow{r}\left)\right)$ and the relative electron density of the interaction site ${\rho }_{\mathrm{rel}}\left(\overrightarrow{s}\right)$, without considering the attenuation effect of the heterogeneity. To take this effect into account, after passing the photon beam through each area of heterogeneity, the exponential attenuation coefficient function of the inhomogeneous voxel was applied on the equation (3). In addition, for each voxel containing heterogeneity in the photon beam path, the ${\rho }_{\mathrm{rel}}^{\alpha }$ coefficient was considered. The α value was tuned to minimize the relative error of the scaling method-calculated EDK. According to factors in the scaling method, inhomogeneous EDKs and their relative errors were recalculated.

Results

Water EDK

The analysis of the energy deposition for water EDK in the radial part for a total 360 and 3.75 degrees is demonstrated in the semi-logarithmic scale (Figure 2). It can be observed that the energy is deposited mainly into a radius of less than 0.3 cm. Therefore, in this study, the degree of compatibility of two different methods for EDK generation was investigated up to this radius.

Figure 2. Radial energy distribution of water energy deposition kernel (EDK) in the semi-logarithmic scale for a total 360 and 3.75 degrees.

Inhomogeneous EDKs

According to Figure 3, different arrangements of the heterogeneities in the beam path (Table 1) cause different unpredictable changes in the energy distribution pattern around the interaction sites. These lead to production of different inhomogeneous media kernels data depending on the location of the interaction and the heterogeneity around it. These kernel data in different heterogeneity conditions (k1 to k6) have various deviations with respect to the water EDK, leading to fluctuation of plots in Figure 3.

Figure 3. Radial energy deposition kernel (EDKs) in semi-logarithmic scale for 3.75 degree for different inhomogeneous arrangements of Table 1 (k1 to k6).

As a result, it is unfeasible to consider all possibilities of inhomogeneous EDK within media with random-heterogeneous patterns. To address this problem, the scaling method was used to convert the water kernel into each heterogeneous kernel based on the distribution of heterogeneities around each interaction point of interest.

Generation of inhomogeneous EDKs using the density scaling method

Based on Figure 4, in general, there is a relatively good agreement between calculated density scaling data and EDKnrc-generated ones before and after the heterogeneous voxels, except for within them, in which a spike has been formed in the curves. By applying the correction factors, the relative differences between EDKnrc-generated as a reference and corrected /uncorrected density scaled kernels were calculated and plotted versus radii at Figure 5.

Figure 4. The relative errors between energy deposition kernel (EDKnrc)-generated as reference data with corrected and uncorrected density scaled kernels for different inhomogeneous arrangement of Table 1 (k1 to k6).

Figure 5. Energy distribution pattern (radius energy deposition kernel (EDK)) generated by EDKnrc user code, density scaling and corrected density scaling for different inhomogeneous arrangements of Table 1 (k1 to k6).

The optimum α values for ${\rho }_{\mathrm{rel}}^{\alpha }$ were determined 0.8 and 1.2 for lung and bone, respectively. It was observed that the performed correction had effectively reduced the errors of the EDK No. 1 to 5, while for EDK No. 6 within lung as an interaction site, these correction factors do not work properly.

Discussion

In this study, we evaluated the consistency of the inhomogeneous EDKs produced by EDKnrc user code and those calculated by the density scaling method. Our result has indicated that although the errors in the theoretical method for calculating inhomogeneous EDKs were relatively acceptable, they could be significantly reduced due to the attenuation coefficient and ${\rho }_{\mathrm{rel}}^{\alpha }$ parameter for the heterogeneous voxels.

EDK illustrates the energy pattern around the single interaction site of the primary photon in the spherical phantom (Figure 1). According to Figure 2, presenting the EDK data in radial parts for two different angles has shown that a significant portion of the energy is distributed up to a radius of 0.3 cm, thus we restricted our studies up to this radius. The studies were also performed at an angle of 3.75 degree for the central cone.

Based on the semi-logarithmic plot in Figure 3, the presence of the different heterogeneities in the water kernel phantom, causes various energy distribution patterns around the interaction site. Furthermore, the difference between the kernels in place of the heterogeneous voxels is more obvious. In (S/C) dose calculations algorithm, the heterogeneous EDKs at any point of the irradiated volume must be specified. Although these EDKs can be calculated using EDKnrc user code for any given set of heterogeneities, it is impractical at a reasonable speed. The solution for this challenge is the use of the density scaling method, in which the heterogeneous kernel can be theoretically calculated at any point in the volume of interest.

To check the accuracy level of the density scaling method, the relative difference between the code-generated and the theoretical EDKs was calculated. For dose calculations purpose, due to the spherical shape of the voxels in the EDK’s phantom, the EDK data must be distributed in Cartesian or spherical coordinates. Therefore, the calculated errors will be distributed at an angle of 360 degree.

Thus, except for kernels No. 4 and 6, in which the lung heterogeneities were close to the interaction site, the error was acceptable. Subsequently, by applying the attenuation coefficients and ${\rho }_{\mathrm{rel}}^{\alpha }$ factors, the relative difference between the theoretical and code-generated EDK was recalculated (dotted line in Figure 4). As seen, the correction results in the significant error reduction, especially at the heterogeneous voxels.

EDKnrc-generated, density scaling and corrected-density scaling method EDKs displayed in a semi-logarithmic scale in Figure 5, demonstrate that the corrected method data are in excellent agreement with data obtained from the EDKnrc user code. This adaptation has been partially reduced in EDK No. 4, in which lung heterogeneities are located near the interaction site. However, the accuracy of the calculations is still valuable. In EDK No. 6, the lungs inhomogeneity is an interaction medium, thus it seems that our correction factors have not worked well as the other EDK inhomogeneities arrangements. Therefore, further research is needed in these cases.

Conclusion

EDK is the most basic component of the (S/C) dose calculation algorithm, especially in real heterogeneous environments. Considering the different points of the volume under irradiation, the spherical EDK phantom will have diverse arrangements of inhomogeneities, and it is not possible to be generated using EDKnrc code within a reasonable time. The fastest practical way to calculate inhomogeneous EDK is using the density scaling method.

The present study indicates that although the density scaling method calculates inhomogeneous EDK with acceptable accuracy, the error values are significant at the location of non-water voxels. Taking into account the aforementioned correction factors, heterogeneous EDKs are calculated with an accuracy of about 5%. If the photon interaction is in lung heterogeneity, these corrections are not as effective as other cases. Nevertheless, it still has acceptable accuracy and works better than the non-corrected scaling method. In order to investigate the importance of the corrections made, it is recommended that their effects on the (S/C) dose calculation algorithm be examined.

References

1. Papanikolaou N, Mackie T R, Meger-Wells C, Gehring M, Reckwerdt P. Investigation of the convolution method for polyenergetic spectra. Medical Physics. 1993; 20(5):1327-36. DOI | PubMed
2. Huang J Y, Eklund D, Childress N L, Howell R M, Mirkovic D, Followill D S, et al. Investigation of various energy deposition kernel refinements for the convolution/superposition method. Medical Physics. 2013; 40(12):121-7. Publisher Full Text | DOI | PubMed
3. Hoban P, Murray D, Round W. Photon beam convolution using polyenergetic energy deposition kernels. Physics in Medicine & Biology. 1994; 39(4):669-85.
4. Pyyry J. Convolution and model-based dose calculation methods in radionuclide and external-beam photon therapy [dissertation]. Helsinki: Department of Neuroscience and Biomedical Engineering; 2016. Available from:https://aaltodoc.aalto.fi/bitstream/handle/123456789/20076/isbn9789526067285.pdf?isAllowed=y&sequence=1.
5. Makrani D S, Hasanzadeh H, Pourfallah T A, Ghasemi A, Jadidi M, Babapour H. Determination of primary electron beam parameters in a Siemens Primus Linac using Monte Carlo simulation. Journal of Paramedical Sciences (JPS). 2015; 6(1):75-9.
6. Tillikainen L, Helminen H, Torsti T, Siljamäki S, Alakuijala J, Pyyry J, et al. A 3D pencil-beam-based superposition algorithm for photon dose calculation in heterogeneous media. Physics in Medicine & Biology. 2008; 53(14):3821-39. DOI | PubMed
7. Mackie T, Bielajew A, Rogers D, Battista J. Generation of photon energy deposition kernels using the EGS Monte Carlo code. Physics in Medicine & Biology. 1988; 33(1):1-20.
8. Ahnesjö A, Saxner M. Characterization of Photon Beams for Kernel-Based Dose Calculation Methods. InTumor Response Monitoring and Treatment Planning; Berlin, Heidelberg: Springer; 1992. p. 479-86.
9. Mackie T, Scrimger J, Battista J. A convolution method of calculating dose for 15-MV x rays. Medical Physics. 1985; 12(2):188-96. DOI | PubMed
10. Bielajew A F, Hirayama H, Nelson W R, Rogers D W O. History, overview and recent improvement of egs4. Report NRC-PIRS-0436; Stanford, California, USA: Stanford Linear Accelerator Center (SLAC); 1994.
11. Rogers D, Seuntjens J, Walters B, Mainegra-Hing E. NRC user codes for EGSnrc. Technical Report PIRS-702 (RevB); Canada: NRC; 2011.
12. Rogers D, Kawrakow I, Seuntjens J, Walters B, Mainegra-Hing E. NRC User Codes for EGSnrc; NRCC Report PIRS-702 (revB); Canada: NRC; 2010.
13. American Brain Tumor Association. About brain tumors, brain tumor education. 2020. Available from: https://www.abta.org/about-brain-tumors/ brain-tumor-education/.
14. O’Connor J. The density scaling theorem applied to lateral electronic equilibrium. Medical Physics. 1984; 11(5):678-80. DOI | PubMed
15. Bjärngard B E. On Fano’s and O’Connor’s theorems. Radiation Research. 1987; 109(2):184-9. DOI | PubMed
16. Ahnesjö A. Collapsed cone convolution of radiant energy for photon dose calculation in heterogeneous media. Medical Physics. 1989; 16(4):577-92. DOI | PubMed
17. O’connor J. The variation of scattered x-rays with density in an irradiated body. Physics in Medicine & Biology. 1957; 1(4):352-69. DOI | PubMed
18. Mohan R, Chui C, Lidofsky L. Differential pencil beam dose computation model for photons. Medical Physics. 1986; 13(1):64-73. DOI | PubMed