A mechanistic relative biological effectiveness model-based biological dose optimization for charged particle radiobiology studies

The general-purpose Monte Carlo toolkit Geant4 (version 10.4) (Agostinelli 2003, Allison et al 2016) was used to perform the particle tracking that subsequently generated the physical and radiobiological quantities used for dose optimizations and RBE modeling. Many different typical physics lists such as 'QBBC', 'QGSP_BERT', 'QGSP_BIC' and 'FTFP_BERT' can be applied in particle therapy simulations which all provide electromagnetic and hadronic physics processes. The dose difference was found below 1% for particles within the therapeutic energy ranges using these physics lists (Geng et al 2018, Guan et al 2018). In this study, the 'FTFP_BERT' physics list was used as a representative. The general particle source (GPS) method was used to define the ion beams traveling along the central z axis to enter the target geometry. An 80  ×  80  ×  40 cm³ water phantom was built as the target for scoring quantities of interest for different ion beams. A scoring cylinder with a radius of 40 cm and thickness of 0.01 cm was built to cover the large size of the beam spot and to maximally capture the laterally scattered particles to subsequently score the so-called integral depth dose (IDD). The production cut of secondary particles (photons, electrons, positrons, protons and all recoil ions) was set to 0.01 cm to match the smallest length (thickness in this case) of the virtual detector, which followed the recommendation made by the Geant4 collaboration (Agostinelli 2003, Allison et al 2016). A derivative class of Geant4's primitive scorer was created to score the quantities of interest for each tracking step processed in the ProcessHits() function of the scorer by filling one-dimensional (1D) ROOT histograms (Brun and Rademakers 1997). In the current work, the abscissa of a 1D ROOT histogram is the depth in water and the ordinate can be dose, radiobiological parameters, or DNA damage yield.

For proton beams, 94 groups of energies that are used at The University of Texas MD Anderson Cancer Center Proton Therapy Center, ranging from 72.5 to 221.8 MeV with the penetration ranges (depths with 90% of peak dose in the distal falloff) of 4.0 to 30.6 cm in water, were selected for Monte Carlo simulations and dose optimizations. For ⁴He ions, the energy and range data were obtained using the National Institute of Standards and Technology ASTAR database (Berger 1995). A total of 161 groups of energies, ranging from 70 to 230 MeV/u with ranges of 4.1 to 33.1 cm in water, were selected. For ¹²C ions, the energy and range data were obtained from the Errata and Addenda: ICRU Report 73 (Sigmund et al 2009). A total of 161 groups of energies, ranging from 120 to 440 MeV/u with ranges of 3.6 to 32.0 cm in water, were selected. A double-layer ripple filter was modeled for carbon ions to broaden the Bragg peak. A total of 10⁷ source particle histories were simulated for each beamlet to make the simulation results, e.g. total dose, meet the statistical uncertainty requirement (relative error of the mean value  <1%) when the dose is larger than 5% of the peak dose.

In the mechanistic RMF model, reproductive cell death by mitotic catastrophe, apoptosis, or other cell death modes is explicitly linked to DNA double-strand break (DSB) induction and processing (Carlson et al 2008, Stewart et al 2018). An independent Monte Carlo Damage Simulation (MCDS, version 3.10A) (Stewart et al 2011) can be used to simulate ion- and energy-dependent DSB yields as input for the RMF model, although other Monte Carlo models (Nikjoo et al 1994, 1999, Friedland et al 1998, 2003, Alloni et al 2013) or experimental measurements could also potentially be used to generate input DSB yields. In the current study, the RMF model was applied to predict the biological effects of different ion beams considering all of the contributions from primary particles and their secondary particles.

Clonogenic cell survival was selected as the primary endpoint for calculating the RBE. Non-small cell lung cancer (NSCLC) H460 cells were selected as the target cells. Gamma rays from a ¹³⁷Cs source were selected as the reference photons. The reference LQ parameters ${{\alpha}_{x}}$ and ${{\beta}_{x}}$ are derived in our previous study (Guan et al 2015a) from cell survival data of H460 cells irradiated with ¹³⁷Cs photons and are used as the reference radiosensitivity parameters for RBE calculations in the current study.

