Calculation of clinical dose distributions in proton therapy from microdosimetry

This article has been accepted for publication and undergone full peer review but has not been through the copyediting, typesetting, pagination and proofreading process, which may lead to differences between this version and the Version of Record. Please cite this article as doi: 10.1002/MP.13861 This article is protected by copyright. All rights reserved Calculation of clinical dose distributions in proton therapy from microdosimetry Bertolet A.1,2, Cortés-Giraldo M.A.2, Souris K.3, Cohilis M.3, Carabe A.1 1Department of Radiation Oncology, Hospital of The University of Pennsylvania, Philadelphia, PA USA 2Department of Atomic, Molecular and Nuclear Physics, Universidad de Sevilla, Seville, Spain 3Université Catholique de Louvain, Center for Molecular Imaging and Experimental Radiotherapy, Louvain, Belgium. Corresponding author: Alejandro Carabe-Fernandez. Department of Radiation Oncology, Hospital of The University of Pennsylvania, Philadelphia, 19104, PA, USA. E-mail: [email protected]; [email protected] Accepted Article This article is protected by copyright. All rights reserved Abstract Purpose: To introduce a new algorithm –MicroCalcfor dose calculation by modeling microdosimetric energy depositions and the spectral fluence at each point of a particle beam. Proton beams are considered as a particular case of the general methodology. By comparing the results obtained against Monte Carlo computations we aim to validate the microdosimetric formalism presented here and in previous works. Methods and material: In previous works, we developed a function on the energy for the average energy imparted to a microdosimetric site per event and a model to compute the energetic spectrum at each point of the patient image. The number of events in a voxel is estimated assuming a model in which the voxel is completely filled by microdosimetric sites. Then, dose at every voxel is computed by integrating the average energy imparted per event multiplied by the number of events per energy beam of the spectral distribution within the voxel. Our method is compared with the Proton Convolution Superposition (PCS) algorithm implemented in Eclipse™ and the fast Monte Carlo code MCsquare, which is here considered the benchmark, for in-water calculations, using in both cases clinically validated beam data. Two clinical cases are considered: a brain and a prostate case. Results: For a SOBP beam in water, the mean difference at the central axis found for MicroCalc is of 0.86% against 1.03% for PCS. Three-dimensional gamma analyses in the PTVs compared with MCsquare for criterion (3%, 3 mm) provide gamma index of 95.07% with MicroCalc versus 94.50% with PCS for the brain case and 99.90% versus 100.00%, respectively, for the prostate case. For selected organs at risk in each case (brainstem and rectum), mean and maximum difference with respect to MCsquare dose are analyzed. In the brainstem, mean differences are 0.25 Gy (MicroCalc) vs 0.56 Gy (PCS), while for the rectum these values are 0.05 Gy (MicroCalc) vs 0.07 Gy (PCS). Conclusions: The accuracy of MicroCalc seems to be, at least, not inferior to PCS, showing similar or better agreement with MCsquare in the considered cases. Additionally, the algorithm enables simultaneous computation of other quantities of interest. These results seem to validate the microdosimetric methodology in which the algorithm is based on. Keywords: microdosimetry, dose calculation, proton therapy, dose algorithms Accepted Article This article is protected by copyright. All rights reserved Introduction Microdosimetry is the study of the energy deposited by charged particles in microscopic structures. This framework was originally developed to evaluate the effects of radiation on cells and sub-structures of cells1–3. Due to the stochastic nature of the interaction between charged particles, microdosimetry theory deals with distributions, rather than deterministic values, of energy imparted to such structures4. Lineal energy (𝑦) is a microdosimetric quantity, defined as the quotient between the energy imparted to a certain structure, called site, and the mean chord length of the particle in that site5. The average of the weighted distribution of this quantity is related to the macroscopic concept of Linear Energy Transfer (LET)6. Finally, LET is considered as one of the factors impacting the relative biological effectiveness (RBE) for each kind of radiation7–9. On the other hand, specific energy (𝑧) is a different stochastic quantity, defined as the energy imparted to the site divided by the site mass. The average of the microdosimetric distribution of 𝑧, coincides with the macroscopic concept of absorbed dose10. Traditionally, distributions of specific energies for different types of radiation have been applied, for instance, to obtain radiobiological parameters defining the effect of a given test radiation11 or to characterize detectors responses12. However, it has never been used to calculate three-dimensional dose distributions in a treatment planning system (TPS) for ion therapy. We present a general methodology to analytically calculate dose based on specific energy, particularized to the proton radiotherapy case. Methods and Materials Starting from previous works13,14, we have created an analytical model for the dependency of the frequency average of the single-event distribution of specific energy 𝑧𝐹 on the proton kinetic energy, i.e., a function 𝑧𝐹= 𝑧𝐹(𝐸). The subscript 𝐹 indicates frequency-average, i.e., the average of the frequencies of appearances for each possible specific energy value in a single event. This subscript is traditionally included in contraposition to the dose-weighted average, 𝑧𝐷, which results from weight each specific energy value for a single event according to the dose imparted during that event. 𝑧𝐹(𝐸) needs to be integrated over the number of events occurred for each proton energy at each Accepted Article This article is protected by copyright. All rights reserved point of the calculation grid. Therefore, a dose distribution can be calculated as long as the information of the number of events per energy is known. Definitions in microdosimetry theory Given a beam of charged particles passing through matter, an event in microdosimetry is the process in which each independent particle track deposits energy by interacting with the atoms of the matter15. The event, by definition, includes the depositions of energy from the secondary particles (mainly electrons) generated along the particle track. A site is a sub-region of the matter traversed by the radiation beam in which energy deposits are scored. Sites can take different shapes and sizes and allows for characterizing the way in which energy is deposited spatially. The sum of all energy deposits in a site per each independent particle track is called single-event energy imparted, 𝜀𝑠, and it is a stochastic quantity, whose actual value for each event depends on several factors16. Therefore, after a number of events, it is possible to obtain the single-event distribution of energy imparted, in terms of its probability density function, 𝑓(𝜀𝑠), by counting the relative frequency of events with each value of energy imparted. The quotient between the energy imparted 𝜀 after a number of events and the mass 𝑚 of the site is called specific energy, 𝑧=𝜀/𝑚. Then, the probability density function of the single-event distribution of specific energy, 𝑓1(𝑧), can be obtained from 𝑓(𝜀𝑠). The average of the single-event distribution of specific energy is noted by 𝑧𝐹, where 𝐹 stands for frequency. For a given number of events, the specific energy for a site is obtained as the sum of the specific energies at each event. Therefore, as 𝑧 is also a stochastic quantity, if we consider different series of events, with average 𝜈, the distribution of the total specific energy can be expressed in terms of its probability density function, 𝑓(𝑧;𝜈). The dose absorbed by the site is not a stochastic quantity and coincides with the average of 𝑓(𝑧;𝜈), i.e. 𝐷=∫𝑑𝑧 𝑓(𝑧;𝜈)=𝑧. (1) It is more interesting to consider the distribution of specific energy for a given dose 𝐷, 𝑓(𝑧;𝐷), rather than in terms of the average number of events 𝜈. This distribution gives the probability of obtaining a specific energy 𝑧 after a series of events when the absorbed dose is 𝐷. This distribution can be obtained from the single-event distribution of specific Accepted Article This article is protected by copyright. All rights reserved energy, 𝑓1(𝑧), since the interactions between particles and matter can be described by a compound Poisson process10. Hence, the absorbed dose can be expressed in terms of the average of the single-event distribution of specific energy as 𝐷=𝜈 𝑧𝐹. (2) Accordingly, if the number of events occurring in a site and the mean specific energy per event are known, it is possible to calculate the corresponding absorbed dose. Spectral decomposition Equation (2) is valid for any beam regardless its particles or spectrum. However, it is challenging to analytically model 𝑧𝐹 for such a broad range of possibilities. A good approach consists of considering the single-event distributions of specific energy for monoenergetic beams of energy 𝐸, with average 𝑧𝐹=𝑧𝐹(𝐸). Then, if the number of events 𝜈(𝐸) from particles with energy 𝐸 is known, the absorbed dose can be calculated as 𝐷=∫𝜈(𝐸) 𝑧𝐹(𝐸) d𝐸. (3) Model for 𝑧𝐹(𝐸) As said above, the single-event specific energy is related to 𝜀𝑠. Equivalently, the average of its single-event distribution is given by 𝑧𝐹(𝐸)=𝜀𝑠(𝐸) 𝑚. (4) The average energy imparted to a site at each particle energy depends on three factors: (a) the stopping power of the particle, (b) how long is the particle track within the site -segment length-, and (c) the probability that secondary electrons generated by the particle leave the site. The trial function for 𝜀𝑠(𝐸) should be, firstly, proportional to the stopping power, which in turn depends on the particle energy. It also should contain a term to consider the losses of energy imparted when the range of secondary electrons is larger than the site dimensions. And finally, it needs to decrease as the particle track becomes shorter. In previous works we have explored how to use Monte Carlo simulations to produce microdosimetric data13 and what functions can be used to characterize microdosimetric quantities14. A function meeting all the referred requirements is: Accepted Article This article is protected by copyright. All rights reserved 𝜀𝑠(𝐸)=𝐶·erf(𝑘𝐸𝑞)·log(ℎ𝑃,𝑚 𝐸+ 𝑏) · 𝛼𝑃,𝑚 ·𝐸(𝑝𝑃,𝑚−1)+𝑒 𝛼𝑃,𝑚·𝐸𝑝𝑃,𝑚 +2𝐷𝑠/3, (5) where 𝐶, 𝑏 and 𝑒 are free parameters of the function whereas ℎ𝑃,𝑚,𝛼𝑃,𝑚 and 𝑝𝑃,𝑚 depend on the considered particle and medium and 𝐷𝑠 depends on the site dimension. The mass of the site can be obtained from its volume and the density of the considered medium. Since dose to water is the quantity normally used in a TPS for dose calculations, the site is considered to be made of liquid water. Number of events per energy, 𝜈(𝐸): a model for the voxel-site relation By definition, the distribution of average number of events per energy 𝜈(𝐸) should coincide with the spectral distribution at the site, i.e., the number of particles of each energy impinging on the site. However, we are normally interested in calculating dose in a voxelized geometry, being the voxel considerably larger than the typical microscopic sites. Therefore, we need to firstly consider the number of events occurring in the voxel. For this, we assume the voxel is represented by an arrangement of sites in a rectangular grid, as Figure 1 shows for the case of spherical sites. Let 𝑁𝑡 be the number of sites in the voxel contained on a layer perpendicular to the particle track, where 𝑡 stands for transversal. Then, if the transversal cross section of a site is 𝜎𝑠, the total cross section in which events are produced is 𝑁𝑡 𝜎𝑠 so that the number of events in the voxel is proportional to Φ𝐸(𝐸)𝑁𝑡 𝜎𝑠, where Φ𝐸(𝐸) is the spectral fluence of particles. Note that this is valid as long as we can assume that the number of indirect events (i.e. produced just by secondary particles) is negligible with respect to the number of direct events (i.e. produced by primary protons), which is expected to be the case for any voxel with dose relevant in radiotherapy. An additional aspect needs to be taken into consideration. Let 𝑑𝑉 be the dimension of the voxel parallel to the direction of the particle tracks and 𝑑𝑠 the maximum dimension of the site in the same direction. Then, the number of sites aligned within the voxel will be 𝑁𝑙=𝑑𝑉/𝑑𝑠, where 𝑙 stands for longitudinal. If the incidental particle energy is such that its range is longer than 𝑑𝑉, then it will exit the voxel producing an event in each site on its way. Therefore, the number of events will be 𝜈=𝑁𝑙. However, if the range of the particle is shorter than 𝑑𝑉, the number of events produced by that particle is fewer since the proton stops before reaching the final sites. We assume, in that case, that the number of events Accepted Article This article is protected by copyright. All rights reserved produced by that particle is proportional to its CSDA range, 𝑅(𝐸), so the number of events in the voxel per energy, 𝜈𝑉(𝐸), is: 𝜈𝑉(𝐸)=Φ𝐸(𝐸)𝑁𝑡 𝜎𝑠 max(𝑑𝑉 𝑑𝑠,⌈𝑅(𝐸) 𝑑𝑠⌉), (6) where ⌈𝑥⌉ means the least integer greater than or equal to 𝑥, since the number of events is an integer quantity. However, according to equation (2), we are interested in the average number of events per site rather than in the complete voxel. As the total number of sites in the voxel is 𝑁𝑡 𝑁𝑙, the average number of events per site is given by 𝜈(𝐸)=𝜈𝑉(𝐸) 𝑁𝑡 𝑁𝑙=1 𝑁𝑙Φ𝐸(𝐸) 𝜎𝑠max(𝑑𝑉 𝑑𝑠,⌈𝑅(𝐸) 𝑑𝑠⌉). (7) A particular case: protons in spherical sites All shown so far remains as valid regardless the particle species, and the site’s shape or dimension considered. Now, we introduce a particular case to illustrate the performance of the method. The analysis beyond this point is focused on proton clinical beams and spherical sites with diameter 𝑑. For this specific case, the longitudinal length of the site equals to the site diameter, 𝑑𝑠=𝑑, and so it does the site-dependent parameter on the expression for 𝜀𝑠 in equation (5): 𝐷𝑠=𝑑.14 In order to obtain the rest of parameters in equation (5), we performed simulations with Geant4-DNA17–20 of mono-energetic protons in spherical sites made of liquid water up to 100 MeV, which is the upper limit the code has for proton transport. To assess the influence of the site dimension in the calculation, we have considered three different cases for the site diameter: 1 μm, 5 μm and 10 μm. To calculate 𝜀𝑠(𝐸), for proton energies higher than 100 MeV we assume the validity of equation (5) can be extrapolated to those energies, as our trial function has a smooth behavior at such energies. Kernel for the spectral fluence of proton beams Once 𝜀𝑠, and therefore 𝑧𝐹, are obtained, the only quantity remaining to be determined is the spectral fluence density, Φ𝐸(𝐸). Again, this means the number of protons of each energy per unit area traversing each voxel in which the CT image of a patient is divided. This information can be obtained through Monte Carlo simulations of the proton Accepted Article This article is protected by copyright. All rights reserved transport in the patient. Some TPS are able to provide that information21,22. However, we have developed a kernel-based algorithm for the spectral fluence at each point given a proton arriving to the patient surface. Here we provide an outline of this, although it is fully presented in a separate work23. Different simulations of the transport of protons in water were carried out by using an application based on the Monte Carlo toolkit Geant424–26. Several values for the energy have been selected in the range of clinical use in proton radiotherapy: from 50 MeV up to 250 MeV. A distinction between primary and secondary protons is made to develop independent models for their transport in water. Primary protons are those which only undergo elastic collisions with the atoms of matter meanwhile secondary protons are those coming from inelastic collisions with nuclei. For each one, the model is composed by three elements: (i) given an initial fluence of protons Φ0 penetrating in water, the fraction of primary Φ𝑃𝑃(𝑑)/Φ0 and secondary protons Φ𝑆𝑃(𝑑)/Φ0 in a whole transversal plane at a given depth 𝑑; (ii) the lateral distribution of primary d𝜙𝑃𝑃(𝑟;𝑑)/d𝑟 and secondary protons d𝜙𝑆𝑃(𝑟;𝑑)/d𝑟 at each depth and each distance 𝑟 to the central axis of the beam; and (iii) the energy spectrum for primary d𝜙𝑃𝑃(𝐸;𝑟,𝑑)/d𝐸 and secondary d𝜙𝑆𝑃(𝐸;𝑟,𝑑)/d𝐸 at each lateral position and each depth. For the first case, we use the function proposed in the work of Ulmer27 for the primary protons, whereas the proportion of secondary protons at each depth Φ𝑆𝑃(𝑑) is modeled by a phenomenological function23: Φ𝑆𝑃(𝑑) Φ0=(𝑎𝑑−𝑏𝑑𝑐)⋅1 2(1+erf(𝑅(𝐸)−𝑑 √2 𝜏 )), (8) where 𝑎,𝑏 and 𝑐 are parameters to be fitted dependent on the beam energy, 𝑅(𝐸) is the CSDA range corresponding to the initial beam energy and 𝜏 is the same parameter than used by Ulmer27, related to the range straggling of the beam. For the lateral distribution, we assume that a single Gaussian centered at the beam axis is enough to model the spread in the fluence of primary protons, but not for the secondary ones. Hence, a Double Gaussian is proposed as a better function to model the latter case. Finally, for the sake of simplicity, the energy spectrum for both primaries and secondaries is modeled as a single Gaussian whose mean energy is given by the residual Accepted Article This article is protected by copyright. All rights reserved range (i.e., 𝑅(𝐸)−𝑑). Then, the kernel of spectral distribution at a voxel with entrance surface d𝑆 at the position (𝑑,𝑟) is given by Φ𝐸 Kernel(𝐸;𝑑,𝑟) d𝑆= ∑ Φ𝑖(𝑑) Φ0d𝜙𝑖(𝑟;𝑑) d𝑟 𝑖=𝑃𝑃,𝑆𝑃 d𝜙𝑖(𝐸;𝑟,𝑑) d𝐸 d𝑆, (9) where 𝑖 refers to protons and secondary protons. Note that both d𝜙𝑖(𝑟;𝑑)/d𝑟 and d𝜙𝑖(𝐸;𝑟,𝑧)/d𝐸 are normalized to the unit so that the argument in the summation provides the probability for a proton to have an energy 𝐸 at the position (𝑑,𝑟). Finally, to obtain the spectral fluence density for a broad beam it is necessary to convolve two-dimensionally the spectral fluence density of the broad beam coming onto the patient surface, Φ𝐸 Surface(𝐸;𝑟), with the kernel presented in equation (9): ΦE(𝐸;𝑑,𝑟)=∫Φ𝐸 Surface(𝐸;𝑟′) Φ𝐸 Kernel(𝐸;𝑑,𝑟−𝑟′)d𝑟′. (10) This model for the proton transport only considers liquid water. Transport through other media is here considered by replacing in equation (9) the depth in water, 𝑑, by the corresponding water-equivalent thickness 𝑑𝑊𝐸𝑇 resulting from the beam raytracing. Dose distributions comparisons All above have been implemented in a script for the Eclipse™ TPS v15.6 (Varian Medical Systems, Palo Alto, CA, USA) for Proton Beam Scanning (PBS) treatments, called hereinafter ‘MicroCalc’ or ‘µCalc’. This allows us for producing actual dose calculations in three-dimensional geometries. The TPS provides the information about the beam lateral width –to be convolved with our lateral model–. On the other hand, the number of protons prior to arrive to the patient surface can be estimated according to the number of monitor units (MU) for each spot21 so that we can obtain Φ𝐸 Surface(𝐸;𝑟). We present comparisons of dose distributions obtained by: (a) the microdosimetric method presented in equations (3-7) together with the kernel for the proton fluence summarized in equation (9); (b) the analytical Proton Convolution-Superposition (PCS) algorithm28 provided by Eclipse™; and (c) Monte Carlo simulations with the open source fast MC code MCsquare29, which is used as a benchmark in this work as it has been validated elsewhere30,31. Accepted Article This article is protected by copyright. All rights reserved scanning Monte-Carlo dose algorithm: Comparisons against measurements and simulations. Phys Med Biol. 2017;62(19):7659-7681. doi:10.1088/1361-6560/aa82a5 23. Bertolet A, Cortés‐Giraldo MA, Souris K, Carabe-Fernandez A. A kernel-based algorithm for fluence and spectral fluence in clinical proton beams to calculate doseaveraged LET and other dosimetric quantities of interest. Med Phys [SUBMITTED]. 2019. 24. Agostinelli S, Allison J, Amako K, et al. Geant4—a simulation toolkit. Nucl Instruments Methods Phys Res Sect A Accel Spectrometers, Detect Assoc Equip. 2003;506(3):250- 303. doi:10.1016/S0168-9002(03)01368-8 25. Allison J, Amako K, Apostolakis J, et al. Geant4 developments and applications. IEEE Trans Nucl Sci. 2006;53(1):270-278. doi:10.1109/TNS.2006.869826 26. Allison J, Amako K, Apostolakis J, et al. Recent developments in Geant4. Nucl Instruments Methods Phys Res Sect A Accel Spectrometers, Detect Assoc Equip. 2016;835:186-225. doi:10.1016/J.NIMA.2016.06.125 27. Ulmer W. Theoretical aspects of energy-range relations, stopping power and energy straggling of protons. Radiat Phys Chem. 2007;76(7):1089-1107. doi:10.1016/j.radphyschem.2007.02.083 28. Varian. Eclipse Proton Algorithms Reference Guide, v13.7. 2015. 29. Souris K, Lee JA, Sterpin E. Fast multipurpose Monte Carlo simulation for proton therapy using multi- and many-core CPU architectures. Med Phys. 2016;43(4):1700- 1712. doi:10.1118/1.4943377 30. Sorriaux J, Testa M, Paganetti H, et al. Experimental assessment of proton dose calculation accuracy in inhomogeneous media. Phys Medica. 2017;38:10-15. doi:10.1016/j.ejmp.2017.04.020 31. Huang S, Kang M, Souris K, et al. Validation and clinical implementation of an accurate Monte Carlo code for pencil beam scanning proton therapy. J Appl Clin Med Phys. 2018;19(5):558-572. doi:10.1002/acm2.12420 32. Bellinzona VE, Ciocca M, Embriaco A, et al. On the parametrization of lateral dose Accepted Article This article is protected by copyright. All rights reserved profiles in proton radiation therapy. Phys Medica. 2015;31(5):484-492. doi:10.1016/j.ejmp.2015.05.004 33. Grassberger C, Paganetti H. Elevated LET components in clinical proton beams. Phys Med Biol. 2011;56(20):6677-6691. doi:10.1088/0031-9155/56/20/011 Accepted Article This article is protected by copyright. All rights reserved FIGURE CAPTIONS Figure 1. Model for the internal structure of a voxel composed by spherical sites. Protons characterized by a given fluence density impinge on the frontal face of the voxel and particles interact with the row of sites along their track. Figure 2. Dose profiles at the central axis for single energy layers beams with energies 110 MeV, 160 MeV and 210 MeV calculated by employing the MicroCalc script, the PCS algorithm in Eclipse™ and the MC code MCsquare. Figure 3. Dose profiles at central axis for a spread-out Bragg peak with maximum range at 20 cm in depth and modulation of 10 cm calculated by employing the MicroCalc script, the PCS algorithm in Eclipse™ and the MC code MCsquare. Figure 4. Dose distributions on an axial plane for a brain case for (a) MCsquare, (b) MicroCalc and (c) PCS algorithm. The yellow contour depicts the PTV, with prescription dose of 54 Gy. On the color map, blue color represents low dose (near 0 Gy) while most intense red color represents 54 Gy. Differences between (d) MicroCalc and MCsquare and (e) PCS and MCsquare are also shown. In these color maps, most intense blue shows MCsquare dose is higher by 5 Gy, most intense red means MCsquare dose is lower by 5 Gy and translucent means no difference. Figure 5. Dose distributions on an axial plane at the middle of the PTV for a prostate case for (a) MCsquare, (b) MicroCalc and (c) PCS algorithm. The yellow contour depicts the PTV, with prescription dose of 68 Gy. On the color map, blue color represents low dose (near 0 Gy) while most intense red color represents 68 Gy. Differences between (d) MicroCalc and MCsquare and (e) PCS and MCsquare are also shown. In these color maps, most intense blue shows MCsquare dose is lower by 3.5 Gy, most intense red means MCsquare dose is higher by 3.5 Gy and white means no difference. LEGENDS Figure 2 µCalc. PCS algorithm. Accepted Article This article is protected by copyright. All rights reserved MCsquare. Figure 3 µCalc. PCS algorithm. MCsquare. Accepted Article mp_13861_f1.eps Thisarticleisprotectedbycopyright.Allrightsreserved Accepted Article Depth (cm) 5 10 15 20 25 30 Dose (a.u.) 0 0.5 1 1.5 2 2.5 3 3.5 110 MeV 160 MeV 210 MeV mp_13861_f2.eps Thisarticleisprotectedbycopyright.Allrightsreserved Accepted Article Depth (cm) 2 4 6 8 10 12 14 16 18 20 22 Dose (a.u.) 0 0.2 0.4 0.6 0.8 1 1.2 mp_13861_f3.eps Thisarticleisprotectedbycopyright.Allrightsreserved Accepted Article mp_13861_f4.tif Thisarticleisprotectedbycopyright.Allrightsreserved Accepted Article mp_13861_f5.tif Thisarticleisprotectedbycopyright.Allrightsreserved Accepted Article