research papers\(\def\hfill{\hskip 5em}\def\hfil{\hskip 3em}\def\eqno#1{\hfil {#1}}\)

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775

An efficient numerical tool for dose deposition prediction applied to synchrotron medical imaging and radiation therapy

CROSSMARK_Color_square_no_text.svg

aDepartment of Physics, Ludwig Maximilians University, Am Coulombwall 1, Munich, Germany, bDepartment of Clinical Radiology, Ludwig Maximilians University, Munich, Germany, cUniversité de Lyon, CREATIS; CNRS UMR5220; Inserm U1044; INSA-Lyon; Université Lyon 1; Centre Léon Bérard, France, and dEuropean Synchrotron Radiation Facility, Grenoble, France
*Correspondence e-mail: alberto.mittone@physik.uni-muenchen.de

(Received 3 May 2013; accepted 21 June 2013; online 5 July 2013)

Medical imaging and radiation therapy are widely used synchrotron-based techniques which have one thing in common: a significant dose delivery to typically biological samples. Among the ways to provide the experimenters with image guidance techniques indicating optimization strategies, Monte Carlo simulation has become the gold standard for accurately predicting radiation dose levels under specific irradiation conditions. A highly important hampering factor of this method is, however, its slow statistical convergence. A track length estimator (TLE) module has been coded and implemented for the first time in the open-source Monte Carlo code GATE/Geant4. Results obtained with the module and the procedures used to validate them are presented. A database of energy-absorption coefficients was also generated, which is used by the TLE calculations and is now also included in GATE/Geant4. The validation was carried out by comparing the TLE-simulated doses with experimental data in a synchrotron radiation computed tomography experiment. The TLE technique shows good agreement versus both experimental measurements and the results of a classical Monte Carlo simulation. Compared with the latter, it is possible to reach a pre-defined statistical uncertainty in about two to three orders of magnitude less time for complex geometries without loss of accuracy.

1. Introduction

Knowledge of the dose delivered during medical imaging and radiation therapy with a synchrotron source represents a considerable issue. In order to preserve the integrity of biological samples (but also other samples), their exposure should be carefully guided following dosimetric criteria which all lead to the basic paradigm that the absorbed dose should be as low as reasonably achievable. In this respect, therapeutic and imaging doses, when combined, cannot be regarded from the same viewpoint because the imaging dose adds to an already high level of therapeutic radiation.

The Monte Carlo (MC) method is often used to compute the distribution of the dose deposited in both imaging and radiotherapy applications. It provides accurate results but still requires long computation times, hindering its use in clinical routine applications. Many MC codes are available to simulate radiation transport (see Table 1[link]).

Table 1
Overview of the main MC codes to simulate radiation transport

MC code Reference
Geant4 Agostinelli et al. (2003[Agostinelli, S. et al. (2003). Nucl. Instrum. Methods Phys. Res. A, 506, 250-303.]), Allison et al. (2006[Allison J. et al. (2006). IEEE Trans. Nucl. Sci. 53, 270-278.])
MCNPX Pelowitz (2005[Pelowitz, D. B. (2005). Editor. MCNPX User's Manual, p. 473. Los Alamos National Laboratory, New Mexico, USA.])
FLUKA Battistoni et al. (2007[Battistoni, G., Cerutti, F., Fassò, A., Ferrari, A., Muraro, S., Ranft, J., Roesler, S. & Sala, P. R. (2007). AIP Conf. Proc. 896, 31-49.]), Ferrari et al. (2005[Ferrari, A., Sala, P. R., Fasso, A. & Ranft, J. (2005). FLUKA Technical Report, CERN-2005-10, INFN/TC 05/11, SLAC-R-773.])
EGSnrc Fragoso et al. (2009[Fragoso, M., Kawrakow, I., Faddegon, B. A., Solberg, T. D. & Chetty, I. J. (2009). Med. Phys. 36, 5451-5466.])
PENELOPE Salvat & Fernández-Varea (2009[Salvat, F. & Fernández-Varea, J. M. (2009). Metrologia, 46, S112-S138.]), Salvat et al. (2011[Salvat, F., Fernández-Varea, J. M. & Sempau, J. (2011). PENELOPE. A Code System for Monte Carlo Simulation of Electron and Photon Transport. Issy-les-Moulineaux: OECD Nuclear Energy Agency.])

They are commonly used for research purposes, either in simulation mode or in combination with various acceleration techniques, such as condensed history methods or variance reduction techniques (VRT), e.g. splitting and Russian roulette, importance sampling, forced detection, etc. (Jenkins et al., 1988[Jenkins, T. M., Nelson, W. R. & Rindi, A. (1988). Monte Carlo Transport of Electrons and Photons. New York: Plenum Press.]). These acceleration techniques provide substantial speed-up but have to be used very cautiously to stay within their limits of validity. Among the VRTs, the track length estimator (TLE) method is known as a very efficient tallying method, well suited for kerma calculations at any given point during photon irradiations (Williamson, 1987[Williamson, J. F. (1987). Med. Phys. 14, 567-576.]).

We investigate the application of TLE to speed up dose calculations for low-energy (from 1 keV to a few hundred keV) photon irradiations, which include radiology, brachytherapy, synchrotron radiation therapy and imaging applications. This article is structured as follows: in §2[link], the TLE method and the related physical issues are first introduced, its implementation in the next release of the open source code GATE (Geant4 Application for Emission Tomography, version 9.5) (Jan et al., 2004[Jan, S. et al. (2004). Phys. Med. Biol. 49, 4543-4561.], 2011[Jan, S. et al. (2011). Phys. Med. Biol. 56, 881-901.]) is described, then the test cases for validation are presented as well as the dose measurement protocol and associated set-up; §3[link] presents the results of the TLE implementation against analogous MC simulations and then their validation against experimental data obtained with a high-resolution computed tomography (CT) test case using monochromatic X-rays (issued from a synchrotron radiation source). Finally, the TLE method is shown to properly cope with very complex geometries with no need for heavy computing resources.

2. Materials and methods

2.1. The TLE method

The use of a TLE technique for calculating particle fluences, kerma and absorbed doses has long been known (Carlsson, 1985[Carlsson, G. A. (1985). The Dosimetry of Ionizing Radiation, pp. 1-75. Orlando: Academic Press.]; Williamson, 1987[Williamson, J. F. (1987). Med. Phys. 14, 567-576.]). It is implemented, for example, in the MCNPX MC code (Demarco et al., 2002[Demarco, J. J., Wallace, R. E. & Boedeker, K. (2002). Phys. Med. Biol. 47, 1321-1332.]; Smans et al., 2010[Smans, K., Zoetelief, J., Verbrugge, B., Haeck, W., Struelens, L., Vanhavere, F. & Bosmans, H. (2010). Med. Phys. 37, 2082-2091.]) and in some dedicated tools, notably for external radiotherapy (van der Zee et al., 2005[Zee, W. van der, Hogenbirk, A. & van der Marck, S. C. (2005). Phys. Med. Biol. 50, 625-641.]) and brachytherapy applications (Chibani & Williamson, 2005[Chibani, O. & Williamson, J. F. (2005). Med. Phys. 32, 3688-3698.]; Taylor et al., 2007[Taylor, R. E., Yegin, G. & Rogers, D. W. (2007). Med. Phys. 34, 445-457.]). Considering a monoenergetic photon beam with energy E, the absorbed dose in charged particle equilibrium (CPE) is given by (Carlsson, 1985[Carlsson, G. A. (1985). The Dosimetry of Ionizing Radiation, pp. 1-75. Orlando: Academic Press.]; Berger et al., 2010[Berger, M. J., Hubbell, J. H., Seltzer, S. M., Chang, J., Coursey, J. S., Sukumar, R., Zucker, D. S. & Olsen, K. (2010). XCOM: Photon Cross Section Database (version 1.5), https://physics.nist.gov/xcom .])

[D=\Phi{E}\,{{{\mu_{\rm{en}}}}/{\rho}},\eqno(1)]

where Φ is the particle fluence and μen/ρ is the mass energy-absorption coefficient. Fluence can be interpreted as the track length density of the particles at a point r in space (Carlsson, 1985[Carlsson, G. A. (1985). The Dosimetry of Ionizing Radiation, pp. 1-75. Orlando: Academic Press.]), i.e.

[\Phi({\bf{r}})={{{\rm{d}}L({\bf{r}})}/{{\rm{d}}V}}.\eqno(2)]

For a photon traversing a voxel of volume V, an estimate of the corresponding fluence is therefore given by (Williamson, 1987[Williamson, J. F. (1987). Med. Phys. 14, 567-576.]; Demarco et al., 2002[Demarco, J. J., Wallace, R. E. & Boedeker, K. (2002). Phys. Med. Biol. 47, 1321-1332.])

[\Phi={L/V},\eqno(3)]

where L is the track length, i.e. the straight-line distance travelled in the voxel between successive collisions. TLE can thus be used to determine the expected value of the dose deposited along every photon trajectory given by (1)[link]. With the TLE scoring method, a photon deposits its energy in all voxels it encounters between successive interaction points, instead of doing so only at interaction points as is the case in an equivalent MC simulation. The major advantage of the use of the TLE is a drastic acceleration of the convergence of the simulation.

2.2. Local energy deposition by secondary electrons

According to the TLE method, secondary electrons are not tracked and their energy is deposited locally (within a single voxel). This approximation is satisfactory when the CPE condition is fulfilled and the electron range remains smaller than either the voxel size or the required spatial accuracy. Comparing the electron range in the continuous-slowing-down approximation (CSDA range) with the voxel size is a conservative criterion, since the projected range is always smaller than the CSDA range (detour factor below 1) (Berger et al., 2005[Berger, M. J. et al. (2005). National Institute of Standards and Technology, Gaithersburg, MD, USA.]). For soft tissues (mass density 1.06 g cm−3), the CSDA range stays below 1 mm up to about 300 keV (see ICRU, 1984[ICRU (1984). ICRU Report 37. International Commission on Radiation Units and Measurements, Bethesda, MD, USA.]).

2.3. Energy-absorption tables

In the energy domain of interest in this study (keV X-rays), the linear energy-absorption (μen) and energy-transfer (μtr) coefficients can be considered to be the same quantity (Attix, 2004[Attix, F. H. (2004). Introduction to Radiological Physics and Radiation Dosimetry. Weinheim: Wiley-VCH.]; Berger et al., 2010[Berger, M. J., Hubbell, J. H., Seltzer, S. M., Chang, J., Coursey, J. S., Sukumar, R., Zucker, D. S. & Olsen, K. (2010). XCOM: Photon Cross Section Database (version 1.5), https://physics.nist.gov/xcom .]; Freud et al., 2008[Freud, N., Letang, J. M., Mary, C., Boudou, C., Ferrero, C., Elleaume, H., Bravin, A., Esteve, F. & Babot, D. (2008). IEEE Trans. Nucl. Sci. 55, 1008-1017.]),

[{\mu_{{\rm{en}}}}\simeq{\mu_{{\rm{tr}}}}.\eqno(4)]

For tissues composed of elements with atomic numbers Z ≤ 20, the relative difference between μtr and μen remains below 1% for energy values up to 1 MeV. The difference can reach a few percent at 1 MeV for elements with Z > 20 (Attix, 2004[Attix, F. H. (2004). Introduction to Radiological Physics and Radiation Dosimetry. Weinheim: Wiley-VCH.]). In the case of a mixture of elements, the mass energy-transfer coefficient μtr/ρ satisfies the additivity rule.

The values of μen/ρ may be obtained in several ways:

(i) Directly from the NIST database for elemental media and a limited number of materials of interest (Berger et al., 2010[Berger, M. J., Hubbell, J. H., Seltzer, S. M., Chang, J., Coursey, J. S., Sukumar, R., Zucker, D. S. & Olsen, K. (2010). XCOM: Photon Cross Section Database (version 1.5), https://physics.nist.gov/xcom .]). In the case of compounds or mixtures, the additivity rule can be used as an approximation.

(ii) Using the EPDL97 database and the approximation by equation (4)[link]. Details about the calculation of the coefficients can be found on the NIST webpage (Berger et al., 2010[Berger, M. J., Hubbell, J. H., Seltzer, S. M., Chang, J., Coursey, J. S., Sukumar, R., Zucker, D. S. & Olsen, K. (2010). XCOM: Photon Cross Section Database (version 1.5), https://physics.nist.gov/xcom .]).

(iii) Using MC simulations [method described by Freud et al. (2008[Freud, N., Letang, J. M., Mary, C., Boudou, C., Ferrero, C., Elleaume, H., Bravin, A., Esteve, F. & Babot, D. (2008). IEEE Trans. Nucl. Sci. 55, 1008-1017.])]. In Geant4, several low-energy electromagnetics physics lists based on Livermore or PENELOPE models (Ivanchenko et al., 2011[Ivanchenko, V. et al. (2011). Prog. Nucl. Sci. Technol. 2, 898-903.]) are available for reference. Both use EPDL97 and EADL to describe photon interactions and atomic relaxation, respectively. As regards Compton scattering, both Livermore and PENELOPE models predict the ionized shell and subsequent atomic relaxation. Note, however, that the fluorescence emission following Compton scattering is a very rare process and has a negligible influence on the μen value.

Fig. 1[link] presents a comparison between μtr/ρ calculated for cortical bone and lung tissue (ICRU, 1989[ICRU (1989). ICRU Report 44. International Commission on Radiation Units and Measurements, Bethesda, MD, USA.]) using the EPDL97 database (Cullen et al., 1997[Cullen, D. E., Hubbell, J. H. & Kissel, L. (1997). EPDL97: The Evaluated Photon Data Library, UCRL-50400, Vol. 6, Rev. 5, https://www-nds.iaea.org/epdl97/document/epdl97.pdf .]), see (ii), μtr/ρ derived from Geant4 (version 9.5) simulations (iii) with the Livermore model, and μen/ρ from NIST (Berger et al., 2010[Berger, M. J., Hubbell, J. H., Seltzer, S. M., Chang, J., Coursey, J. S., Sukumar, R., Zucker, D. S. & Olsen, K. (2010). XCOM: Photon Cross Section Database (version 1.5), https://physics.nist.gov/xcom .]), (i). No significant difference is detected between these three sets of data.

[Figure 1]
Figure 1
Comparison between energy absorption-to-attenuation coefficient ratios calculated for cortical bone and lung tissue (ICRU, 1989[ICRU (1989). ICRU Report 44. International Commission on Radiation Units and Measurements, Bethesda, MD, USA.]) using the EPDL97 database (Cullen et al., 1997[Cullen, D. E., Hubbell, J. H. & Kissel, L. (1997). EPDL97: The Evaluated Photon Data Library, UCRL-50400, Vol. 6, Rev. 5, https://www-nds.iaea.org/epdl97/document/epdl97.pdf .]), as derived from Geant4 (version 9.5) simulations with the Livermore model [according to the method described by Freud et al. (2008[Freud, N., Letang, J. M., Mary, C., Boudou, C., Ferrero, C., Elleaume, H., Bravin, A., Esteve, F. & Babot, D. (2008). IEEE Trans. Nucl. Sci. 55, 1008-1017.])], and as given by NIST (Berger et al., 2010[Berger, M. J., Hubbell, J. H., Seltzer, S. M., Chang, J., Coursey, J. S., Sukumar, R., Zucker, D. S. & Olsen, K. (2010). XCOM: Photon Cross Section Database (version 1.5), https://physics.nist.gov/xcom .]).

In the case of iodine (Fig. 2[link]), slight differences are found between the μtr/ρ calculated via EPDL97 and μen/ρ from NIST, up to about 4% at 1 MeV, corresponding to the average fraction of the kinetic energy of secondary charged particles escaping through radiative processes (Berger et al., 2010[Berger, M. J., Hubbell, J. H., Seltzer, S. M., Chang, J., Coursey, J. S., Sukumar, R., Zucker, D. S. & Olsen, K. (2010). XCOM: Photon Cross Section Database (version 1.5), https://physics.nist.gov/xcom .]). Moreover, no significant difference (below about 1%) appears between the μtr/ρ values calculated through EPDL97 and estimated by Geant4, whichever model (PENELOPE or Livermore) is used.

[Figure 2]
Figure 2
Comparison between energy absorption-to-attenuation coefficient ratios calculated for iodine using the EPDL97 database (Cullen et al., 1997[Cullen, D. E., Hubbell, J. H. & Kissel, L. (1997). EPDL97: The Evaluated Photon Data Library, UCRL-50400, Vol. 6, Rev. 5, https://www-nds.iaea.org/epdl97/document/epdl97.pdf .]), simulated using Geant4 (version 9.5) in conjunction with the PENELOPE and Livermore physics models, and as provided by NIST (Berger et al., 2010[Berger, M. J., Hubbell, J. H., Seltzer, S. M., Chang, J., Coursey, J. S., Sukumar, R., Zucker, D. S. & Olsen, K. (2010). XCOM: Photon Cross Section Database (version 1.5), https://physics.nist.gov/xcom .]). The curves corresponding to the Geant4, PENELOPE and Livermore models can barely be distinguished. The relative difference between EPDL97 and Geant4 stays below about 1%. When energy increases, the difference between μtr/ρ (either from EPDL97 or Geant4) and μen/ρ from NIST increases up to about 4% at 1 MeV. In the EPDL97 library the Doppler broadening effect is not considered. It has been shown (Ye et al., 2006[Ye, S. J., Ove, R. & Naqvi, S. A. (2006). Health Phys. 91, 361-366.]) that the usage of this approximation compared with a more advanced model (impulse approximation) does not produce significant differences in the results, at least for the considered energy range.

The database we generated from the EPDL97 library, included in the next release of GATE, incorporates the μen coefficients mentioned above for all elements with 1 ≤ Z ≤ 100 in the energy range [1 keV, 1 MeV] and has been used in all TLE calculations reported here.

2.4. Implementation of the TLE in GATE

During the initialization of a simulation, the pre-generated database, presented in §2.3[link], is loaded. The coefficients are tabulated for a set of 117–868 energies and for each Z value.

The discontinuities are treated by insertion of two identical energies with different μen values. When the energy of the photon is equivalent to the absorption edge energy the higher values of μen are considered for the simulation.

Materials are defined as mixtures of n elements with corresponding mass fractions wn. When an energy value used by the simulation does not appear in the database, the value of the corresponding μen is calculated on the fly by means of a logarithmic interpolation between the μen values referring to the two nearest energies. We implemented this method in GATE (version 6.2) and it will be available for use in the next release of GATE. In order to use it, it is sufficient to call TLEDoseActor in a way similar to the pre-existing DoseActor (Jan et al., 2011[Jan, S. et al. (2011). Phys. Med. Biol. 56, 881-901.]). The distribution maps of the energy deposited in the volume V of interest are created if the actor is attached to V. The volume must be a matrix of voxels named NestedParametrizedVolume. By definition, one element of the dose score matrix corresponds to one dosel. The implementation of a look-up table to store the tabulated values μen of the already defined materials is foreseen to accelerate the computations.

2.5. Benchmarking of the TLE implementation in GATE against MC simulation

The test case chosen for the comparison between TLE and equivalent classical MC simulations is a chest radiography consisting of a CT model (obtained during a real CT acquisition and provided in the set of examples distributed with GATE) of a thorax phantom (512 × 512 × 53 voxels of 0.602 mm × 0.602 mm × 3.000 mm; see Fig. 3[link]) irradiated by an X-ray tube, simulated as a point source with a 81 kVp polychromatic spectrum. The thorax phantom is the Dynamic Thorax Phantom from CIRS Inc. It is composed of several different materials (in terms of both density and linear attenuation) included between air and metal implant. The distance from the source to the phantom centre was set to 100 cm. Simulations were carried out using Geant4 (version 9.5) with the PENELOPE physics models. The voxel size of dose and energy deposited maps was set to 3.08 mm × 3.08 mm × 9.9375 mm.

[Figure 3]
Figure 3
CT model of the thorax phantom used in the chest radiography simulation test case. (a) Transverse slice; (b) sagittal slice.

2.6. Experimental validation of the TLE

Measurements were performed at the biomedical beamline (ID17) of the European Synchrotron Radiation Facility (ESRF, Grenoble, France), pursuing a twofold objective: (i) evaluating the photon flux, i.e. the photon fluence per unit time and per storage-ring electron current necessary to estimate the number of histories to be used in the simulations, and (ii) determining dose values in different experimental configurations to be compared with dosimetric results simulated by our TLE code.

2.6.1. Photon flux measurement protocol

In the simulations, the number of histories nh needed to compute the total deposited dose is given by nh = φItS, where S is the beam cross-sectional surface area, I is the storage ring electron current, t is the exposure time and φ is the photon flux (photons s−1 mm−2 mA−1).

The photon flux was derived from the dose absorbed in water dosimetry measurements, carried out using a ionization chamber connected to a PTW dosimeter (PTW 31010 calibrated by PTB Freiburg, Germany) at different X-ray energies in the range 33.7–72 keV. The active volume of the ionization chamber is 0.125 mm3. As the height of the synchrotron beam is smaller than the chamber size, measurements were performed by vertically scanning the ionization chamber with the photon beam (Prezado et al., 2011[Prezado, Y., Vautrin, M., Martínez-Rovira, I., Bravin, A., Estève, F., Elleaume, H., Berkvens, P. & Adam, J. F. (2011). Med. Phys. 38, 1709-1717.]). The chamber was positioned perpendicularly to the direction of the beam: in this case no correction concerning the directional dependence of the device had to be applied. A correction taking into account the energy behaviour of the chamber was applied. The correction factor curves are provided directly by the manufacturer of the chamber. The measurements were repeated twice and the mean value was taken as the reference estimate. The photon flux at the point of measurement was evaluated as

[{\varphi_{{\rm{PTW}}}}={{D\rho}\over{E{\mu_{{\rm{en}}}}tI}},\eqno(5)]

where D is the dose recorded by the ionization chamber. Note that equation (5)[link] derives from equation (1)[link] (keeping the same notations), and is valid in the case of a monochromatic X-ray beam with negligible scattering, which is a satisfactory approximation when the ionization chamber is placed in air.

A high-efficiency germanium detector (HEGD) was used to cross-check the results delivered by the PTW dosimeter. The HEGD consists of a high-purity germanium detector operating at liquid-nitrogen temperature. It is segmented into two rows of 10 mm height and 432 parallel strips with 0.35 mm pitch. The number of counts was determined using the detector-specific energy-dependent photon-counts relation, averaging over repeated acquisitions. From the mean value b of counts (expressed in bits), the photon flux at the point of measurement was calculated as

[{\varphi_{{\rm{HEGD}}}}={{kb}/{{S_{{\rm{pix}}}}\,tI}}, \eqno(6)]

where k is the bit–photon conversion factor (corresponding to the number of photons integrated by the HEGD per bit) and Spix is the irradiated pixel surface area (Coan et al., 2006[Coan, P., Peterzol, A., Fiedler, S., Ponchut, C., Labiche, J. C. & Bravin, A. (2006). J. Synchrotron Rad. 13, 260-270.]).

2.6.2. Dose measurement protocol

Measurements of the dose deposited at different positions within a cylindrical plastic phantom were performed using a CT ionization chamber (Radcal 10X6-3CT, Monrovia, USA) with a read-out connection to a Radcal dosimeter. The Radcal chamber is designed for CT dose measurements and provides dose values in water. The phantom consists of an external layer of polymethylmethacrylate (PMMA) with an outer diameter D = 100 mm and height of 150 mm, and an internal cylinder of polyethylene (PE) of 60 mm diameter (Fig. 4[link]). Three holes, 13 mm in diameter, make it possible to insert the ionization chamber at different positions (0, 18, 40 mm from the centre). The experiment was carried out with a monochromatic X-ray beam of 60 keV. The phantom was rotated by 360° around its central axis (perpendicular to the beam direction, see arrow in Fig. 4[link]) at a speed of 5° s−1. The X-ray beam used had a parallelepipedal shape with a cross section equal to 2 mm (V) × 100 mm (H). The holes were filled by cylinders of PE, except for the one containing the chamber.

[Figure 4]
Figure 4
Phantom geometry and experimental configuration used in the CT measurements.
2.6.3. Simulation parameters

Simulations were carried out using Geant4 (version 9.5) with the PENELOPE physics models. An ad hoc electron cut-off value of 100 keV was used in order to avoid tracking any secondary electron. The distance between the source and the sample was set to 50 cm. The measurements performed with the Radcal chamber on the CT cylindrical phantom were simulated with the active volume of the ionization chamber, 3 cm3, being replaced by a water volume (the chamber provides dose-to-water values). For the simulation, a voxel size of 0.2 mm × 0.2 mm × 1 mm was chosen and the dose estimates were computed by averaging the values obtained in all the voxels contained in the volume of the chamber's active region. The number of histories chosen for the simulations was 109 in order to obtain a mean statistical uncertainty below 1% in the central irradiated slice

2.7. Performances comparison between TLE and MC dose deposition simulations in anatomically complex cases

In order to assess the potential of the TLE method to produce dose distributions in complex cases and to evaluate the necessary computing time, two simulations were performed using experimental CT images as inputs for defining the geometry and composition of the volume in which the dose deposition would be evaluated. We considered the CT data of an excised and formalin-fixed human knee joint slice (444 × 451 × 1 voxels, with voxel size 0.3 mm × 0.3 mm × 10 mm) and of a tumour-bearing human breast slice (321 × 320 × 1 voxels, with voxel size 0.3 mm × 0.3 mm × 10 mm) acquired at 60 keV via high-resolution CT imaging (Snigirev et al., 1995[Snigirev, A., Snigireva, I., Kohn, V., Kuznetsov, S. & Schelokov, I. (1995). Rev. Sci. Instrum. 66, 5486-5492.]; Förster et al., 1980[Förster, E., Goetz, K. & Zaumseil, P. (1980). Krist Technik, 15, 937-945.]). The CT images were acquired experimentally using a detector with a pixel size of 50 µm.

The parallel X-ray beam had a rectangular cross section of 10 mm × 136 mm in the knee joint measurement and of 10 mm × 97 mm in the breast measurement. The object was rotated around an axis perpendicular to the beam direction and 200 angular projections over 360° were simulated. The CT volume of the knee joint was segmented using a simple thresholding method in order to attribute each voxel to a specific material (solution of 4% formalin in water, muscle and cortical bone). The breast was segmented into multiple components using a marker-controlled watershed viscous transform (Vachier & Meyer, 2005[Vachier, C. & Meyer, F. (2005). J. Math. Imaging Vis. 22, 251-267.]). By assuming that a three-dimensional image is composed of reliefs, this method is well adapted to the intrinsic nature of the CT images which are characterized by a strong signal on the borders of each feature. The medical a priori knowledge of the sample composition was used to assign the following materials to the different regions: adipose tissue, glandular tissue, air and 4% formalin solution in water. In both cases, 108 histories were numerically tracked. The simulation parameters were the same as those described in §2.6.3[link].

In order to test the validity of the local deposition assumption of all energy of the produced secondary electrons, one of the simulations for the breast case was repeated with 107 events by activating the tracking of the secondary electrons. The resulting deposited dose is presented in the next section.

3. Results

3.1. Benchmarking against MC simulation

In the test case described in §2.5[link], we first verified that our TLE implementation in GATE leads to the same dose distribution computed by the equivalent MC dose actor available in GATE. This is in fact the case within the statistical fluctuations, as shown in Figs. 5[link] and 6[link]. The integrated deposited energy values calculated with the two methods are in agreement within the given uncertainty (Table 2[link]). The uncertainties are calculated considering the standard deviation of the integrated value of the deposited energy using 20 simulations.

Table 2
Integrated energy with relative uncertainty in the thorax model, computed with the TLE and the MC method, respectively, under the same simulation conditions (4 × 108 incident photons)

MC TLE Relative difference
3.04 × 106 ± 0.009% meV 3.04 × 106 ± 0.008% meV 0.064%
[Figure 5]
Figure 5
Dose distribution maps obtained in the thorax model with the MC method (a) and the TLE method (b). 4 × 108 incident photons were used. The black oval line in (a) delimits the region considered for the comparison of the total deposited energy shown in Table 2[link].
[Figure 6]
Figure 6
Energy deposition profiles in the thorax model simulated with both TLE and MC methods (4 × 108 incident photons). The profiles correspond to the horizontal section marked by the green line in Fig. 5[link].

In Fig. 5[link] and in the profiles in Fig. 6[link] it can be seen that the uncertainty is much larger in the case of the MC simulation because of the slower convergence, as compared with the TLE method.

3.2. Accuracy test of the TLE method against measurements

The flux values measured with the PTW chamber and the HEGD detector, reported in Table 3[link], are in satisfactory agreement. From now on, only the values measured with the PTW are taken into account to assess the photon flux necessary for the simulations.

Table 3
Photon flux values obtained from the measurements at different energies with the HEGD and the PTW chamber

Energy (keV) HEGD [photons (mm2 mA s)−1] PTW [photons (mm2 mA s)−1]
33.7 (1.82 ± 0.14) × 107 (1.84 ± 0.02) × 107
45.0 (1.27 ± 0.06) × 108 (1.26 ± 0.01) × 108
50.5 (1.74 ± 0.09) × 108 (1.71 ± 0.02) × 108
65.0 (1.69 ± 0.09) × 108 (1.60 ± 0.02) × 108
72.0 (1.41 ± 0.07) × 108 (1.33 ± 0.01) × 108

The uncertainties on the HEGD values were calculated using the statistical uncertainties on repeated measurements, while for the PTW the uncertainties on both the measurements and the calibration were considered.

In Table 4[link], the results of the comparison between measured and simulated dose for the CT phantom are reported. The error bars on the simulated values are calculated propagating the uncertainty on the fluence and the statistical error of the simulations. For the measured values, uncertainties are obtained again by propagating the statistical spread of repeated measurements and that of the ionization chamber calibration (this value is provided by the chamber constructor). In Fig. 7[link], the dose distributions calculated in one transversal slice of the CT phantom are shown, with the CT ionization chamber successively positioned in the three holes. The brighter circular area in the phantom is the volume filled with water (corresponding to the PTW ionization chamber) for which the dose was both measured and simulated.

Table 4
Comparison of simulated and measured doses (in Gy) for the CT geometries of Fig. 7[link]

The ionization chamber positions refer to Fig. 4[link]. For each position the computations and the measurements were repeated three times because for each measurement the irradiation time of the chamber is slightly different, thus leading to a different number of photons that have to be entered in the simulations.

Position Measured (Gy) Simulated (Gy) % Difference
1 0.262 ± 0.010 0.259 ± 0.007 1.1
1 0.261 ± 0.010 0.259 ± 0.007 0.8
1 0.261 ± 0.010 0.259 ± 0.007 0.8
2 0.264 ± 0.010 0.265 ± 0.007 −0.3
2 0.264 ± 0.010 0.265 ± 0.007 −0.3
2 0.263 ± 0.010 0.265 ± 0.007 −0.7
3 0.261 ± 0.010 0.266 ± 0.007 −1.9
3 0.260 ± 0.010 0.266 ± 0.007 −2.3
3 0.261 ± 0.010 0.266 ± 0.007 −1.9
[Figure 7]
Figure 7
Simulated dose distribution maps of a transversal section of the CT phantom corresponding to three different positions of the ionization chamber (the water-filled volume is the white circular area) as reported in Fig. 4[link]. The grey levels correspond to dose values in Gy. The holes not hosting the ionization chamber are filled with PE.

3.3. Performances comparison in anatomically complex cases

The results of the simulations carried out with the high-resolution CT volumes of the human knee joint and breast are presented in Fig. 8[link]. Relatively smooth dose maps (∼0.3% mean statistical error per voxel in the whole knee joint volume and ∼0.2% in the breast volume) are obtained using the TLE (with 108 histories) whereas, for the same number of incident photons, equivalent MC simulations give results affected by very high uncertainty. In fact, for the MC case a large number of voxels are left with no dose information at all and this leads to a mean statistical uncertainty per voxel of ∼6.8% in the knee joint and ∼5.4% in the breast. Contrary to the standard MC, the use of TLE makes it possible to obtain meaningful dose maps also with a very small number of events (as shown in Fig. 8[link]). For the sake of comparison, integration over all the voxels composing the sample has been performed to work out the average values of the deposited dose. The difference between the average dose values simulated with the standard MC and the TLE method is ∼1.2% in the case of the knee joint image and ∼0.6% for the breast test case.

[Figure 8]
Figure 8
Dose distributions obtained in the human knee joint and breast with the MC (centre) and the TLE method (right), with 108 incident photon counts. (a), (b) and (c): knee joint; (d), (e) and (f): breast. (a), (d) are the sample models derived from experimental CT data, (b), (e) are the dose maps obtained with the TLE method, while (c), (f) are the maps computed with the equivalent MC method. In (a), the region marked (1) has been assigned to the formalin solution, (2) to the adipose tissue, (3) to muscle and (4) to the cortical and trabecular bone. In (d), (5) has been assigned to glandular tissue, (6) represents the adipose tissue, (7) is the tumour and (8) is the formalin solution. The dose scales in (b), (c), (e) and (f) are normalized to unity.

The total computation time was about 11 h for the knee and 9 h for the breast, using ten 3.3 GHz CPUs for the simulations based on the TLE method, whereas, using the MC method, 10 h are required for the knee and 8 h for the breast. These differences in time are due to several additional operations required within the TLE and related to its coefficients database (see §2.3[link]).

It is possible to estimate the number of events, Nf, required to obtain a statistical error of about σf, by means of Nf = (σi/σf)2Ni, where σi and σf are the standard deviations with Ni and Nf events, respectively. It is possible to estimate this number in order to recover the same average statistical error as a TLE calculation for a considered region of interest (i.e. the entire volume) (Table 5[link]).

Table 5
Performance comparison in anatomically complex cases

Sample TLE uncertainty MC uncertainty Nf/Ni
Knee 0.3% 6.8% ∼460
Breast 0.2% 5.4% ∼650
†The value of Nf/Ni has been corrected considering that the TLE method requires about 10% more computation time with respect to the MC with the same number of histories. This value can be easily related to the gain in speed of the simulation because of the linear correlation between the computation time and the number of histories used in the simulation.

The reason for this very long time is the great number and the small size of voxels in the two presented cases.

In Table 6[link] the values of the integrated deposited energy calculated with and without forcing the secondary electrons to deposit locally their energy are reported for the breast case scoring 107 events: the two estimates are in good agreement within the statistical uncertainty.

Table 6
Integrated energy deposition calculated with (asterisk) and without switching the tracking option of the secondary electrons

Values are expressed in MeV.

TLE* TLE Relative difference
7.98 × 105 ± 0.015% 7.98 × 105 ± 0.015% 0.033%

4. Concluding remarks

In the present study the implementation of the TLE method in the next release of the open source GATE/Geant4 9.5 simulation platform is described. As a result, a new software module named TLEDoseActor was created, which uses tables of energy-absorption coefficients generated from the EPDL97 database. It is also shown that TLE provides a powerful tool for simulating the dose deposited by photon beams in the keV energy region, when secondary electrons can be deemed to deposit all their energy locally, using a small number of events and without loss of accuracy. Potential applications of this new tool range from CT imaging to nuclear medicine as well as low-energy X-ray external radiotherapy.

The GATE implementation of the TLE was benchmarked against analogous MC dose calculations. No significant difference was found between the results of the two methods referring to the total deposited energy. In addition, the spatial distributions of deposited energy were also found to be in agreement within statistical fluctuations. However, the TLE method provides the substantial advantage of exhibiting a strong variance reduction in comparison with equivalent MC simulations, with no significant loss of accuracy. In an anatomically complex case (cf. §3.3[link]) we found that MC computations required about 500 times more events to reach the same averaged statistical uncertainty per voxel as the TLE method, using the same parameters.

The TLE dose calculations were also validated experimentally in the case of a CT dosimetry experiment with a cylindrical phantom irradiated by a 60 keV monochromatic synchrotron beam. Differences between experiment and simulation lower than 2.3% were found.

The applicability of the TLE method is, however, currently limited to an energy range between 1 keV and 1 MeV owing to the limited range of the pre-calculated mass energy-absorption coefficients in the database.

Most importantly, the high-energy limit for applying the TLE approach is fixed by the assumption that the electron transport can be neglected (local energy deposit), with no significant energy escape in the form of secondary radiation, i.e. energy-transfer and energy-absorption coefficients are assumed to coincide, this limitation becoming more stringent in the case of high-Z elements.

In summary, the TLE method makes it possible to simulate dose distributions in very complex geometries, like organ volumes recorded during high-resolution CT scans (with the identification of the atomic composition of the tissue assigned to each voxel) in several hours, depending on the parameters used, whereas analogous MC simulations requiring equivalent amounts of processing time suffer from much slower statistical convergence. The detailed statistical study of the variance reduction behaviour, which depends on geometric (voxel size), beam configuration and material parameters (density and composition), will be the subject of a future investigation.

Acknowledgements

The authors acknowledge the support through the DFG, Cluster of Excellence Munich-Centre for Advanced Photonics (EXE158) and the provision of beam time and laboratories by the ESRF as well as the partial support of the Rhône-Alpes Research Program in hadron therapy (PRRH Etoile) and the European collaboration Envision (grant agreement no. 241851). They also would like to thank Dr H. Requardt, T. Brochard and Dr C. Nemoz for the technical support and electronic, computing and software assistance at the ESRF biomedical beamline.

References

First citationAgostinelli, S. et al. (2003). Nucl. Instrum. Methods Phys. Res. A, 506, 250–303.  Web of Science CrossRef CAS Google Scholar
First citationAllison J. et al. (2006). IEEE Trans. Nucl. Sci. 53, 270–278.  Google Scholar
First citationAttix, F. H. (2004). Introduction to Radiological Physics and Radiation Dosimetry. Weinheim: Wiley-VCH.  Google Scholar
First citationBattistoni, G., Cerutti, F., Fassò, A., Ferrari, A., Muraro, S., Ranft, J., Roesler, S. & Sala, P. R. (2007). AIP Conf. Proc. 896, 31–49.  CrossRef CAS Google Scholar
First citationBerger, M. J. et al. (2005). National Institute of Standards and Technology, Gaithersburg, MD, USA.  Google Scholar
First citationBerger, M. J., Hubbell, J. H., Seltzer, S. M., Chang, J., Coursey, J. S., Sukumar, R., Zucker, D. S. & Olsen, K. (2010). XCOM: Photon Cross Section Database (version 1.5), https://physics.nist.gov/xcomGoogle Scholar
First citationCarlsson, G. A. (1985). The Dosimetry of Ionizing Radiation, pp. 1–75. Orlando: Academic Press.  Google Scholar
First citationChibani, O. & Williamson, J. F. (2005). Med. Phys. 32, 3688–3698.  Web of Science CrossRef PubMed Google Scholar
First citationCoan, P., Peterzol, A., Fiedler, S., Ponchut, C., Labiche, J. C. & Bravin, A. (2006). J. Synchrotron Rad. 13, 260–270.  Web of Science CrossRef IUCr Journals Google Scholar
First citationCullen, D. E., Hubbell, J. H. & Kissel, L. (1997). EPDL97: The Evaluated Photon Data Library, UCRL-50400, Vol. 6, Rev. 5, https://www-nds.iaea.org/epdl97/document/epdl97.pdfGoogle Scholar
First citationDemarco, J. J., Wallace, R. E. & Boedeker, K. (2002). Phys. Med. Biol. 47, 1321–1332.  Web of Science CrossRef PubMed Google Scholar
First citationFerrari, A., Sala, P. R., Fasso, A. & Ranft, J. (2005). FLUKA Technical Report, CERN-2005–10, INFN/TC 05/11, SLAC-R-773.  Google Scholar
First citationFörster, E., Goetz, K. & Zaumseil, P. (1980). Krist Technik, 15, 937–945.  Google Scholar
First citationFragoso, M., Kawrakow, I., Faddegon, B. A., Solberg, T. D. & Chetty, I. J. (2009). Med. Phys. 36, 5451–5466.  Web of Science CrossRef CAS PubMed Google Scholar
First citationFreud, N., Letang, J. M., Mary, C., Boudou, C., Ferrero, C., Elleaume, H., Bravin, A., Esteve, F. & Babot, D. (2008). IEEE Trans. Nucl. Sci. 55, 1008–1017.  Web of Science CrossRef Google Scholar
First citationICRU (1984). ICRU Report 37. International Commission on Radiation Units and Measurements, Bethesda, MD, USA.  Google Scholar
First citationICRU (1989). ICRU Report 44. International Commission on Radiation Units and Measurements, Bethesda, MD, USA.  Google Scholar
First citationIvanchenko, V. et al. (2011). Prog. Nucl. Sci. Technol. 2, 898–903.  Google Scholar
First citationJan, S. et al. (2004). Phys. Med. Biol. 49, 4543–4561.  Web of Science CrossRef PubMed CAS Google Scholar
First citationJan, S. et al. (2011). Phys. Med. Biol. 56, 881–901.  Web of Science CrossRef CAS PubMed Google Scholar
First citationJenkins, T. M., Nelson, W. R. & Rindi, A. (1988). Monte Carlo Transport of Electrons and Photons. New York: Plenum Press.  Google Scholar
First citationPelowitz, D. B. (2005). Editor. MCNPX User's Manual, p. 473. Los Alamos National Laboratory, New Mexico, USA.  Google Scholar
First citationPrezado, Y., Vautrin, M., Martínez-Rovira, I., Bravin, A., Estève, F., Elleaume, H., Berkvens, P. & Adam, J. F. (2011). Med. Phys. 38, 1709–1717.  Web of Science CrossRef CAS PubMed Google Scholar
First citationSalvat, F. & Fernández-Varea, J. M. (2009). Metrologia, 46, S112–S138.  CrossRef CAS Google Scholar
First citationSalvat, F., Fernández-Varea, J. M. & Sempau, J. (2011). PENELOPE. A Code System for Monte Carlo Simulation of Electron and Photon Transport. Issy-les-Moulineaux: OECD Nuclear Energy Agency.  Google Scholar
First citationSmans, K., Zoetelief, J., Verbrugge, B., Haeck, W., Struelens, L., Vanhavere, F. & Bosmans, H. (2010). Med. Phys. 37, 2082–2091.  Web of Science CrossRef PubMed Google Scholar
First citationSnigirev, A., Snigireva, I., Kohn, V., Kuznetsov, S. & Schelokov, I. (1995). Rev. Sci. Instrum. 66, 5486–5492.  CrossRef CAS Web of Science Google Scholar
First citationTaylor, R. E., Yegin, G. & Rogers, D. W. (2007). Med. Phys. 34, 445–457.  Web of Science CrossRef PubMed CAS Google Scholar
First citationVachier, C. & Meyer, F. (2005). J. Math. Imaging Vis. 22, 251–267.  Web of Science CrossRef Google Scholar
First citationWilliamson, J. F. (1987). Med. Phys. 14, 567–576.  CrossRef CAS PubMed Web of Science Google Scholar
First citationYe, S. J., Ove, R. & Naqvi, S. A. (2006). Health Phys. 91, 361–366.  Web of Science CrossRef PubMed CAS Google Scholar
First citationZee, W. van der, Hogenbirk, A. & van der Marck, S. C. (2005). Phys. Med. Biol. 50, 625–641.  Web of Science CrossRef PubMed Google Scholar

© International Union of Crystallography. Prior permission is not required to reproduce short quotations, tables and figures from this article, provided the original authors and source are cited. For more information, click here.

Journal logoJOURNAL OF
SYNCHROTRON
RADIATION
ISSN: 1600-5775
Follow J. Synchrotron Rad.
Sign up for e-alerts
Follow J. Synchrotron Rad. on Twitter
Follow us on facebook
Sign up for RSS feeds