The details of the RMF model and its application to calculating particle RBE can be found in Carlson et al (2008), Frese et al (2012) and Kamp et al (2015), and references therein. Here only some key steps are listed in applying the RMF mechanisms in Monte Carlo physics simulations to generate the spatial distributions of the radiobiological parameters needed for the next-step biological dose optimization.

• (1)  
  Estimating the cell-specific model constants $\kappa$ and $\theta$ using the physical and biophysical parameters for the specified cell type with the reference photons.
• (2)  
  Calculating radiosensitivity parameters ${{\alpha}_{i}}$ and ${{\beta}_{i}}$ for a charged particle (with the index i) of a given energy during the Monte Carlo particle tracking process.
• (3)  
  Calculating dose-averaged radiosensitivity parameters ${{\alpha}_{D}}$ and ${{\beta}_{D}}$ for a mixed field of radiation at a given spatial location (with the dose D).

The detailed calculation methods are described in the supplementary material (stacks.iop.org/PMB/64/015008/mmedia). After the Monte Carlo simulation, the spatial distributions of dose, DSB yield, ${{\alpha}_{D}}$, and ${{\beta}_{D}}$ in the water phantom for each beamlet can be obtained. In the current study, each of the above quantities is scored in an individual 1D ROOT histogram along the z axis of the phantom respectively. With these datasets, the RBE at different locations from a radiation field composed of multiple beamlets can be calculated.

The Python (version 3.4.3) programming language was used to develop the biological dose optimization algorithms using the NumPy and SciPy libraries. The purpose of biological dose optimization is to deliver a specified uniform RBE-weighted dose to the target region, which means that both the physical dose and the RBE should be considered. However, RBE is also a function of the physical dose and other physical and biological factors. Therefore, RBE-weighted dose optimization is usually a dynamic and complex calculation process involving multiple iterations until the predefined objective function is minimized. The developed Python code can be used to implement three functions:

• (1)  
  Solving beam weights and calculating the RBE in an iterative process.
• (2)  
  Predicting cell surviving fraction using the LQ cell survival model.
• (3)  
  Predicting the DSB yield.

The details of the optimization algorithm and the calculation process of biological effect are described in the supplementary material.

The biological dose optimization results for proton beams are shown in figure 1 as different subplots. In addition to the physical and biological dose distributions, figure 1(A) shows the spatial variation of RBE; figure 1(B) shows the dose-averaged α and β values; figure 1(C) shows the induced DSB including dose-normalized values in DSB/Gy/Gbp and absolute values in DSB/Gbp at the current dose level for each location; figure 1(D) shows the predicted cell survival.

Zoom In Zoom Out Reset image size

In figure 1(A), from the entrance to the proximal rise of the SOBP, the RBE shows a slightly decreasing trend from 1.12 to 1.10. This observation can be explained using equation (16) in the supplementary material. In figure 1(B), the dose-averaged α and β values are approximately constant; therefore, with all other parameters fixed in equation (16), RBE decreases with an increase in dose, which is the dose variation trend from the entrance to the proximal edge of the SOBP. Then, the RBE increases to 1.25 until the end of the SOBP and sharply increases at the distal edge to a maximum value of 1.72 (at the depth of 10.42 cm). In figure 1(B), the dose-averaged α and β values increase slightly within the SOBP but increase sharply at the end of the range. According to the equation (16), RBE increases with the increase of dose-averaged α and β values when other parameters are fixed. The RBE values and the variation trend in the current work are consistent with the original data reported by Frese et al using an analytical method (Frese et al 2012).

