- Email: [email protected]
3D dose distribution calculation in a voxelized human phantom by means of Monte Carlo method
ARTICLE IN PRESS Applied Radiation and Isotopes 68 (2010) 709–713
Contents lists available at ScienceDirect
Applied Radiation and Isotopes journal homepage: www.elsevier.com/locate/apradiso
3D dose distribution calculation in a voxelized human phantom by means of Monte Carlo method V. Abella, R. Miro´ , B. Juste, G. Verdu´ ´cnica de Vale ncia, Camı´ de Vera, s/n. 46022 Vale ncia, Spain Departamento de Ingenierı´a Quı´mica y Nuclear, Universidad Polite
a r t i c l e in fo
Keywords: Biomedical applications of radiation Zubal phantom MCNP5 Voxelization
abstract The aim of this work is to provide the reconstruction of a real human voxelized phantom by means of a MatLabs program and the simulation of the irradiation of such phantom with the photon beam generated in a Theratron 780s (MDS Nordion) 60Co radiotherapy unit, by using the Monte Carlo transport code MCNP (Monte Carlo N-Particle), version 5. The project results in 3D dose mapping calculations inside the voxelized antropomorphic head phantom. The program provides the voxelization by first processing the CT slices; the process follows a twodimensional pixel and material identification algorithm on each slice and three-dimensional interpolation in order to describe the phantom geometry via small cubic cells, resulting in an MCNP input deck format output. Dose rates are calculated by using the MCNP5 tool FMESH, superimposed mesh tally, which gives the track length estimation of the particle flux in units of particles/cm2. Furthermore, the particle flux is converted into dose by using the conversion coefficients extracted from the NIST Physical Reference Data. The voxelization using a three-dimensional interpolation technique in combination with the use of the FMESH tool of the MCNP Monte Carlo code offers an optimal simulation which results in 3D dose mapping calculations inside anthropomorphic phantoms. This tool is very useful in radiation treatment assessments, in which voxelized phantoms are widely utilized. & 2009 Elsevier Ltd. All rights reserved.
1. Introduction Monte Carlo techniques play an essential role in the field of medical physics, especially in radiotherapy treatment planning. The present work reports on a methodology that has been developed using MCNP, a Monte Carlo based computer code, in order to calculate the dose delivered by a photon beam generated in a Theratron 780s (MDS Nordion) 60Co radiotherapy unit to a voxelized human phantom constructed from pre-segmented CT slices. Previous works (Miro´ et al., 2006) validated the cobalt therapy unit model by irradiating a simple heterogeneous water cube-shaped phantom model and comparing simulated and experimental results. The reference phantom model utilized in this work is the Zubal et al. (1994) phantom a creation of Dr George Zubal from Yale University, which consists of a series of pre-segmented CT slices of a human torso and head. The first approach of this project is to load these CT slices into the Matlabs code, which reconstructs the voxelization describing the phantom geometry via small cubic cells. Each voxel defines the mixture of the different materials that
Corresponding author. Tel.: + 34 963879631; fax: + 34 963877639.
E-mail address: [email protected] (R. Miro´). 0969-8043/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.apradiso.2009.10.016
compose it. The output of this code follows the MCNP input deck format and is integrated in a full input file including the specification of the 60Co radiotherapy unit. The second goal of the project is to simulate the whole system. The radiotherapy unit is included in the system as a phase-space file by means of the surface source SWW/SRR card, and the particle flux obtained at each voxel is converted into absolute dose via conversion coefficients from the NIST database. This work aims to the application of the voxelizing techniques to our own pre-segmented CT slices and the computation of absorbed dose in real patients-based phantoms, utilizing Monte Carlo techniques, and approaching simulation CPU times that are reasonable for radiotherapy treatment plans.
2. Simulation approach The following steps were crucial for the successful development of the simulation: The description of the cobalt 60Co radiotherapy unit and validation of the MCNP model, the voxelization of the Zubal phantom, the simulation of the phantom irradiation with a photon beam by MCNP, and the calculation of the dose from the particle flux via photon cross sections and electron stopping-power.
ARTICLE IN PRESS 710
V. Abella et al. / Applied Radiation and Isotopes 68 (2010) 709–713
2.1. Description of the radiotherapy unit In a previous work (Miro´ et al., 2006) the validation of a Theratron 780s (MDS Nordion) 60Co radiotherapy unit using the MCNP model was performed, by calculating the energy deposition in a homogeneous as well as in a heterogeneous phantom, consisting of a water tank with a polystyrene heterogeneity inside, and comparing the obtained simulation results with experimental data. Cobalt units use a 60Co radioactive source which is placed in the treatment head. To deliver dose to patients the radiation beam from the source is collimated by jaws. One of the reasons of using the Theratron 780s is the fact that the 60Cobalt gamma spectrum is well known and therefore easier to describe for modeling purposes. The model comprises a Theraton Radiotherapy irradiator positioned to give a beam focused on the axis origin of the phantom. The cobalt unit has a leaf collimator which provides rectangular fields of different sizes, in this case a field of 10 cm 10 cm is used. The 60Co source is made up of small pellets placed in a cylindrical capsule of stainless steel and surrounded by different materials, such as tungsten and brass. The cylindrical source capsule and its housing geometry as well as the cobalt activity with its 6 different gamma-ray lines spectrum have been modelled. A volumetric type of source has been specified, modeling a realistic cylindrical 60Co source capsule which radiates uniformly into a 4p solid angle. 2.2. Voxelization and the Zubal phantom The first approach to the voxelization of a human phantom consisted of an attempt to build up such phantom from series of pre-segmented CT slices, so that a simple algorithm created with Matlab would divide the slices into squares, evaluate the materials contained in these squares and make a 3-dimensional interpolation in order to define a mixture of the material contained in each voxel. The process of evaluation of the materials contained in the squares dividing each CT slice is carried out by pixel intensity identification and by labeling each square with the proportion of different pixel intensities that it contains. Subsequent to the 2dimensional calculation, a 3-dimensional interpolation is performed in order to build up the 3-dimensional voxelized phantom, which consists of a description of the mixture of pixel intensities in each voxel that composes it. Afterwards, a body tissue composition material and a density is assigned to each pixel intensity, resulting in voxels that depict a mixture of materials, and the whole phantom is described in a format compatible with the usual MCNP input deck format.
In order to validate this methodology, the pre-segmented Zubal et al. (1994) phantom was utilized. This phantom consists of 243 pre-segmented slices, 42 of which correspond to the ones defining the head. The Zubal phantom is a reconstruction from X-ray CT images of 512 512 pixels with a resolution of 1 mm in the XY plane and 0.5 cm (from the neck to the crown of the head) and 1 cm (for the rest of the body) in the Z-axis. The final phantom is interpolated to create a 128 128 243 byte volume with isotropic dimensions of 4 mm, and the portion used for this paper is of 128 128 42 byte that corresponds to the head of the phantom. Each pixel of the volume contains an index number designating it as belonging to a given organ or internal structure; 31 of such index numbers are assigned to the head. The process of obtaining the MCNP input model starts by inputting the 42 slices corresponding to the head of the Zubal phantom in the Matlab program. The program processes the slices throughout several subroutines, assigning an index number to each voxel of the desired dimensions (this case is 4 mm). Afterwards the program makes a material by index number identification and passes the information to the last routine which transforms it into a MCNP input. The 31 index numbers assigned to the head result in 13 different body tissue composition materials and densities, according to the specifications of the ICRU REPORT 46 (1992) and ICRU REPORT 44 (1989). The densities of each material are also stored in the program and related with the index numbers (Table 1). The resulting voxelized head model turns out to have too many voxels with an air material assigned, and therefore a reduction of the sides, up and down limiting voxels has been applied in the Matlab algorithm, in order to reduce the computing time for the MCNP simulation. The resulting MCNP model is a 177,828 voxels or cells (excluding the universe) input. 2.3. MCNP calculations The resulting input from the Matlab program was used for the MCNP simulation of the phantom irradiation. The 60Co radiotherapy unit model has been validated in previous works, as mentioned in Section 2.1. These results were obtained by making use of the SSW/SSR card, which allows writing a surface source file that corresponds to the starting surface for the particles to be simulated, in our case together with the human phantom. The number of particles registered at the surface is 940,700. This whole setting makes a surface source distance of 80 cm, which is the recommended for a 60Co therapy unit. The radiation transport in the phantom is performed following individual photon and electron histories through the whole geometry. A detailed photon physics treatment, including
Table 1 Recommended material compositions by ICRU REPORT 46 (1992) and ICRU REPORT 44 (1989). Tissue
Air Skin Brain Skull Vertebral column Skeletal muscle Adipose tissue Blood Bone marrow Cartilage Mandible Eye lens
r
Elemental composition (percentage by mass) H
C
N
O
Others
10.0 10.7 5.0 6.3 10.2 11.4 10.2 11.5 9.6 4.6 9.6
20.4 14.5 21.2 26.1 14.3 59.8 11.0 64.4 9.9 19.9 19.5
4.2 2.2 4.0 3.9 3.4 0.7 3.3 0.7 2.2 4.1 5.7
64.5 71.2 43.5 43.6 71.0 27.8 74.5 23.1 74.4 43.5 64.6
0.2Na,0.1P,0.2S,0.3Cl,0.1K 0.2Na,0.4P,0.2S,0.3Cl,0.3K 0.1Na,0.2Mg,8.1P,0.3S,17.6Ca 0.1Na,0.1Mg,6.1P,0.3S,0.1Cl, 0.1K,13.3Ca 0.1Na,0.2P,0.3S,0.1Cl,0.4K 0.1Na,0.1S,0.1Cl 0.1Na,0.1P,0.2S,0.3Cl,0.2K,0.1Fe 0.1Na,0.1S,0.1Cl 0.5Na,2.2P,0.9S,0.3Cl 0.1Na,0.2Mg,8.6P,0.3S,18.7Ca 0.1Na,0.1P,0.3S,0.1Cl
1090 1040 1610 1420 1050 950 1060 980 1100 1680 1070
ARTICLE IN PRESS V. Abella et al. / Applied Radiation and Isotopes 68 (2010) 709–713
photoelectric effect with fluorescence production, incoherent and coherent scattering and pair production, has been considered in the energy range between 0.001 and 2.6 MeV. In this case, a ‘photon and electron’ model (MODE PE card of MCNP) has been used, in which all photon collisions except coherent scattering can create electrons which are banked for later transport. Tally cards are used in MCNP to specify what type of information the user wants to obtain from the Monte Carlo calculation. In this work, an FMESH Tally is necessary in order to define a mesh tally superimposed to the voxel geometry, and to obtain information from each voxel, so that dose maps are easily reproduced. By default, the mesh tally calculates the track length estimate of the particle flux, averaged over a mesh cell (which corresponds to a voxel in our problem), in units of particles/cm2. The card has been organized to give the particle flux in each voxel divided in 25 energy bins from 0 to 1.6 MeV, using the EMESH
711
option, so that two separate MESHTAL files are created, one for photons and another for electrons, containing the results for each Z bin (voxels in the XY plane) and for each energy bin and the associated relative uncertainty. In order to accelerate the calculations, the MCNP code has been parallelized in as SGI Altix 3700 platform, using the MPI parallel protocol, in this case with 6 processors. This way, the problem is divided in 6 problems reducing the computerized time in the same proportion. Furthermore, the MCNP code has been adapted to the problem geometry by increasing the maximum number of cells allowed in the input geometry description to 200,000 with the Intel Fortran Compiler 9.0, on the Linux parallel computing machine (X-5 Monte Carlo Team, 2003). The final total computing time for the tallied photons is of 703.75 min (11.72 h), and for the tallied electrons is 713.66 min (11.89 h).
Fig. 1. Sketch of the head phantom provided by SABRINA. Full head and inside cut.
x 10-15 4 10
3.5 3
Y axis (cm)
5
2.5 0
2
-5
1.5 1
-10
0.5 -10 -8
-6
-4
-2
0
2
4
6
8
10
X axis (cm) Fig. 2. Relative dose rate isodose curves for an XY voxel plane situated at the center of the phantom, in Grays/s.
ARTICLE IN PRESS 712
V. Abella et al. / Applied Radiation and Isotopes 68 (2010) 709–713
2.4. Fluence to dose conversion The results of the FMESH tally are given in units of particles/ cm2 and, due to the complexity of the different material mixtures for the calculation of the photon cross sections and electron stopping power coefficients, the conversion to dose was managed separately using a Matlab program. According to this, both photon cross sections and electron stopping power were obtained from the NIST physical reference data web page http://physics.nist.gov/ PhysRefData/. This information is stored in matrixes to be used by
x 10-15 4 10
3.5
5 Y axis (cm)
the dose calculation algorithm. The other valuable information to be stored is the group of matrixes defining the particle flux at each plane and energy bin, taken from the MESHTAL file. The Matlab program starts its computing by first locating the voxel plane (z bin) at which the calculation is desired, and consequently identifying the different materials in each voxel. The material of the voxel is associated with a cross section coefficient for the photons and with the stopping power for the electrons on one hand, and with the corresponding particle fluence results for photons and electrons, respectively, on the other. A matrix entries multiplication is carried out following the formula for fluence to relative dose rate conversion. Finally, both dose rates for photon and electron are added and converted from MeV/g part s to Grays/part s by a factor of 1.69 10 9. The results are presented in the form of relative dose rate profile curves and relative depth dose curves, for the desired plane.
0
3
3. Results
2.5
The final results obtained by the whole process are presented in two manners. First, relative dose rate graph for the central XY plane (from 0.4 to 0 cm) is presented, in the form of isodose curves and an additional graph presenting the relative dose rate in each voxel. Finally, a relative depth dose graph is shown. The results, in Grays/s per particle emitted, are displayed in Figs. 1–4. An average relative error of 0.0623 for the photons is estimated at the central Z bin with X= 26. For the electrons, the average relative errors for the same location are estimated to be 0.123.
2 1.5
-5
1 -10
0.5 -10
-5
0
5
10
X axis (cm) Fig. 3. Relative dose rate in each voxel of the XY voxel plane situated at the center of the phantom, in Gray/s.
4. Conclusions This work provides a very useful procedure for the ongoing optimization of the radiotherapy treatment planning. The results
Fig. 4. Relative depth dose for the central XY plane and a central line X= 26, in %.
ARTICLE IN PRESS V. Abella et al. / Applied Radiation and Isotopes 68 (2010) 709–713
prove that a computerized voxelized human phantom is a valid tool for the calculations of absorbed dose inirradiated patients. Depth dose curves show that photon beams can be focused on patient tumors with a relative error of around 5%, for a 4 mm voxelized phantom. Bigger voxel sizes would lead to smaller relative errors in the calculations, although geometry description would not be so accurate. The future development of this methodology will consist in the voxelization of non-pre-segmented CT images, and coupling such program with the MCNP simulations and with the fluence-to-dose conversion algorithm. This work is a step forward in the increasing necessity for a user friendly radiotherapy treatment planning program based on Monte Carlo and for specific patientbased anthropomorphic phantoms in the medical physics field.
713
References International commission on radiation units and measurements, photon, electron, proton and neutron interaction data for body tissues, ICRU REPORT 46, 1992. International commission on radiation units and measurements, tissue substitutes in radiation dosimetry and measurement, ICRU REPORT 44, 1989. Miro´, R., Juste, B., Gallardo, S., Santos, A., Verdu´, G., 2006. Cobalt therapy dosimetric calculations over a voxelized heterogeneous phantom: validation of different Monte Carlo models and methodologies against experimental data. IEEE Transactions of Nuclear Science 53 (6), 3808–3817. National Institute of Standards and Technology, Physical reference data, available: /http://physics.nist.gov/PhysRefData/S. X-5 Monte Carlo Team, MCNP—a general Monte Carlo N-particle transport code version 5, LA_UR03-1987, Los Alamos National Labotratory, April 2003. Zubal, I.G., Harrell, C.R., Smith, E.O., Rattner, Z., Gindi, G., Hoffer, P.B., 1994. Computerized three-dimensional segmented human anatomy. Medical Physics 21 (2), 299–302.
Contents lists available at ScienceDirect
Applied Radiation and Isotopes journal homepage: www.elsevier.com/locate/apradiso
3D dose distribution calculation in a voxelized human phantom by means of Monte Carlo method V. Abella, R. Miro´ , B. Juste, G. Verdu´ ´cnica de Vale ncia, Camı´ de Vera, s/n. 46022 Vale ncia, Spain Departamento de Ingenierı´a Quı´mica y Nuclear, Universidad Polite
a r t i c l e in fo
Keywords: Biomedical applications of radiation Zubal phantom MCNP5 Voxelization
abstract The aim of this work is to provide the reconstruction of a real human voxelized phantom by means of a MatLabs program and the simulation of the irradiation of such phantom with the photon beam generated in a Theratron 780s (MDS Nordion) 60Co radiotherapy unit, by using the Monte Carlo transport code MCNP (Monte Carlo N-Particle), version 5. The project results in 3D dose mapping calculations inside the voxelized antropomorphic head phantom. The program provides the voxelization by first processing the CT slices; the process follows a twodimensional pixel and material identification algorithm on each slice and three-dimensional interpolation in order to describe the phantom geometry via small cubic cells, resulting in an MCNP input deck format output. Dose rates are calculated by using the MCNP5 tool FMESH, superimposed mesh tally, which gives the track length estimation of the particle flux in units of particles/cm2. Furthermore, the particle flux is converted into dose by using the conversion coefficients extracted from the NIST Physical Reference Data. The voxelization using a three-dimensional interpolation technique in combination with the use of the FMESH tool of the MCNP Monte Carlo code offers an optimal simulation which results in 3D dose mapping calculations inside anthropomorphic phantoms. This tool is very useful in radiation treatment assessments, in which voxelized phantoms are widely utilized. & 2009 Elsevier Ltd. All rights reserved.
1. Introduction Monte Carlo techniques play an essential role in the field of medical physics, especially in radiotherapy treatment planning. The present work reports on a methodology that has been developed using MCNP, a Monte Carlo based computer code, in order to calculate the dose delivered by a photon beam generated in a Theratron 780s (MDS Nordion) 60Co radiotherapy unit to a voxelized human phantom constructed from pre-segmented CT slices. Previous works (Miro´ et al., 2006) validated the cobalt therapy unit model by irradiating a simple heterogeneous water cube-shaped phantom model and comparing simulated and experimental results. The reference phantom model utilized in this work is the Zubal et al. (1994) phantom a creation of Dr George Zubal from Yale University, which consists of a series of pre-segmented CT slices of a human torso and head. The first approach of this project is to load these CT slices into the Matlabs code, which reconstructs the voxelization describing the phantom geometry via small cubic cells. Each voxel defines the mixture of the different materials that
Corresponding author. Tel.: + 34 963879631; fax: + 34 963877639.
E-mail address: [email protected] (R. Miro´). 0969-8043/$ - see front matter & 2009 Elsevier Ltd. All rights reserved. doi:10.1016/j.apradiso.2009.10.016
compose it. The output of this code follows the MCNP input deck format and is integrated in a full input file including the specification of the 60Co radiotherapy unit. The second goal of the project is to simulate the whole system. The radiotherapy unit is included in the system as a phase-space file by means of the surface source SWW/SRR card, and the particle flux obtained at each voxel is converted into absolute dose via conversion coefficients from the NIST database. This work aims to the application of the voxelizing techniques to our own pre-segmented CT slices and the computation of absorbed dose in real patients-based phantoms, utilizing Monte Carlo techniques, and approaching simulation CPU times that are reasonable for radiotherapy treatment plans.
2. Simulation approach The following steps were crucial for the successful development of the simulation: The description of the cobalt 60Co radiotherapy unit and validation of the MCNP model, the voxelization of the Zubal phantom, the simulation of the phantom irradiation with a photon beam by MCNP, and the calculation of the dose from the particle flux via photon cross sections and electron stopping-power.
ARTICLE IN PRESS 710
V. Abella et al. / Applied Radiation and Isotopes 68 (2010) 709–713
2.1. Description of the radiotherapy unit In a previous work (Miro´ et al., 2006) the validation of a Theratron 780s (MDS Nordion) 60Co radiotherapy unit using the MCNP model was performed, by calculating the energy deposition in a homogeneous as well as in a heterogeneous phantom, consisting of a water tank with a polystyrene heterogeneity inside, and comparing the obtained simulation results with experimental data. Cobalt units use a 60Co radioactive source which is placed in the treatment head. To deliver dose to patients the radiation beam from the source is collimated by jaws. One of the reasons of using the Theratron 780s is the fact that the 60Cobalt gamma spectrum is well known and therefore easier to describe for modeling purposes. The model comprises a Theraton Radiotherapy irradiator positioned to give a beam focused on the axis origin of the phantom. The cobalt unit has a leaf collimator which provides rectangular fields of different sizes, in this case a field of 10 cm 10 cm is used. The 60Co source is made up of small pellets placed in a cylindrical capsule of stainless steel and surrounded by different materials, such as tungsten and brass. The cylindrical source capsule and its housing geometry as well as the cobalt activity with its 6 different gamma-ray lines spectrum have been modelled. A volumetric type of source has been specified, modeling a realistic cylindrical 60Co source capsule which radiates uniformly into a 4p solid angle. 2.2. Voxelization and the Zubal phantom The first approach to the voxelization of a human phantom consisted of an attempt to build up such phantom from series of pre-segmented CT slices, so that a simple algorithm created with Matlab would divide the slices into squares, evaluate the materials contained in these squares and make a 3-dimensional interpolation in order to define a mixture of the material contained in each voxel. The process of evaluation of the materials contained in the squares dividing each CT slice is carried out by pixel intensity identification and by labeling each square with the proportion of different pixel intensities that it contains. Subsequent to the 2dimensional calculation, a 3-dimensional interpolation is performed in order to build up the 3-dimensional voxelized phantom, which consists of a description of the mixture of pixel intensities in each voxel that composes it. Afterwards, a body tissue composition material and a density is assigned to each pixel intensity, resulting in voxels that depict a mixture of materials, and the whole phantom is described in a format compatible with the usual MCNP input deck format.
In order to validate this methodology, the pre-segmented Zubal et al. (1994) phantom was utilized. This phantom consists of 243 pre-segmented slices, 42 of which correspond to the ones defining the head. The Zubal phantom is a reconstruction from X-ray CT images of 512 512 pixels with a resolution of 1 mm in the XY plane and 0.5 cm (from the neck to the crown of the head) and 1 cm (for the rest of the body) in the Z-axis. The final phantom is interpolated to create a 128 128 243 byte volume with isotropic dimensions of 4 mm, and the portion used for this paper is of 128 128 42 byte that corresponds to the head of the phantom. Each pixel of the volume contains an index number designating it as belonging to a given organ or internal structure; 31 of such index numbers are assigned to the head. The process of obtaining the MCNP input model starts by inputting the 42 slices corresponding to the head of the Zubal phantom in the Matlab program. The program processes the slices throughout several subroutines, assigning an index number to each voxel of the desired dimensions (this case is 4 mm). Afterwards the program makes a material by index number identification and passes the information to the last routine which transforms it into a MCNP input. The 31 index numbers assigned to the head result in 13 different body tissue composition materials and densities, according to the specifications of the ICRU REPORT 46 (1992) and ICRU REPORT 44 (1989). The densities of each material are also stored in the program and related with the index numbers (Table 1). The resulting voxelized head model turns out to have too many voxels with an air material assigned, and therefore a reduction of the sides, up and down limiting voxels has been applied in the Matlab algorithm, in order to reduce the computing time for the MCNP simulation. The resulting MCNP model is a 177,828 voxels or cells (excluding the universe) input. 2.3. MCNP calculations The resulting input from the Matlab program was used for the MCNP simulation of the phantom irradiation. The 60Co radiotherapy unit model has been validated in previous works, as mentioned in Section 2.1. These results were obtained by making use of the SSW/SSR card, which allows writing a surface source file that corresponds to the starting surface for the particles to be simulated, in our case together with the human phantom. The number of particles registered at the surface is 940,700. This whole setting makes a surface source distance of 80 cm, which is the recommended for a 60Co therapy unit. The radiation transport in the phantom is performed following individual photon and electron histories through the whole geometry. A detailed photon physics treatment, including
Table 1 Recommended material compositions by ICRU REPORT 46 (1992) and ICRU REPORT 44 (1989). Tissue
Air Skin Brain Skull Vertebral column Skeletal muscle Adipose tissue Blood Bone marrow Cartilage Mandible Eye lens
r
Elemental composition (percentage by mass) H
C
N
O
Others
10.0 10.7 5.0 6.3 10.2 11.4 10.2 11.5 9.6 4.6 9.6
20.4 14.5 21.2 26.1 14.3 59.8 11.0 64.4 9.9 19.9 19.5
4.2 2.2 4.0 3.9 3.4 0.7 3.3 0.7 2.2 4.1 5.7
64.5 71.2 43.5 43.6 71.0 27.8 74.5 23.1 74.4 43.5 64.6
0.2Na,0.1P,0.2S,0.3Cl,0.1K 0.2Na,0.4P,0.2S,0.3Cl,0.3K 0.1Na,0.2Mg,8.1P,0.3S,17.6Ca 0.1Na,0.1Mg,6.1P,0.3S,0.1Cl, 0.1K,13.3Ca 0.1Na,0.2P,0.3S,0.1Cl,0.4K 0.1Na,0.1S,0.1Cl 0.1Na,0.1P,0.2S,0.3Cl,0.2K,0.1Fe 0.1Na,0.1S,0.1Cl 0.5Na,2.2P,0.9S,0.3Cl 0.1Na,0.2Mg,8.6P,0.3S,18.7Ca 0.1Na,0.1P,0.3S,0.1Cl
1090 1040 1610 1420 1050 950 1060 980 1100 1680 1070
ARTICLE IN PRESS V. Abella et al. / Applied Radiation and Isotopes 68 (2010) 709–713
photoelectric effect with fluorescence production, incoherent and coherent scattering and pair production, has been considered in the energy range between 0.001 and 2.6 MeV. In this case, a ‘photon and electron’ model (MODE PE card of MCNP) has been used, in which all photon collisions except coherent scattering can create electrons which are banked for later transport. Tally cards are used in MCNP to specify what type of information the user wants to obtain from the Monte Carlo calculation. In this work, an FMESH Tally is necessary in order to define a mesh tally superimposed to the voxel geometry, and to obtain information from each voxel, so that dose maps are easily reproduced. By default, the mesh tally calculates the track length estimate of the particle flux, averaged over a mesh cell (which corresponds to a voxel in our problem), in units of particles/cm2. The card has been organized to give the particle flux in each voxel divided in 25 energy bins from 0 to 1.6 MeV, using the EMESH
711
option, so that two separate MESHTAL files are created, one for photons and another for electrons, containing the results for each Z bin (voxels in the XY plane) and for each energy bin and the associated relative uncertainty. In order to accelerate the calculations, the MCNP code has been parallelized in as SGI Altix 3700 platform, using the MPI parallel protocol, in this case with 6 processors. This way, the problem is divided in 6 problems reducing the computerized time in the same proportion. Furthermore, the MCNP code has been adapted to the problem geometry by increasing the maximum number of cells allowed in the input geometry description to 200,000 with the Intel Fortran Compiler 9.0, on the Linux parallel computing machine (X-5 Monte Carlo Team, 2003). The final total computing time for the tallied photons is of 703.75 min (11.72 h), and for the tallied electrons is 713.66 min (11.89 h).
Fig. 1. Sketch of the head phantom provided by SABRINA. Full head and inside cut.
x 10-15 4 10
3.5 3
Y axis (cm)
5
2.5 0
2
-5
1.5 1
-10
0.5 -10 -8
-6
-4
-2
0
2
4
6
8
10
X axis (cm) Fig. 2. Relative dose rate isodose curves for an XY voxel plane situated at the center of the phantom, in Grays/s.
ARTICLE IN PRESS 712
V. Abella et al. / Applied Radiation and Isotopes 68 (2010) 709–713
2.4. Fluence to dose conversion The results of the FMESH tally are given in units of particles/ cm2 and, due to the complexity of the different material mixtures for the calculation of the photon cross sections and electron stopping power coefficients, the conversion to dose was managed separately using a Matlab program. According to this, both photon cross sections and electron stopping power were obtained from the NIST physical reference data web page http://physics.nist.gov/ PhysRefData/. This information is stored in matrixes to be used by
x 10-15 4 10
3.5
5 Y axis (cm)
the dose calculation algorithm. The other valuable information to be stored is the group of matrixes defining the particle flux at each plane and energy bin, taken from the MESHTAL file. The Matlab program starts its computing by first locating the voxel plane (z bin) at which the calculation is desired, and consequently identifying the different materials in each voxel. The material of the voxel is associated with a cross section coefficient for the photons and with the stopping power for the electrons on one hand, and with the corresponding particle fluence results for photons and electrons, respectively, on the other. A matrix entries multiplication is carried out following the formula for fluence to relative dose rate conversion. Finally, both dose rates for photon and electron are added and converted from MeV/g part s to Grays/part s by a factor of 1.69 10 9. The results are presented in the form of relative dose rate profile curves and relative depth dose curves, for the desired plane.
0
3
3. Results
2.5
The final results obtained by the whole process are presented in two manners. First, relative dose rate graph for the central XY plane (from 0.4 to 0 cm) is presented, in the form of isodose curves and an additional graph presenting the relative dose rate in each voxel. Finally, a relative depth dose graph is shown. The results, in Grays/s per particle emitted, are displayed in Figs. 1–4. An average relative error of 0.0623 for the photons is estimated at the central Z bin with X= 26. For the electrons, the average relative errors for the same location are estimated to be 0.123.
2 1.5
-5
1 -10
0.5 -10
-5
0
5
10
X axis (cm) Fig. 3. Relative dose rate in each voxel of the XY voxel plane situated at the center of the phantom, in Gray/s.
4. Conclusions This work provides a very useful procedure for the ongoing optimization of the radiotherapy treatment planning. The results
Fig. 4. Relative depth dose for the central XY plane and a central line X= 26, in %.
ARTICLE IN PRESS V. Abella et al. / Applied Radiation and Isotopes 68 (2010) 709–713
prove that a computerized voxelized human phantom is a valid tool for the calculations of absorbed dose inirradiated patients. Depth dose curves show that photon beams can be focused on patient tumors with a relative error of around 5%, for a 4 mm voxelized phantom. Bigger voxel sizes would lead to smaller relative errors in the calculations, although geometry description would not be so accurate. The future development of this methodology will consist in the voxelization of non-pre-segmented CT images, and coupling such program with the MCNP simulations and with the fluence-to-dose conversion algorithm. This work is a step forward in the increasing necessity for a user friendly radiotherapy treatment planning program based on Monte Carlo and for specific patientbased anthropomorphic phantoms in the medical physics field.
713
References International commission on radiation units and measurements, photon, electron, proton and neutron interaction data for body tissues, ICRU REPORT 46, 1992. International commission on radiation units and measurements, tissue substitutes in radiation dosimetry and measurement, ICRU REPORT 44, 1989. Miro´, R., Juste, B., Gallardo, S., Santos, A., Verdu´, G., 2006. Cobalt therapy dosimetric calculations over a voxelized heterogeneous phantom: validation of different Monte Carlo models and methodologies against experimental data. IEEE Transactions of Nuclear Science 53 (6), 3808–3817. National Institute of Standards and Technology, Physical reference data, available: /http://physics.nist.gov/PhysRefData/S. X-5 Monte Carlo Team, MCNP—a general Monte Carlo N-particle transport code version 5, LA_UR03-1987, Los Alamos National Labotratory, April 2003. Zubal, I.G., Harrell, C.R., Smith, E.O., Rattner, Z., Gindi, G., Hoffer, P.B., 1994. Computerized three-dimensional segmented human anatomy. Medical Physics 21 (2), 299–302.