- Email: [email protected]
A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation
Radiation Physics and Chemistry ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Contents lists available at ScienceDirect
Radiation Physics and Chemistry journal homepage: www.elsevier.com/locate/radphyschem
A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation O. Koncek a,n, J. Krivonoska b a b
Department of Dosimetry and Application of Ionizing Radiation, Faculty of Nuclear Sciences and Physical Engineering, CTU, Prague, Czech Republic UJP PRAHA, Nad Kaminkou 1345, Prague, Zbraslav, Czech Republic
H I G H L I G H T S
Monte Carlo simulation of photon and electron spectra for a 60Co therapy unit. Monte Carlo simulation of a polyenergetic Pencil Beam Kernel (PBK). Inhomogeneity correction based on Monte Carlo simulation. Implementation of dose calculation algorithm on the CT matrix with 3D inhomogeneity correction.
art ic l e i nf o
a b s t r a c t
Article history: Received 22 June 2013 Accepted 19 March 2014
The MCNP Monte Carlo code was used to simulate the collimating system of the 60Co therapy unit to calculate the primary and scattered photon fluences as well as the electron contamination incident to the isocentric plane as the functions of the irradiation field size. Furthermore, a Monte Carlo simulation for the polyenergetic Pencil Beam Kernels (PBKs) generation was performed using the calculated photon and electron spectra. The PBK was analytically fitted to speed up the dose calculation using the convolution technique in the homogeneous media. The quality of the PBK fit was verified by comparing the calculated and simulated 60Co broad beam profiles and depth dose curves in a homogeneous water medium. The inhomogeneity correction coefficients were derived from the PBK simulation of an inhomogeneous slab phantom consisting of various materials. The inhomogeneity calculation model is based on the changes in the PBK radial displacement and on the change of the forward and backward electron scattering. The inhomogeneity correction is derived from the electron density values gained from a complete 3D CT array and considers different electron densities through which the pencil beam is propagated as well as the electron density values located between the interaction point and the point of dose deposition. Important aspects and details of the algorithm implementation are also described in this study. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Dose calculation Pencil beam Monte Carlo simulation Inhomogeneity correction
1. Introduction Dose calculation algorithms are essential for the treatment plan preparation in radiation therapy. The goal of these algorithms is to accurately predict the dose delivered to the patient during radiation therapy in order to decide if a given treatment plan is clinically acceptable. The most accurate dose prediction method considered is the Monte Carlo simulation and even though some treatment planning systems provide this method (Fix et al., 2010; Fippel, 1999), the calculations are still very time consuming. Due to
n
Corresponding author. Tel.: þ 420 725 476 230. E-mail address: [email protected] (O. Koncek).
the time demanding nature of these calculations, the dose engines based on the convolution techniques are still widely used for radiotherapy treatment planning. In the last decade efforts were made to make these algorithms either faster (Gu et al., 2009 and 2011; Hissoiny et al., 2009) or more accurate (Lu and Chen, 2010; Tillikainen et al., 2008) as documented in the scientific literature. An integral part of the development of these algorithms is verification of their accuracy under clinical conditions. The accuracy can be verified by comparisons of the calculation with measurements in a phantom (Gray et al., 2009; Xu et al., 2009; Breitman et al., 2007; Gagné and Zavgorodni, 2007; Esch et al., 2006; Bragg and Conway, 2006) or with the results of Monte Carlo simulations (Sterpin et al., 2007; Krieger and Sauer, 2005; Chow et al., 2003).
http://dx.doi.org/10.1016/j.radphyschem.2014.03.019 0969-806X/& 2014 Elsevier Ltd. All rights reserved.
Please cite this article as: Koncek, O., Krivonoska, J., A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation. Radiat. Phys. Chem. (2014), http://dx.doi.org/10.1016/j.radphyschem.2014.03.019i
O. Koncek, J. Krivonoska / Radiation Physics and Chemistry ∎ (∎∎∎∎) ∎∎∎–∎∎∎
2
bi 40. The value of the parameter n depends on the type of the simulated particle (n ¼8 for photon PBK, n¼ 4 for electron PBK) !!2 n r j2 ai ðzÞ exp min Hða; bÞ ¼ ∑ PBKnorm ðr j ; zÞ ∑ ð2Þ 2 2 bi ðzÞ j i ¼ 1 π bi ðzÞ Fitting the PBK using the Gaussian curves is a common practice for electrons (Hogstrom et al., 1981) but can also be used for photons (Ahnesjo, 1992; Ulmer and Nardet, 1995; Ulmer et al., 2004). The advantage of such fitting procedure is that the calculation of the dose convolution integral (Eq. (3), φ denotes fluence) can be performed analytically—Eq. (4), for parallel beam; parameters u1, u2, v1, and v2 denote field size defined by all four jaws, x mat=w and y denote lateral positions, z is depth, ρel is the relative electron density of material (where index mat means material, index w relates to water), SSD is the source to surface distance, and erf is the error function: Dðx; y; zÞ ¼ ∬ φðu; vÞPBKðx u; y v; zÞdu dv Aðz; ρ
mat=w
Fig. 1. Radial profiles of the photon pencil beam kernels obtained by the Monte Carlo simulation for three different depths; RSD…relative standard deviation.
Þφ
el Dðx; y; zÞ ¼ 4ð1 þ z=SSDÞ 2 "
bi ðz; ρel ! y v2
" 2.1. Treatment head model
erf 60
A parametric model of a Co unit head, manufactured by UJP PRAHA, was created in the Matlab code (the Mathworks, USA). The model includes a detailed description of the radioactive capsule and the movable primary and secondary collimation systems. After entering the desired field size, the Matlab routine calculates the necessary transformations of the model and generates a file used for the MCNP5 Monte Carlo simulation. The correctness of the model was verified by comparing the calculated and measured beam profiles and depth dose curves in a water phantom. A simulation of photon and electron spectra in the air at the distance of 80 cm below the source was performed for field sizes from 5 5 up to 36 36 cm2 (minimum and maximum field sizes). Fluences were scored on the beam axis and also off-axis in the x and y directions with a step-size of 1 cm. 2.2. Pencil Beam Kernel simulation Pencil Beam Kernel represents the dose contribution of a very narrow beam incident to a water-like medium. Using the calculated photon and electron spectra, separate PBK simulations for photons and electrons were performed. Energy deposition (tally *F8) of a point monodirectional source was scored in a homogeneous water phantom having radius and depth of 60 cm. The phantom was divided by horizontal planes and coaxial cylinders of various radii (with minimum spacing of 0.01 cm immediately surrounding the coordinate system origin). The total of 5 108 particles was simulated with the energy cut-off for photon and electron transport set to 1 keV and 5 keV, respectively. The photon and electron PBKs were normalized (Ahnesjo and Aspradakis, 1999, Ulmer and Nardet, 1995) for each depth using the area integral A(z): ∬ PBKðx; y; zÞdxdy ¼ AðzÞ
ð1Þ
Normalized photon PBK profiles are shown in Fig. 1. In the next step, the normalized PBK (PBKnorm ¼ PBK/A(z)) values were fitted by a sum of Gaussian curves, using a routine written in Matlab. This routine minimizes the difference between the simulated PBK and the sum of Gaussians with respect to the set of parameters (amplitude) a and (radial displacement) b (see Eq. (2)) and satisfies the compulsory conditions Σai ¼1, ai Z0, and
!
x u2
∑ai ðzÞ erf
2. Experimental
ð3Þ
mat=w
bi ðz; ρel
mat=w
Þ
Þ
þ erf
x u1
erf
bi ðz; ρel !# y v1
mat=w
bi ðz; ρel
mat=w
Þ
!# Þ ð4Þ
Moreover the error function can also be expressed analytically. In this calculation, the fit published by Gautschi (1972) with maximum deviation of 1.5 10 7 was used. The normalization parameter A(z) [Gy cm3/photon] according to the depth z was analytically described as AðzÞ ¼ Amax f½p1 þ p2 ðp3 zÞp4 expð p5 zÞ p6 expð p7 z2 Þ g 13
ð5Þ
3
[Gy cm /photon], p1 ¼1.007, p2 ¼0.225, where Amax ¼4.4313 10 p3 ¼0.2218 cm 1, p4 ¼ 0.9667, p5 ¼0.06064 cm 1, p6 ¼ 0.576, and p7 ¼26.93 cm 2. 2.3. Inhomogeneities In order to account for the influence of inhomogeneities on the calculation a similar PBK simulation, as described in Section 2.2, was performed using various inhomogeneity layers (lung, bone, adipose, skeleton, etc.) in the water phantom. The material composition was adopted from Cassell et al. (1981) and McConn et al. (2011). This simulation was focused primarily on the effects of the PBK shape and amplitude changes at the interfaces of the different density materials. 2.4. Superposition computational method The computational algorithm considers a fan line beam impacting on a rectangular CT matrix. The coordinates in the dose matrix are calculated due to the beam divergence as x' ¼x(1 þz/SSD), y' ¼y (1 þz/SSD) and z'¼z[(x2 þ y2)/SSD2 þ1]1/2. In order to calculate the dose distribution, the intersection points of fan beams and the rectangular 3D CT matrix have to be determined. To calculate the intersection points of one ray, the calculation of the intersections with all equidistant x-, y-, and z-planes was performed. The final order of the intersection points with the given planes (as the beam traverses the 3D matrix) is reached by arranging the individual values (a group of the intersections with x-, y-, and z-planes) according to their distance from the source. The distance between two neighboring intersection points gives the length of the beam trajectory in the voxel. Intersections are considered separately for each plane; therefore
Please cite this article as: Koncek, O., Krivonoska, J., A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation. Radiat. Phys. Chem. (2014), http://dx.doi.org/10.1016/j.radphyschem.2014.03.019i
O. Koncek, J. Krivonoska / Radiation Physics and Chemistry ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Fig. 2. Relative photon energy spectra on the central axis reaching the scoring plane at distance 80 cm for five different field sizes. Spectra are normalized to the total photon fluence for each field size.
Fig. 3. Relative photon fluence versus field size of various components reaching a plane at distance 80 cm. The number of photos is normalized to the total photon fluence for the maximum field size.
the value of the electron density of the voxel gained from the electron density matrix can be directly matched to the chord length of this voxel, since only the changes in x, y and z coordinates caused by the traversing beam in a voxel matrix can be registered. In the next step, a similar calculation process is applied, though in the calculation plane, and the result is used to calculate the radial broadening of the PBK between the interaction point and the calculation point.
3. Results and discussion The relative photon fluence on the central axis obtained from the Monte Carlo simulation is shown in Fig. 2. The total photon fluence for the field size of 10 10 cm2 is 1.256 10 5 cm 2/history. Fig. 3 shows the relative photon fluence on the central axis for different photon components reaching the isocentric plane. The number of primary photons remains relatively constant as the field size increases (from 58.5% to 58.7% of the total fluence). The contribution from the adjustable collimator (primary and
3
Fig. 4. Relative electron energy spectra (normalized to the total electron fluence for each field size) on the central axis reaching the scoring plane at 80 cm for five different field sizes. Mean energy is 0.643 MeV for field 5×5 cm2, 0.606 MeV for field 10×10 cm2, 0.564 MeV for field 20×20 cm2, 0.537 MeV for field 30×30 cm2, and 0.532 MeV for field 36×36 cm2.
Fig. 5. Variation in electron contamination with field size. The red line represents a linear fit described by the equation: y ¼ 0.0295*x – 0.0431; R2 ¼ 0.998.
secondary) to the total fluence lies in the interval from 1.5% for the field size of 4 4 cm2 to 7.8% for the field size of 36 36 cm2. The changes in the photon spectrum with regard to the off-axis distance were neglected for the PBK simulation (the maximum deviation of 1.3% was observed for the field size of 36 36 cm2 in the 16 cm distance from the central axis and for the energy of 1.332 MeV). The electron spectra for different field sizes are shown in Fig. 4. The mean electron energy decreases with increasing field size from 0.643 MeV to 0.532 MeV. This is due to the fact that the major source of electrons for small field sizes is the capsule, while for larger fields the collimator system becomes dominant. The fluence of electron contamination increases with increasing field size by a factor of 9 from the smallest to the largest field (as presented in Fig. 5). For the maximum field size, the total electron fluence represents only 0.66% of the total photon fluence, but it could not be neglected in dose calculation because the dose delivered per unit fluence of electrons is significantly higher than that of photons.
Please cite this article as: Koncek, O., Krivonoska, J., A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation. Radiat. Phys. Chem. (2014), http://dx.doi.org/10.1016/j.radphyschem.2014.03.019i
4
O. Koncek, J. Krivonoska / Radiation Physics and Chemistry ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Fig. 7. Result of 2D gamma evaluation (ΔD¼ 3%, DTA ¼ 3 mm) between the Monte Carlo simulation and the convolution algorithm with the GBPL inhomogeneity correction.
Fig. 6. Depth dose curve in a phantom with slab inhomogeneity between 1cm o z o 3cm with relative electron density 0.288. Red samples refer to Monte Carlo simulation (error bar represents 1 standard deviation), solid curves to calculation with superposition (blue) or convolution algorithm with GBPL correction (black). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
The PBK, which was derived in a water medium, is converted for the use in the inhomogeneous medium according the relative electron density of the material Ahnesjo and Aspradakis, 1999 and mat=w mat w the lateral PBK displacement changes according to bi ¼bi ρel , mat=w where ρel is the relative electron density determined from the CT number. The depth parameter z from Eq. (5) has to be adjusted for other materials as z(ρ)¼z ρq, where ρ is the relative electron density and the value of the parameter q varies according to the type of tissue (e.g. 1.57 for bone tissue and 6.15 for lung tissue). The shape of the radial profile of the PBK due to the changes in backward and forward scattering near the inhomogeneity layer can be omitted (for 60Co beam) and the correction is proportional only to the difference of A(z) on the inhomogeneity border. The decrease (or increase) in forward or backward scattering along the beam line can be expressed as an exponential function δA(δz, ρ) exp( ρ μ δz) where δz is the distance between the calculation point and the boundary of the inhomogeneous material, ρ is the relative electron density, and μ is the linear attenuation coefficient. To test the accuracy of the superposition algorithm, the calculated dose profiles and depth dose curves were compared to the Monte Carlo simulation of a broad parallel beam incident to a water phantom with various inhomogeneities. Moreover, the Generalized Batho Power Law (GBPL) correction (El-Khatib, and Battista, 1983; Siddon, 1984; Sonntag and Cunningham, 1977; Thomas 1991; Wong and Henkelman, 1982) was implemented to the convolution algorithm. The comparison of the dose calculation (for the case of lung tissue) is shown in Fig. 6. The maximum deviation of the superposition algorithm is 0.51%, of the GBPL 0.55%. Therefore, in this case there is no significant difference between these two algorithms, although the GBPL represents only a 1D correction not considering the electron transport.
Finally, the scan2mcnp (White Rock Science, USA) program was used to create a voxel phantom using the CT images. The Monte Carlo simulation of the dose distribution in the voxel phantom (including the treatment head model) was performed. This distribution was compared to the dose calculated with the superposition and convolution (with the GBPL correction) algorithms. The comparison was performed using the 2D gamma evaluation (with the criteria ΔD¼3%, DTA ¼3 mm) proposed by Low et al. (1998); see Fig. 7. The gamma value for the superposition algorithm is less than 1 in all pixels (with maximum dose deviation of 2.1%) whereas for the convolution algorithm γ 4 1 in 1.5% of all pixels (from 512 512 matrix). Voxels from the voxel phantom with low statistics (relative standard deviation worse than 5%) were excluded from the comparison. 4. Conclusion In this paper the design and implementation of the superposition algorithm were described based on the Pencil Beam Kernel for the 60Co source and the 3D inhomogeneity correction derived from various Monte Carlo simulations of the PBK on the inhomogeneity borders. The presented algorithm is believed to enhance the accuracy of the calculated dose distributions and although it is primarily derived for the cobalt beam, some principles can be used in future calculations for the photon megavoltage beams. References Ahnesjo, A., 1992. A pencil beam model for photon dose calculation. Med. Phys. 19, 263–273. Ahnesjo, A., Aspradakis, M.M., 1999. Dose calculation for external photon beams in radiotherapy. Phys. Med. Biol. 44, R99–R155. Bragg, C.M., Conway, J., 2006. Dosimetric verification of the anisotropic analytical algorithm for radiotherapy treatment planning. Ratiother. Oncol. 81, 315–323. Breitman, K., Rathee, S., Newcomb, C., Murray, B., Robinson, D., et al., 2007. Experimental validation of the eclipse AAA algorithm. J. Appl. Clin. Med. Phys. 8, 76–92. Cassell, K.J., Hobday, P.A., Parker, R.P., 1981. The implementation of a generalised Batho inhomogeneity correction for radiotherapy planning with direct use of CT numbers. Phys. Med. Biol. 26, 825–833. Chow, J.C.L., Wong, E., Chen, J.Z., Dyk, J.V., 2003. Comparison of dose calculation algorithms with Monte Carlo methods for photon arcs. Med. Phys. 30, 2686–2694. El-Khatib, E., Battista, J.J., 1983. Improved lung dose calculation using tissiuemaximum ratios in the Batho correction. Med. Phys. 11, 279–286.
Please cite this article as: Koncek, O., Krivonoska, J., A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation. Radiat. Phys. Chem. (2014), http://dx.doi.org/10.1016/j.radphyschem.2014.03.019i
O. Koncek, J. Krivonoska / Radiation Physics and Chemistry ∎ (∎∎∎∎) ∎∎∎–∎∎∎ Esch, A.V., Tillikainen, L., Pyykkonen, J., Tenhuen, M., et al., 2006. Testing of the analytical anisotropic algorithm for photon dose calculation. Med. Phys. 33, 4130–4148. Fippel, M., 1999. Fast Monte Carlo dose calculation for photon beams based on the VMC electron algorithm. Med. Phys. 26, 1466–1475. Fix, M.K., Frei, D., Volken, W., Neuenschwander, H., Born, E.J., Manser, P., 2010. Monte Carlo dose calculation improvements for low energy electron beams using eMC. Phys. Med. Biol. 55, 4577–4588. Gagné, I.M., Zavgorodni, S., 2007. Evaluation of the analytical anisotropic algorithm in an extreme water–lung interface phantom using Monte Carlo dose calculations. J. Appl. Clin. Med. Phys. 8, 33–46. Gautschi, W., 1972. Error Function and Fresnel Integrals. In: Abramowitz, M., Stegun, I.A. (Eds.), Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover Publications, Mineola, pp. 295–330. Gray, A., Oliver, L.D., Johnston, P.N., 2009. The accuracy of the pencil beam convolution and anisotropic analytical algorithms in predicting the dose effects due to attenuation from immobilization devices and large air gaps. Med. Phys. 36, 3181–3191. Gu, X., Choi, D., Men, C., Pan, H., Majumdar, A., Jiang, S.B., 2009. GPU-based ultrafast dose calculation using a finite size pencil beam model. Phys. Med. Biol. 54, 6287–6297. Gu, X., Jelen, U., Li, J., Jia, X., Jiang, S.B., 2011. A GPU-based finite-size pencil beam algorithm with 3D-density correction for radiotherapy dose calculation. Phys. Med. Biol. 56, 3337–3350. Hissoiny, S., Ozell, B., Després, P., 2009. Fast convolution-superposition dose calculation on graphics hardware. Med. Phys. 36, 1998–2005. Hogstrom, K.R., Mills, M.D., Almond, P.R., 1981. Electron beam dose calculation. Phys. Med. Biol. 26, 445–459. Krieger, T., Sauer, O.A., 2005. Monte Carlo-versus pencil-beam-/collapsed-conedose calculation in a heterogeneous multi-layer phantom. Phys. Med. Biol. 50, 859–868.
5
Low, D.A., Harms, W.B., Mutic, S., Purdy, A., 1998. A technique for the quantitative evaluation of dose distribution. Med. Phys. 25, 656–661. Lu, W., Chen, M., 2010. Fluence-convolution broad-beam (FCBB) dose calculation. Phys. Med. Biol. 55, 7211–7229. McConn Jr., R.J., Gesh, C.J., Pagh, R.T., Rucker, R.A., Williams III, R., 2011. Compendium of Material Composition Data for Radiation Transport Modeling. Pacific Northwest National Laboratory, Richland, WA (PNNL-15870 Rev. 1). Siddon, R.L., 1984. Generalised Batho correction factor. Phys. Med. Biol. 29, 1575–1579. Sonntag, M.R., Cunningham, J.R., 1977. Corrections to absorbed dose calculations for tissue inhomogeneities. Med. Phys. 4, 431–436. Sterpin, E., Tomsej, M., Smedt, B.D., Reynaert, N., Vynckier, S., 2007. Monte Carlo evaluation of the AAA treatment planning algorithm in a heterogeneous multilayer phantom and IMRT clinical treatments for an Elekta SL25 linear akcelerátor. Med. Phys. 34, 1665–1677. Thomas, S.J., 1991. A modified power-law formula for inhomogeneity corrections in beams of high-energy x rays. Med. Phys. 18, 719–723. Tillikainen, L., Helminen, H., Torsti, T., Siljamäki, S., Alakuijala, J., Pyyry, J., Ulmer, W., 2008. Phys. Med. Biol. 53, 3821–3839. Ulmer, W., Nardet, D., 1995. A tripple Gaussian pencil beam model for photon beam treatment planning. Z. Med. Phys. 5, 25–30. Ulmer, W., Pyyry, J., Kaissl, W., 2004. A 3D photon superposition/convolution algorithm and its foundation on results of Monte Carlo calculations. Phys. Med. Biol. 50, 1767–1790. Wong, J.W., Henkelman, R.M., 1982. Reconsideration of the power-law (Batho) equation for inhomogeneity corrections. Med. Phys. 9, 521–530. Xu, Z., Walsh, S.E., Telivala, T.P., Meek, A.G., Yang, G., 2009. Evaluation of the eclipse electron Monte Carlo dose calculation for small fields. J. Appl. Clin. Med. Phys. 10, 75–85.
Please cite this article as: Koncek, O., Krivonoska, J., A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation. Radiat. Phys. Chem. (2014), http://dx.doi.org/10.1016/j.radphyschem.2014.03.019i
Contents lists available at ScienceDirect
Radiation Physics and Chemistry journal homepage: www.elsevier.com/locate/radphyschem
A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation O. Koncek a,n, J. Krivonoska b a b
Department of Dosimetry and Application of Ionizing Radiation, Faculty of Nuclear Sciences and Physical Engineering, CTU, Prague, Czech Republic UJP PRAHA, Nad Kaminkou 1345, Prague, Zbraslav, Czech Republic
H I G H L I G H T S
Monte Carlo simulation of photon and electron spectra for a 60Co therapy unit. Monte Carlo simulation of a polyenergetic Pencil Beam Kernel (PBK). Inhomogeneity correction based on Monte Carlo simulation. Implementation of dose calculation algorithm on the CT matrix with 3D inhomogeneity correction.
art ic l e i nf o
a b s t r a c t
Article history: Received 22 June 2013 Accepted 19 March 2014
The MCNP Monte Carlo code was used to simulate the collimating system of the 60Co therapy unit to calculate the primary and scattered photon fluences as well as the electron contamination incident to the isocentric plane as the functions of the irradiation field size. Furthermore, a Monte Carlo simulation for the polyenergetic Pencil Beam Kernels (PBKs) generation was performed using the calculated photon and electron spectra. The PBK was analytically fitted to speed up the dose calculation using the convolution technique in the homogeneous media. The quality of the PBK fit was verified by comparing the calculated and simulated 60Co broad beam profiles and depth dose curves in a homogeneous water medium. The inhomogeneity correction coefficients were derived from the PBK simulation of an inhomogeneous slab phantom consisting of various materials. The inhomogeneity calculation model is based on the changes in the PBK radial displacement and on the change of the forward and backward electron scattering. The inhomogeneity correction is derived from the electron density values gained from a complete 3D CT array and considers different electron densities through which the pencil beam is propagated as well as the electron density values located between the interaction point and the point of dose deposition. Important aspects and details of the algorithm implementation are also described in this study. & 2014 Elsevier Ltd. All rights reserved.
Keywords: Dose calculation Pencil beam Monte Carlo simulation Inhomogeneity correction
1. Introduction Dose calculation algorithms are essential for the treatment plan preparation in radiation therapy. The goal of these algorithms is to accurately predict the dose delivered to the patient during radiation therapy in order to decide if a given treatment plan is clinically acceptable. The most accurate dose prediction method considered is the Monte Carlo simulation and even though some treatment planning systems provide this method (Fix et al., 2010; Fippel, 1999), the calculations are still very time consuming. Due to
n
Corresponding author. Tel.: þ 420 725 476 230. E-mail address: [email protected] (O. Koncek).
the time demanding nature of these calculations, the dose engines based on the convolution techniques are still widely used for radiotherapy treatment planning. In the last decade efforts were made to make these algorithms either faster (Gu et al., 2009 and 2011; Hissoiny et al., 2009) or more accurate (Lu and Chen, 2010; Tillikainen et al., 2008) as documented in the scientific literature. An integral part of the development of these algorithms is verification of their accuracy under clinical conditions. The accuracy can be verified by comparisons of the calculation with measurements in a phantom (Gray et al., 2009; Xu et al., 2009; Breitman et al., 2007; Gagné and Zavgorodni, 2007; Esch et al., 2006; Bragg and Conway, 2006) or with the results of Monte Carlo simulations (Sterpin et al., 2007; Krieger and Sauer, 2005; Chow et al., 2003).
http://dx.doi.org/10.1016/j.radphyschem.2014.03.019 0969-806X/& 2014 Elsevier Ltd. All rights reserved.
Please cite this article as: Koncek, O., Krivonoska, J., A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation. Radiat. Phys. Chem. (2014), http://dx.doi.org/10.1016/j.radphyschem.2014.03.019i
O. Koncek, J. Krivonoska / Radiation Physics and Chemistry ∎ (∎∎∎∎) ∎∎∎–∎∎∎
2
bi 40. The value of the parameter n depends on the type of the simulated particle (n ¼8 for photon PBK, n¼ 4 for electron PBK) !!2 n r j2 ai ðzÞ exp min Hða; bÞ ¼ ∑ PBKnorm ðr j ; zÞ ∑ ð2Þ 2 2 bi ðzÞ j i ¼ 1 π bi ðzÞ Fitting the PBK using the Gaussian curves is a common practice for electrons (Hogstrom et al., 1981) but can also be used for photons (Ahnesjo, 1992; Ulmer and Nardet, 1995; Ulmer et al., 2004). The advantage of such fitting procedure is that the calculation of the dose convolution integral (Eq. (3), φ denotes fluence) can be performed analytically—Eq. (4), for parallel beam; parameters u1, u2, v1, and v2 denote field size defined by all four jaws, x mat=w and y denote lateral positions, z is depth, ρel is the relative electron density of material (where index mat means material, index w relates to water), SSD is the source to surface distance, and erf is the error function: Dðx; y; zÞ ¼ ∬ φðu; vÞPBKðx u; y v; zÞdu dv Aðz; ρ
mat=w
Fig. 1. Radial profiles of the photon pencil beam kernels obtained by the Monte Carlo simulation for three different depths; RSD…relative standard deviation.
Þφ
el Dðx; y; zÞ ¼ 4ð1 þ z=SSDÞ 2 "
bi ðz; ρel ! y v2
" 2.1. Treatment head model
erf 60
A parametric model of a Co unit head, manufactured by UJP PRAHA, was created in the Matlab code (the Mathworks, USA). The model includes a detailed description of the radioactive capsule and the movable primary and secondary collimation systems. After entering the desired field size, the Matlab routine calculates the necessary transformations of the model and generates a file used for the MCNP5 Monte Carlo simulation. The correctness of the model was verified by comparing the calculated and measured beam profiles and depth dose curves in a water phantom. A simulation of photon and electron spectra in the air at the distance of 80 cm below the source was performed for field sizes from 5 5 up to 36 36 cm2 (minimum and maximum field sizes). Fluences were scored on the beam axis and also off-axis in the x and y directions with a step-size of 1 cm. 2.2. Pencil Beam Kernel simulation Pencil Beam Kernel represents the dose contribution of a very narrow beam incident to a water-like medium. Using the calculated photon and electron spectra, separate PBK simulations for photons and electrons were performed. Energy deposition (tally *F8) of a point monodirectional source was scored in a homogeneous water phantom having radius and depth of 60 cm. The phantom was divided by horizontal planes and coaxial cylinders of various radii (with minimum spacing of 0.01 cm immediately surrounding the coordinate system origin). The total of 5 108 particles was simulated with the energy cut-off for photon and electron transport set to 1 keV and 5 keV, respectively. The photon and electron PBKs were normalized (Ahnesjo and Aspradakis, 1999, Ulmer and Nardet, 1995) for each depth using the area integral A(z): ∬ PBKðx; y; zÞdxdy ¼ AðzÞ
ð1Þ
Normalized photon PBK profiles are shown in Fig. 1. In the next step, the normalized PBK (PBKnorm ¼ PBK/A(z)) values were fitted by a sum of Gaussian curves, using a routine written in Matlab. This routine minimizes the difference between the simulated PBK and the sum of Gaussians with respect to the set of parameters (amplitude) a and (radial displacement) b (see Eq. (2)) and satisfies the compulsory conditions Σai ¼1, ai Z0, and
!
x u2
∑ai ðzÞ erf
2. Experimental
ð3Þ
mat=w
bi ðz; ρel
mat=w
Þ
Þ
þ erf
x u1
erf
bi ðz; ρel !# y v1
mat=w
bi ðz; ρel
mat=w
Þ
!# Þ ð4Þ
Moreover the error function can also be expressed analytically. In this calculation, the fit published by Gautschi (1972) with maximum deviation of 1.5 10 7 was used. The normalization parameter A(z) [Gy cm3/photon] according to the depth z was analytically described as AðzÞ ¼ Amax f½p1 þ p2 ðp3 zÞp4 expð p5 zÞ p6 expð p7 z2 Þ g 13
ð5Þ
3
[Gy cm /photon], p1 ¼1.007, p2 ¼0.225, where Amax ¼4.4313 10 p3 ¼0.2218 cm 1, p4 ¼ 0.9667, p5 ¼0.06064 cm 1, p6 ¼ 0.576, and p7 ¼26.93 cm 2. 2.3. Inhomogeneities In order to account for the influence of inhomogeneities on the calculation a similar PBK simulation, as described in Section 2.2, was performed using various inhomogeneity layers (lung, bone, adipose, skeleton, etc.) in the water phantom. The material composition was adopted from Cassell et al. (1981) and McConn et al. (2011). This simulation was focused primarily on the effects of the PBK shape and amplitude changes at the interfaces of the different density materials. 2.4. Superposition computational method The computational algorithm considers a fan line beam impacting on a rectangular CT matrix. The coordinates in the dose matrix are calculated due to the beam divergence as x' ¼x(1 þz/SSD), y' ¼y (1 þz/SSD) and z'¼z[(x2 þ y2)/SSD2 þ1]1/2. In order to calculate the dose distribution, the intersection points of fan beams and the rectangular 3D CT matrix have to be determined. To calculate the intersection points of one ray, the calculation of the intersections with all equidistant x-, y-, and z-planes was performed. The final order of the intersection points with the given planes (as the beam traverses the 3D matrix) is reached by arranging the individual values (a group of the intersections with x-, y-, and z-planes) according to their distance from the source. The distance between two neighboring intersection points gives the length of the beam trajectory in the voxel. Intersections are considered separately for each plane; therefore
Please cite this article as: Koncek, O., Krivonoska, J., A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation. Radiat. Phys. Chem. (2014), http://dx.doi.org/10.1016/j.radphyschem.2014.03.019i
O. Koncek, J. Krivonoska / Radiation Physics and Chemistry ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Fig. 2. Relative photon energy spectra on the central axis reaching the scoring plane at distance 80 cm for five different field sizes. Spectra are normalized to the total photon fluence for each field size.
Fig. 3. Relative photon fluence versus field size of various components reaching a plane at distance 80 cm. The number of photos is normalized to the total photon fluence for the maximum field size.
the value of the electron density of the voxel gained from the electron density matrix can be directly matched to the chord length of this voxel, since only the changes in x, y and z coordinates caused by the traversing beam in a voxel matrix can be registered. In the next step, a similar calculation process is applied, though in the calculation plane, and the result is used to calculate the radial broadening of the PBK between the interaction point and the calculation point.
3. Results and discussion The relative photon fluence on the central axis obtained from the Monte Carlo simulation is shown in Fig. 2. The total photon fluence for the field size of 10 10 cm2 is 1.256 10 5 cm 2/history. Fig. 3 shows the relative photon fluence on the central axis for different photon components reaching the isocentric plane. The number of primary photons remains relatively constant as the field size increases (from 58.5% to 58.7% of the total fluence). The contribution from the adjustable collimator (primary and
3
Fig. 4. Relative electron energy spectra (normalized to the total electron fluence for each field size) on the central axis reaching the scoring plane at 80 cm for five different field sizes. Mean energy is 0.643 MeV for field 5×5 cm2, 0.606 MeV for field 10×10 cm2, 0.564 MeV for field 20×20 cm2, 0.537 MeV for field 30×30 cm2, and 0.532 MeV for field 36×36 cm2.
Fig. 5. Variation in electron contamination with field size. The red line represents a linear fit described by the equation: y ¼ 0.0295*x – 0.0431; R2 ¼ 0.998.
secondary) to the total fluence lies in the interval from 1.5% for the field size of 4 4 cm2 to 7.8% for the field size of 36 36 cm2. The changes in the photon spectrum with regard to the off-axis distance were neglected for the PBK simulation (the maximum deviation of 1.3% was observed for the field size of 36 36 cm2 in the 16 cm distance from the central axis and for the energy of 1.332 MeV). The electron spectra for different field sizes are shown in Fig. 4. The mean electron energy decreases with increasing field size from 0.643 MeV to 0.532 MeV. This is due to the fact that the major source of electrons for small field sizes is the capsule, while for larger fields the collimator system becomes dominant. The fluence of electron contamination increases with increasing field size by a factor of 9 from the smallest to the largest field (as presented in Fig. 5). For the maximum field size, the total electron fluence represents only 0.66% of the total photon fluence, but it could not be neglected in dose calculation because the dose delivered per unit fluence of electrons is significantly higher than that of photons.
Please cite this article as: Koncek, O., Krivonoska, J., A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation. Radiat. Phys. Chem. (2014), http://dx.doi.org/10.1016/j.radphyschem.2014.03.019i
4
O. Koncek, J. Krivonoska / Radiation Physics and Chemistry ∎ (∎∎∎∎) ∎∎∎–∎∎∎
Fig. 7. Result of 2D gamma evaluation (ΔD¼ 3%, DTA ¼ 3 mm) between the Monte Carlo simulation and the convolution algorithm with the GBPL inhomogeneity correction.
Fig. 6. Depth dose curve in a phantom with slab inhomogeneity between 1cm o z o 3cm with relative electron density 0.288. Red samples refer to Monte Carlo simulation (error bar represents 1 standard deviation), solid curves to calculation with superposition (blue) or convolution algorithm with GBPL correction (black). (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)
The PBK, which was derived in a water medium, is converted for the use in the inhomogeneous medium according the relative electron density of the material Ahnesjo and Aspradakis, 1999 and mat=w mat w the lateral PBK displacement changes according to bi ¼bi ρel , mat=w where ρel is the relative electron density determined from the CT number. The depth parameter z from Eq. (5) has to be adjusted for other materials as z(ρ)¼z ρq, where ρ is the relative electron density and the value of the parameter q varies according to the type of tissue (e.g. 1.57 for bone tissue and 6.15 for lung tissue). The shape of the radial profile of the PBK due to the changes in backward and forward scattering near the inhomogeneity layer can be omitted (for 60Co beam) and the correction is proportional only to the difference of A(z) on the inhomogeneity border. The decrease (or increase) in forward or backward scattering along the beam line can be expressed as an exponential function δA(δz, ρ) exp( ρ μ δz) where δz is the distance between the calculation point and the boundary of the inhomogeneous material, ρ is the relative electron density, and μ is the linear attenuation coefficient. To test the accuracy of the superposition algorithm, the calculated dose profiles and depth dose curves were compared to the Monte Carlo simulation of a broad parallel beam incident to a water phantom with various inhomogeneities. Moreover, the Generalized Batho Power Law (GBPL) correction (El-Khatib, and Battista, 1983; Siddon, 1984; Sonntag and Cunningham, 1977; Thomas 1991; Wong and Henkelman, 1982) was implemented to the convolution algorithm. The comparison of the dose calculation (for the case of lung tissue) is shown in Fig. 6. The maximum deviation of the superposition algorithm is 0.51%, of the GBPL 0.55%. Therefore, in this case there is no significant difference between these two algorithms, although the GBPL represents only a 1D correction not considering the electron transport.
Finally, the scan2mcnp (White Rock Science, USA) program was used to create a voxel phantom using the CT images. The Monte Carlo simulation of the dose distribution in the voxel phantom (including the treatment head model) was performed. This distribution was compared to the dose calculated with the superposition and convolution (with the GBPL correction) algorithms. The comparison was performed using the 2D gamma evaluation (with the criteria ΔD¼3%, DTA ¼3 mm) proposed by Low et al. (1998); see Fig. 7. The gamma value for the superposition algorithm is less than 1 in all pixels (with maximum dose deviation of 2.1%) whereas for the convolution algorithm γ 4 1 in 1.5% of all pixels (from 512 512 matrix). Voxels from the voxel phantom with low statistics (relative standard deviation worse than 5%) were excluded from the comparison. 4. Conclusion In this paper the design and implementation of the superposition algorithm were described based on the Pencil Beam Kernel for the 60Co source and the 3D inhomogeneity correction derived from various Monte Carlo simulations of the PBK on the inhomogeneity borders. The presented algorithm is believed to enhance the accuracy of the calculated dose distributions and although it is primarily derived for the cobalt beam, some principles can be used in future calculations for the photon megavoltage beams. References Ahnesjo, A., 1992. A pencil beam model for photon dose calculation. Med. Phys. 19, 263–273. Ahnesjo, A., Aspradakis, M.M., 1999. Dose calculation for external photon beams in radiotherapy. Phys. Med. Biol. 44, R99–R155. Bragg, C.M., Conway, J., 2006. Dosimetric verification of the anisotropic analytical algorithm for radiotherapy treatment planning. Ratiother. Oncol. 81, 315–323. Breitman, K., Rathee, S., Newcomb, C., Murray, B., Robinson, D., et al., 2007. Experimental validation of the eclipse AAA algorithm. J. Appl. Clin. Med. Phys. 8, 76–92. Cassell, K.J., Hobday, P.A., Parker, R.P., 1981. The implementation of a generalised Batho inhomogeneity correction for radiotherapy planning with direct use of CT numbers. Phys. Med. Biol. 26, 825–833. Chow, J.C.L., Wong, E., Chen, J.Z., Dyk, J.V., 2003. Comparison of dose calculation algorithms with Monte Carlo methods for photon arcs. Med. Phys. 30, 2686–2694. El-Khatib, E., Battista, J.J., 1983. Improved lung dose calculation using tissiuemaximum ratios in the Batho correction. Med. Phys. 11, 279–286.
Please cite this article as: Koncek, O., Krivonoska, J., A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation. Radiat. Phys. Chem. (2014), http://dx.doi.org/10.1016/j.radphyschem.2014.03.019i
O. Koncek, J. Krivonoska / Radiation Physics and Chemistry ∎ (∎∎∎∎) ∎∎∎–∎∎∎ Esch, A.V., Tillikainen, L., Pyykkonen, J., Tenhuen, M., et al., 2006. Testing of the analytical anisotropic algorithm for photon dose calculation. Med. Phys. 33, 4130–4148. Fippel, M., 1999. Fast Monte Carlo dose calculation for photon beams based on the VMC electron algorithm. Med. Phys. 26, 1466–1475. Fix, M.K., Frei, D., Volken, W., Neuenschwander, H., Born, E.J., Manser, P., 2010. Monte Carlo dose calculation improvements for low energy electron beams using eMC. Phys. Med. Biol. 55, 4577–4588. Gagné, I.M., Zavgorodni, S., 2007. Evaluation of the analytical anisotropic algorithm in an extreme water–lung interface phantom using Monte Carlo dose calculations. J. Appl. Clin. Med. Phys. 8, 33–46. Gautschi, W., 1972. Error Function and Fresnel Integrals. In: Abramowitz, M., Stegun, I.A. (Eds.), Handbook of Mathematical Functions with Formulas, Graphs, and Mathematical Tables. Dover Publications, Mineola, pp. 295–330. Gray, A., Oliver, L.D., Johnston, P.N., 2009. The accuracy of the pencil beam convolution and anisotropic analytical algorithms in predicting the dose effects due to attenuation from immobilization devices and large air gaps. Med. Phys. 36, 3181–3191. Gu, X., Choi, D., Men, C., Pan, H., Majumdar, A., Jiang, S.B., 2009. GPU-based ultrafast dose calculation using a finite size pencil beam model. Phys. Med. Biol. 54, 6287–6297. Gu, X., Jelen, U., Li, J., Jia, X., Jiang, S.B., 2011. A GPU-based finite-size pencil beam algorithm with 3D-density correction for radiotherapy dose calculation. Phys. Med. Biol. 56, 3337–3350. Hissoiny, S., Ozell, B., Després, P., 2009. Fast convolution-superposition dose calculation on graphics hardware. Med. Phys. 36, 1998–2005. Hogstrom, K.R., Mills, M.D., Almond, P.R., 1981. Electron beam dose calculation. Phys. Med. Biol. 26, 445–459. Krieger, T., Sauer, O.A., 2005. Monte Carlo-versus pencil-beam-/collapsed-conedose calculation in a heterogeneous multi-layer phantom. Phys. Med. Biol. 50, 859–868.
5
Low, D.A., Harms, W.B., Mutic, S., Purdy, A., 1998. A technique for the quantitative evaluation of dose distribution. Med. Phys. 25, 656–661. Lu, W., Chen, M., 2010. Fluence-convolution broad-beam (FCBB) dose calculation. Phys. Med. Biol. 55, 7211–7229. McConn Jr., R.J., Gesh, C.J., Pagh, R.T., Rucker, R.A., Williams III, R., 2011. Compendium of Material Composition Data for Radiation Transport Modeling. Pacific Northwest National Laboratory, Richland, WA (PNNL-15870 Rev. 1). Siddon, R.L., 1984. Generalised Batho correction factor. Phys. Med. Biol. 29, 1575–1579. Sonntag, M.R., Cunningham, J.R., 1977. Corrections to absorbed dose calculations for tissue inhomogeneities. Med. Phys. 4, 431–436. Sterpin, E., Tomsej, M., Smedt, B.D., Reynaert, N., Vynckier, S., 2007. Monte Carlo evaluation of the AAA treatment planning algorithm in a heterogeneous multilayer phantom and IMRT clinical treatments for an Elekta SL25 linear akcelerátor. Med. Phys. 34, 1665–1677. Thomas, S.J., 1991. A modified power-law formula for inhomogeneity corrections in beams of high-energy x rays. Med. Phys. 18, 719–723. Tillikainen, L., Helminen, H., Torsti, T., Siljamäki, S., Alakuijala, J., Pyyry, J., Ulmer, W., 2008. Phys. Med. Biol. 53, 3821–3839. Ulmer, W., Nardet, D., 1995. A tripple Gaussian pencil beam model for photon beam treatment planning. Z. Med. Phys. 5, 25–30. Ulmer, W., Pyyry, J., Kaissl, W., 2004. A 3D photon superposition/convolution algorithm and its foundation on results of Monte Carlo calculations. Phys. Med. Biol. 50, 1767–1790. Wong, J.W., Henkelman, R.M., 1982. Reconsideration of the power-law (Batho) equation for inhomogeneity corrections. Med. Phys. 9, 521–530. Xu, Z., Walsh, S.E., Telivala, T.P., Meek, A.G., Yang, G., 2009. Evaluation of the eclipse electron Monte Carlo dose calculation for small fields. J. Appl. Clin. Med. Phys. 10, 75–85.
Please cite this article as: Koncek, O., Krivonoska, J., A 3D superposition pencil beam dose calculation algorithm for a 60Co therapy unit and its verification by MC simulation. Radiat. Phys. Chem. (2014), http://dx.doi.org/10.1016/j.radphyschem.2014.03.019i