After the penetration range, the physical dose from proton beams is negligible. Therefore, it is not necessary to report the RBE values. It should be noted that the reported RBE at each location is calculated at its corresponding dose level, rather than at the identical cell survival rate everywhere. Therefore, these RBE values should not simply be compared because only at identical surviving fractions will the comparison of RBE with different beam qualities be indicative of the RBE for cell survival. Nevertheless, these reported RBE values can still be used in a qualitative way to show the increasing trend of RBE along the beam path.

The RMF model directly models underlying biological mechanisms such as DSB induction and processing to predict the RBE for cell survival, rather than relying on empirical relationships between experimental data and LET_d, for example, as done in some other phenomenological RBE models (Wilkens and Oelfke 2004, McNamara et al 2015, Polster et al 2015). The depth profile of DSB yield is also generated using equation (20) in the supplementary material. In figure 1(C), the dose-normalized DSB yield is shown to vary slowly along the beam path but increases sharply at the end of the SOBP and in the dose distal edge, which actually behaves similarly to the variation trend of LET_d (Paganetti 2014, Guan et al 2015b). This variation trend is consistent with the observations in our previous study using the ɣH2AX foci staining at the end of the range of a pristine proton beam (Guan et al 2015a).

In figure 1(D), the predicted surviving fraction within the SOBP is 0.1 as expected. The optimization results have thus provided a useful beam delivery strategy to investigate the biological effects along the SOBP formed by scanned beams. For example, an irradiation device can be designed to place the 12 columns of a 96-well cell culture plate at different locations so that the spatial variation of the biological effect along the SOBP can be sampled in a single irradiation run. The ten middle columns of the 96-well plate (12 columns  ×  8 rows) can be placed within the SOBP, the first column proximal and the last column distal to the SOBP, as marked in figure 1(D), to investigate biological effects both inside and outside of the target region. The principle of this experimental design using the high-throughput strategy to rapidly obtain radiobiological data has been described in our previous studies (Guan et al 2015a, Geng et al 2018).

Figure 1(E) shows the beam compositions used to form the proton SOBP. 36 groups of beam energies from 81.4 MeV (range  =  5 cm) to 118.6 MeV (range  =  10.1 cm) were used but the data from the 81.4 MeV beam are unseen owing to their relatively low contribution. In figure 1(F), the relative difference (calculated using equation (18) in the supplementary material) between the delivered biological dose after optimization and the prescribed value varies from  −1.0% to 0.7% for the whole target region but except for the target boundaries the relative difference is much lower and nearly within  ±0.3%, showing the high accuracy of the biological dose optimization.

In figures 1(B) and (C), the dose averaged α and β values, as well as the DSB yield beyond the dose distal edge, have non-zero values with large data fluctuations. They are actually the corresponding contributions from secondary particles. The dose compositions from different particle species can be analyzed to interpret these non-zero values beyond the distal edge of depth dose profile. In figure 2, the Geant4 calculated dose contributions from different species of a high-energy (207 MeV) proton beam are plotted. Beyond the dose distal falloff, the main dose contributors are secondary neutrons (essentially their recoil protons). The dose is at approximately 0.1% of the peak dose, which is negligible. These low-fluence recoil protons can result in large uncertainties of the calculated α and β values and DSB yield. Because the energy of the recoil protons is relatively low, their high-LET characteristics can result in large α and β values and DSB yield. However, it should be noted that these high α and β values and DSB yield do not indicate the absolute biological effect is high because the absolute dose level is negligible therein.

Zoom In Zoom Out Reset image size

In addition to the physical dose and biological dose, RBE, α and β, DSB yield, and cell survival of H460 cells from ⁴He ions with an SOBP ranging from 5 to 10 cm in water are shown in figures 3(A)–(D). The dose fragmentation tail beyond the SOBP is clearly evident. The mixed radiation field makes a quantitative description of the biological effect of heavier ion beams more complicated. In figure 3(A), RBE for ⁴He beams is higher than for proton beams. Proximal to the SOBP, the RBE varies slowly and approximately maintains a value of 1.2. Within the SOBP, the RBE increases from 1.3 to 1.8. In the distal edge, the RBE increases up to 2.3 and then decreases. In the fragmentation tail, the RBE continues to increase with the depth. The spatial variation trend of RBE values in figure 3(A) is consistent with the published data by Mairani et al (2016), while the numerical values are different because different target RWD values and different cell lines were used in these two studies.

Zoom In Zoom Out Reset image size

The above observations can be explained using dose and LQ parameters α and β as in equation (16) in the supplementary material, from which it is apparent that RBE is a monotonically decreasing function of dose but an increasing function of both α and β for the specific ion if all other parameters are fixed. As shown in figure 3(B), both α and β increase in the fragmentation tail but the dose continues to drop. Therefore, an increase in the RBE in the fragmentation tail is evident. Large fluctuations are evident for RBE and α in the fragmentation tail owing to the low particle fluence, but these fluctuations do not occur for other quantities of interest. It should be noted that the increased α and β values in the fragmentation tail are caused by the decrease in energy (increase in LET) of nuclear fragments along the penetration depth. The experimental data of Furusawa et al have demonstrated the varying trend of α and β with LET: they first increase and then decrease (Furusawa et al 2000). The values predicted using the RMF model are consistent with the gathered experimental data (Frese et al 2012) when LET is lower than 150–200 keV µm⁻¹. Independent investigations are underway to account for the spatial proximity and complexity of individual DSBs generated by low-energy ions in the MCDS model to account for the correction for an 'overkill' or 'proximity' effect at high ${{\left({{Z}_{ef\,f}}/\beta \right)}^{2}}$ values (Butkus et al 2018). For the DSB induction data in figure 3(C), the normalized DSB yield (per Gy per Gbp) increases within the SOBP, but the absolute DSB yield (per Gbp) at the given dose level is unchanged across the SOBP.

The cell survival curve for ⁴He ions is shown in figure 3(D). It is similar to the result obtained by using protons, except that beyond the SOBP, the surviving fraction is not 1.0 because of the non-zero dose in the fragmentation tail of ⁴He ions.

Figure 3(E) shows the beam compositions used to form the ⁴He SOBP. 39 groups of beam energies ranging from 78.0 to 116.0 MeV/u were used but the data from group #9 with 78.0 MeV/u and group #46 with 115.0 MeV/u are unseen owing to their relatively low contributions. In figure 3(F), the relative difference between the delivered biological dose after the optimization and the prescribed value varies from  −2.0% to 1.0% for the whole target region. However, except at the target boundaries, the relative difference is much lower and nearly within  ±0.3%, showing the high accuracy of the biological dose optimization.

To interpret the observations in the fragmentation tail, a 206 MeV/u ⁴He ion beam was selected as the representative to analyze the effects; specifically, the dose contributions from primary and different secondary particles were examined using Geant4. Our results, shown in figure 4, reveal that the tail's main constituents are species with Z  =  1, and more specifically, by secondary protons and deuterons (details not shown). The dose fragmentation tail is one of the major concerns for using ions heavier than protons because this tail can deliver unnecessary dose to surrounding normal tissues. Our simulation results also reveal that secondary ions have already started to contribute significantly to total dose even from the beam entrance to the Bragg peak as shown in figure 4.

Zoom In Zoom Out Reset image size

The Bragg peak of a ¹²C ion beam is so sharp that ripples may appear in the formed SOBP (Grevillot et al 2015). Therefore, a ridge filter, also called a ripple filter, is usually used to broaden the Bragg peak in clinical carbon ion therapy. A pseudo ripple filter made of ABS resin (density  =  1.04 g cm⁻³) has been modelled in Monte Carlo simulations of each pristine carbon ion beam to remove the ripples in the later optimized SOBP.

The target value was still set to 3.8 Gy (RBE) in the biological dose optimization procedure to expect a cell survival rate of 10% for H460 cells. In addition to the physical dose and biological dose, RBE, α and β, DSB yield, and cell survival rate of H460 cells from ¹²C ions with an SOBP ranging from 5 to 10 cm in water are shown in figures 5(A)–(D). The biological dose in the fragmentation tail of ¹²C beams is higher than for ⁴He beams (figure 3(A)) but the biological dose of ¹²C beams at the entrance is lower. No large data fluctuations (small ripples and spikes) were observed for RBE and α values in the fragmentation of ¹²C beams as were found with protons and ⁴He beams. This improved precision is possibly owing to the relatively high fluence of secondary particles.

Zoom In Zoom Out Reset image size

The high complexity of the radiation field generated by carbon ions has made the analysis of the biological effect of ¹²C ion beams more complicated than that for protons and ⁴He ions. As shown in figure 5(A), the RBE for ¹²C beams is much higher than for ⁴He and proton beams. Proximal to the SOBP, the RBE varies slowly from 1.4 to 1.9. Within the SOBP, the RBE increases up to 3.9. In the distal edge, the RBE increases up to 4.9 and then decreases. In the fragmentation tail, the RBE first decreases rapidly and then starts to increase slowly with depth, which is different from the continually increasing trend found for the ⁴He beams.

Equation (16) in the supplementary material can still be used to explain the biological observations. The RBE variation trend in figure 5(A) is the same as the trend for α and β in figure 5(B). The physical dose continues to drop slowly in the tail and the RBE can be qualitatively analyzed in the edge distal to the Bragg peak and the associated tail knowing that α and β are determined by the energy (or LET) of the ions. The characteristics of both of primary and secondary particles should be considered in analyzing the biological effects. In the distal edge, the decrease of primary beam energy (LET, α and β increase) and the decrease of dose can all result in the increase of the RBE. Later, all of the primary particles stop inside the medium, so only secondary particles can contribute to the dose. Without primary particles, there will not be more secondary ions generated and the fluence of existing secondary ions will continue to decrease. Some low energy recoils may stop as well, so the average energy of recoils may increase (similar to the beam hardening effect), resulting in a decrease of the LET and a decrease of α and β. Later, the secondary recoils continue to lose energy so the LET will start to increase, resulting in an increase of α and β and an increase of the RBE. The reported variation trend of α along the beam path using the RMF mechanism is consistent with the published work by Inaniwa et al (2010) using a modified MKM mechanism for human salivary gland tumor cells (Inaniwa et al 2010). The reported RBE values and spatial variation trend are consistent with the previously published data by Frese et al (2012) and Kamp et al (2015) (Frese et al 2012, Kamp et al 2015).

For the DSB induction data in figure 5(C), the normalized DSB yield (per Gy per Gbp) increases within the SOBP, but the absolute DSB yield (per Gbp) at a given dose level decreases with depth along the SOBP, which is also different from the constant absolute DSB yield observed in ⁴He beams. Although the absolute induced DSB decreases with depth along the SOBP, the frequency of intra-track exchange-type chromosome aberrations increases due to the very high LET at deeper depths. After DNA damage processing, the total lethal DNA lesions may be uniform across the SOBP, so the final cell survival rate is uniform within the SOBP.

The cell survival curve for ¹²C ions is shown in figure 5(D). In the fragmentation tail, the surviving fraction starts increasing from 0.86 to 0.96 at z  =  15 cm in water. Compared with proton and ⁴He beams, the cell survival in the fragmentation tail of ¹²C ions is lower because of the relatively higher biological dose mainly contributed by nuclear fragments.

Figure 5(E) shows the beam compositions used to form the ¹²C SOBP. 38 groups of beam energy from 115.0 to 224.0 MeV/u were used but the data from group #16 with 150.0 MeV/u and group #18 with 154.0 MeV/u are unseen owing to their relatively low contributions. In figure 5(F), the relative difference between the delivered biological dose after the optimization and the prescribed value varies from  −1.5% to 1.0% for the whole target region. However, except at the target boundaries, the relative difference is much lower and nearly within  ±0.3%, showing the high accuracy of the biological dose optimization.

Similar to the study for ⁴He ions, Geant4 Monte Carlo simulations were also performed to investigate the dose contributions from primary and secondary particles for a mono-energetic ¹²C ion beam (400 MeV/u), shown in figure 6. More nuclear interactions are involved between the ¹²C beam and the medium so the dose contributions from nuclear fragments are very high. Our results also reveal that the main nuclear fragments are protons and ⁴He ions (details not shown). The contributions from other secondary particles such as lithium, beryllium, and boron species are also considerable.

Zoom In Zoom Out Reset image size

The spatially varied RBE values for different ion beams have been reported above. Some concepts should be reiterated for appropriately describing and reporting the biological effects of particle therapy. One should be cautious in applying these quantities to avoid erroneously extending these results in ways that the data do not support.

First, the RBE for cell survival is defined as the dose ratio when the test radiation and reference radiation result in the same cell surviving fraction for the specified cell type. At different surviving fractions, the RBE values are dose and radiosensitivity dependent. Therefore, the basic requirement for reporting RBE and comparing RBE among different beam qualities is that the same cell surviving fraction criterion must be met (Hall and Giaccia 2012).

Second, the RBE is used to describe the relative level of the biological effect of the test radiation to the reference photons, rather than the absolute level. Therefore, a high RBE at a region does not indicate that the absolute biological effect therein, such as the cell kill, is high as well. For example, in the fragmentation tail of ⁴He ions, the RBE is high according to RMF model predictions, but the cell kill effect is low (as described in section 3.2). By contrast, the biological dose, i.e. the RBE-weighted dose to produce the same biological effect as the reference photons, is the appropriate radiobiological quantity to estimate the absolute biological effect. Moreover, an identical biological dose for the same cell line means having the same cell survival rate if the LQ cell survival model is applied.

Although the dose optimization strategy developed in the current study is very flexible, there are some limitations for our method of generating basic datasets. For example, the microdosimetric quantity y_F is not accurately calculated in Monte Carlo simulations. LET is simply used instead of y_F. This simplification is valid only when the energy of the charged particle is high enough so that it can pass through the target cell without changing the LET along the particle path through the cell (ICRU 1983, Lindborg and Waker 2017). When the particle range becomes comparable to the diameter of the target cell or cell nucleus, this simplification can generate large uncertainties because the particle's LET varies greatly along its path and it may stop inside the target. These effects are especially important for very low energy particles around the Bragg peak. These calculations of microdosimetric quantities will be refined in our future work to improve the accuracy of calculation results.

Another simplification in the current work is the assumption that lung cancer cells are distributed uniformly in the whole phantom. In practice, outside the target tumor are normal tissues. The biological response of cancer cells and normal cells to radiation can be quite different; therefore, the biological dose calculations for tumor cells and normal cells should be treated distinctly. For complete biological optimization, normal tissue tolerances and responses must also be understood and quantified in the appropriate biological context for effective modelling (Jones 1999). Thus, the tissue heterogeneities in real patients increase the complexity in making bio-IMIT plans. The experimental results of our ongoing studies in normal tissue response using the 3D culture technique will be reported in future. Even for a given cell type, the heterogeneity still exists in the radiobiological response. One practical way to study the impact of heterogeneity in radiosensitivity is to perform a sensitivity analysis of the modelling results over the expected distribution of radiosensitivity, which is currently under investigation by Carlson et al.

In addition, only 1D longitudinal data were used in the optimization process in the current study. To extend this work to more complex two- and three-dimensional heterogeneous patient geometry, more accurate dose calculation algorithms should be implemented in the optimization process. Extending our Python based repository to generate full biological dose IMIT optimizations for complex and heterogeneous patient geometries will be one of our next steps.

In addition to the mechanism-based RMF model demonstrated in this study, other biophysical models, e.g. MKM, LEM IV, and a recently developed α and β calculation formalism by Vassiliev et al (2017), can also be invoked in the optimization process. The data generation and optimization algorithms for these models will be developed in our future work.