First benchmarking of the BIANCA model for cell survival prediction in a clinical hadron therapy scenario

Up to now, the BIANCA model (version BIANCA II) has been applied to monochromatic beams of protons, C-ions and He ions, finding good agreement with experimental survival curves for V79 and AG01522 cells (e.g. Carante et al (2018)). In this work the approach was extended to O-ions, which are also of interest for hadron therapy (Tommasino et al 2015), and simulated survival curves were compared with V79 data taken from the literature. The value of the unrejoining fraction, f , was left unchanged with respect to the 0.08 value already used for protons, He-ions and C-ions. Thus, the CL yield was the only adjustable parameter. Figure 1 shows the results obtained for five O-ion beams, having LET values of 18, 46, 238, 276 and 754 keV μm⁻¹; the experimental data were taken from Stoll et al (1995). Since the data were derived from single experiments and thus do not have any associated uncertainty, a quantitative evaluation of the agreement is not feasible. However, the simulations are in line with the data for all considered LET values. Concerning the simulation outputs, here and in the following figures the associated errors, which were always within 5%, were not reported because they would be too small to be clearly visible.

Zoom In Zoom Out Reset image size

The mean number of CLs per unit track length, expressed as CL μm⁻¹, was found to increase with LET, reflecting the increasing energy-deposition clustering. However, the mean number of CLs per unit dose and DNA mass, expressed as CL · Gy⁻¹ · cell⁻¹, showed a peak at 238 keV μm⁻¹ and then decreased, reflecting the behaviour of the RBE (see below). This difference in the LET-dependence between CL μm⁻¹ and CL · Gy⁻¹ · cell⁻¹, which has been observed also for C-ions (Carante et al 2018), can be explained by considering the following relationship:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm CL}}{{\rm Gy}\cdot {\rm cell}}=6.25\cdot \frac{{\rm CL}}{\mu {\rm m}}\cdot \frac{V}{L}\nonumber \end{align} \tag{ 1 }$$

where L is the LET in keV μm⁻¹ and V is the cell nucleus volume in μm³. At very high LET, where we found that the yield of CL · μm⁻¹ increases in a sub-linear fashion, the yield of CL · Gy⁻¹ · cell⁻¹ starts decreasing due to the presence of the L⁻¹ term.

The yield of CLs found in this work for O-ions was very similar to that found in Carante et al (2018) for C-ions with the same LET, suggesting that these two ion types have a similar LET-dependence. More specifically, at 18 keV μm⁻¹ the CL yield (expressed as mean number of CLs per μm) was 0.009 for O-ions and .007 for C-ions; at 46 keV μm⁻¹ it was 0.031 for O-ions and 0.034 for C-ions; at 238 keV μm⁻¹ it was 0.389 for O-ions and 0.343 for C-ions; at 276 keV μm⁻¹ it was 0.384 for O-ions and 0.385 for C-ions; finally, at 754 keV μm⁻¹ it was 0.678 for O-ions and 0.657 for C-ions. The consequences of this issue will be discussed in section 3.2.

After adjusting the CL yield for survival curves obtained at different LET values of a given ion type, one can fit the LET-dependence of the CL yield for that ion type. Then, by taking the CL yield from the fit, it is possible to perform a full prediction of cell survival at any LET value. This is important in view of RBE calculations, for which the cell survival dose-response must be known for as many LET values as possible and for each involved particle.

Therefore, based on previous simulations on V79 cells exposed to protons, He-ions and C-ions of different LET values (Carante et al 2018), survival curves were predicted at many different LET values for these three ion types. More specifically, 12 survival curves were produced for protons (LET range: 2.5–30 keV μm⁻¹), 18 curves for He-ions (range: 5–90 keV μm⁻¹), and 29 curves for C-ions (range: 10–500 keV μm⁻¹).

Each of these (predicted) survival curves was then fitted by the 'usual' linear-quadratic expression S(D)  =  exp(−αD  −  βD²). The linear and quadratic coefficients derived from these fits allowed calculating RBE values for different ions at different LET values, according to the following expression (e.g. Mairani et al (2010)):

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm RBE}=2{{\beta }_{i}}\left[ -{{\alpha }_{X}}+\sqrt{\alpha _{X}^{2}-4{{\beta }_{X}}{\rm ln}S} \right]/2{{\beta }_{X}}\left[ -{{\alpha }_{i}}+\sqrt{\alpha _{i}^{2}-4{{\beta }_{i}}{\rm ln}S} \right]\nonumber \end{align} \tag{ 2 }$$

where α_X and β_X are the photon coefficients, α_i and β_i are the ion coefficients (at a given LET), and S is the considered survival level.

In figure 2, panel (a), the lines represent the RBE at 10% survival calculated for protons and He-ions. The RBE values calculated for protons were compared with experimental data taken from Belli et al (1998) (LET values: 7.7, 11.0, 20.0 and 30.5 keV μm⁻¹) and Folkard et al (1996) (LET values: 10.1, 17.8 and 27.6 keV μm⁻¹), whereas those for He-ions were compared with ³He data taken from Furusawa et al (2000) (20 points in the LET range 18.6–90.8 keV μm⁻¹), as well as a ⁷Li point at 60 keV μm⁻¹ taken from Pathak et al (2007). The model calculations were in good agreement with the data, showing that BIANCA can be used for RBE predictions for V79 cells exposed to protons or He-ions. In particular for protons, for which experimental errors are available, all the simulation outcomes were within 3σ with respect to the data. Concerning the discrepancy between simulations and data, only for the point at 7.7 keV μm⁻¹ the discrepancy was 19.8%, whereas for the other points it was less than 15%. The fact that the ⁷Li point is higher than the He value at the same LET can be explained by considering that, in the ⁷Li experiment, ⁶⁰Co gamma-rays were used as a reference radiation instead of x-rays. This suggests that the BIANCA predictions for He-ions may also be applied to ions with Z  =  3, although further comparisons are desirable.

Zoom In Zoom Out Reset image size

Concerning heavier ions, the 10% RBE predictions for carbon were compared with C-ion data taken from Furusawa et al (2000) (24 points in the LET range 22.5–502.0 keV μm⁻¹), as well as five O-ion points taken from Stoll et al (1995) (LET range: 18–754 keV μm⁻¹), two boron points taken from Cox et al (1977) (LET: 110 and 200 keV μm⁻¹), and one nitrogen point also taken from Cox et al (1977) (LET: 470 keV μm⁻¹). As shown in figure 2, panel (b), also for these ion types the model calculations were in good agreement with the data. In particular for B-, O- and N-ions, for which the errors were available, the (reduced) chi-square was 1.5. Except for 46 keV μm⁻¹ O-ions, for which the discrepancy between simulations and data was 30%, for the other points it was  ⩽15%. This shows that BIANCA can be used for RBE predictions for V79 cells exposed to C-ions, and suggests that the C-ion predictions may also be applied to ions with Z  =  5, Z  =  7 and Z  =  8.

The comparisons performed at 10% survival were repeated at 50% survival. The results obtained are reported in figure 3 (panel (a) for light ions and panel (b) for heavy ions). Concerning light ions, the proton results were compared with data taken from Folkard et al (1996) and Belli et al (1998), whereas the He-ion results were compared with ³He data taken from Furusawa et al (2000), as well as a ⁷Li point taken from Pathak et al (2007). Except for an underestimation of proton effectiveness around 30 keV μm⁻¹, which will be discussed below, the light-ion calculations were in good agreement with the data. In particular for the first five proton points the (reduced) chi-square was 1.6, and the maximum discrepancy between simulations and data was 13%. This showed that BIANCA can predict the 50% survival RBE for V79 cells exposed to protons or He-ions, and possibly ions with Z  =  3.

Zoom In Zoom Out Reset image size

An underestimation of the effectiveness of higher LET protons was expected, since for these particles the model tends to overestimate the surviving fraction at low doses; as discussed previously (Carante et al 2018), this is possibly related to the choice of the function that describes the chromosome fragment rejoining probability. This also explains why this RBE underestimation was not found for 10% survival RBE, which involves higher doses. In any case, it should be taken into account that these LET values are very unlikely to play a relevant role in proton therapy, not even in the SOBP distal region. For instance, for the proton beam used at INFN-LNS in Catania (Italy) to treat ocular tumours, the (average) LET in the most distal SOBP region, just before the dose fall-off, is 21.7 keV μm⁻¹ (Chaudhary et al 2014) Even lower values, below 15 keV μm⁻¹, are reported in Giovannini et al (2016) for two examples of cranial irradiation.

Concerning heavy ions, the C-ion predictions for RBE at 50% survival were compared with C-ion data taken from Furusawa et al (2000), as well as five O-ion points taken from Stoll et al (1995), two Boron points taken from Cox et al (1977), and one N point taken from Cox et al (1977). For B-, O- and N-ions, for which the errors were available, the (reduced) chi-square was 1.8, and all simulation outcomes were within 3 σ with respect to the data. Concerning the discrepancy between simulations and data, except for 18 keV μm⁻¹ O-ions, for which the discrepancy was 27%, for the other points it was  ⩽19%. The good agreement between simulations and data showed that BIANCA can be used for predicting the 50% survival RBE for V79 cells exposed to C-ions, and suggested that the C-ion predictions may also be applied to ions with Z  =  5 and Z  =  8, as well as Z  =  7 that is expected to have an intermediate behaviour between Z  =  6 and Z  =  8.

According to the model features, in principle the CL yield for a given radiation quality (i.e. particle type and energy) should be adjusted separately for each considered cell line, since the CL yield depends not only on radiation quality, but also on the target cell features. However, recently we have developed a formula to fully predict the ion-survival of the cell line of interest based on the ion-survival of a reference cell line, as well as the photon response of both (Carante et al 2018). More specifically, the CL yield to predict the survival of the cell line of interest following exposure to a given radiation quality (i.e. a given ion type and energy) can be derived as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \frac{{\rm CL}}{\mu {\rm m}}={{\left( \frac{{\rm CL}}{\mu {\rm m}} \right)}_{ref}}\cdot \left[ {\left( \frac{{\rm CL}}{{\rm Gy}\cdot {\rm cell}} \right)}\bigg /{{{\left( \frac{{\rm CL}}{{\rm Gy}\cdot {\rm cell}} \right)}_{ref}}} \right]\cdot \frac{{{V}_{ref}}}{V}.\nonumber \end{align} \tag{ 3 }$$

Here, (CL/μm)_ref is the CL yield used for the reference cell line exposed to the same radiation quality, (CL · Gy⁻¹ · cell⁻¹) and (CL · Gy⁻¹ · cell⁻¹)_ref are the CL yields used to simulate photon exposure for the cell line of interest and the reference cell line, respectively, and finally V_ref and V is the nucleus volume used for the reference cell line and the cell line of interest, respectively. In a previous work (Carante et al 2018), this approach has been successfully applied to AG01522 and U87 cells exposed to proton beams of different LET.

In the present work we adopted the same approach for CHO cells irradiated at HIT by two opposing fields of protons or C-ions mimicking a typical patient treatment scenario, as described in Elsässer et al (2010). In that study the cells were grown on several plastic slices that, immediately before irradiation, were inserted into an acrylic glass container (5  ×  10  ×  16 cm³), allowing for a 5 mm spatial resolution along the beam direction. The cells were then irradiated by two opposing beams with energy in the range 90–120 MeV/u (protons) or 175–230 MeV/u (carbon), delivering 1.5 Gy in the entrance channel (for both protons and carbon) and 5.3 Gy (protons) or 3.9 Gy (carbon) in the center of the target region. The irradiations were planned by the treatment planning software TRiP98 (Kramer and Scholz 2000) and optimized by LEM to achieve a homogeneous cell survival level in the target, simulating a 4 cm target located between 6 and 10 cm water-equivalent depth. Further details can be found in Elsässer et al (2010), which also reports the cell survival predictions performed by the LEM model.

As a first step, to fix the model parameters for photons, the photon CHO survival curve reported in Elsässer et al (2010) (α  =  0.105 Gy⁻¹ and β  =  0.025 Gy⁻²) was reproduced by BIANCA, using a CL yield of 1.2 CL · Gy⁻¹ · cell⁻¹ and an un-rejoining parameter f   =  0.08. The 1.2 value was inserted in equation (3) to derive the CL yields for simulating the survival of CHO cells exposed to different ion types and energies; the same nucleus volume was used for CHO and V79 cells. Afterwards, these CL yields were used as code input values to produce CHO survival curves for many different (monochromatic) beams (12 proton beams having LET in the range 2.5–30.0 keV μm⁻¹, 18 He-ion beams with LET in the range 5–90 keV μm⁻¹, and 29 C-ion beams with LET in the range 10–500 keV μm⁻¹). Each of these survival curves was then fitted by a function of the form S(D)  =  exp(−αD  −  βD²), which allowed constructing a table consisting of the linear and quadratic coefficients for different ion types and LET values, analogous to the V79 table described in section 3.2. As mentioned in section 3.2, as a practical solution, for Z  =  3 we used the coefficients derived for He ions, whereas for Z  =  4 or Z  =  5 or Z  =  7, 8 we used the coefficients derived for C-ions.

Afterwards, this table was read by the FLUKA radiation transport code. The FLUKA development version 2018.1 was used. The settings corresponding to the SDUM option PRECISIOn in the DEFAULTS card were implemented. For carbon ions, the external event generator relativistic quantum molecular dynamic (rQMD-2.4) (Sorge et al 1989) was linked to assure accurate modelling of nucleus-nucleus collisions above 125 MeV/u. Moreover, coalescence and evaporation of heavy ions were activated by means of two PHYSICS cards. The nominal beam energies available at HIT (e.g. Parodi et al (2012)) for carbon ions in the range 175–230 MeV/u, and for protons in the range 90–120 MeV/u (as in the experiment reported in Elsässer et al (2010)) were used in FLUKA to model the irradiation field. The weights for each single Bragg peak contributing to the 4 cm extended SOBP were calculated by means of an in-house written routine. Two opposing fields were used like in the reference experiment. Simulations were performed on a water target, and the necessary information (particle type and energy, and absorbed dose) were calculated by FLUKA on a voxel-by-voxel basis. More specifically, whenever according to FLUKA a certain amount of energy was deposited in a target voxel by a given particle type of given energy (and thus given LET), FLUKA read the corresponding alpha and beta coefficients from the table produced by BIANCA. Afterwards, to account for the effects of mixed fields, the mean value of alpha and the mean value of beta in that voxel were calculated according to the theory of dual radiation action (Zaider and Rossi 1980) as described in Mairani et al (2010), i.e.:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha =\frac{\mathop{\sum }_{i=1}^{n}{{\alpha }_{i}}{{D}_{i}}}{\mathop{\sum }_{i=1}^{n}{{D}_{i}}}\nonumber \end{align} \tag{ 4 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \sqrt{\beta }=\frac{\mathop{\sum }_{i=1}^{n}\sqrt{{{\beta }_{i}}}{{D}_{i}}}{\mathop{\sum }_{i=1}^{n}{{D}_{i}}}.\nonumber \end{align} \tag{ 5 }$$

Here, D_i is the absorbed dose due to the ith particle according to FLUKA, α_i and β_i are the corresponding radiobiological coefficients provided by BIANCA, and α and β are their mean values. After calculating α and β as described above, the corresponding survival level, the RBE-weighted dose and the RBE were calculated as follows:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle -{\rm ln}S=\alpha D+\beta {{D}^{2}}\nonumber \end{align} \tag{ 6 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {{D}_{{\rm RBE}}}=\left[ \sqrt{\alpha _{X}^{2}-4{{\beta }_{X}}{\rm ln}S} \right]/2{{\beta }_{X}}\nonumber \end{align} \tag{ 7 }$$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle {\rm RBE}=\frac{{{D}_{{\rm RBE}}}}{D}\nonumber \end{align} \tag{ 8 }$$

where α_X and β_X are the photon parameters provided by BIANCA.

The results obtained for carbon are reported in figure 4 (panel (a)), where the solid line represents the prediction performed by BIANCA interfaced to FLUKA, the points are the experimental data taken from Elsässer et al (2010) by using a software available online called WebPlotDigitizer (https://automeris.io/WebPlotDigitizer), and the dashed line is the prediction performed by LEM (version LEM IV) integrated into TRiP98, which was also taken from Elsässer et al (2010) by means of WebPlotDigitizer. The predictions of the two models (BIANCA and LEM) are rather similar and are both in good agreement with the data. In particular, the (reduced) chi-square for the BIANCA simulations was 2.2, and the discrepancy between BIANCA and the data was  ⩽17%.

Zoom In Zoom Out Reset image size

This suggests that also BIANCA, following interface to a radiation transport code like FLUKA, can be used to predict cell survival along C-ion hadron therapy beams in a clinical scenario. Panel b reports the same kind of results for protons. Although BIANCA tends to predict a lower survival—and thus a higher biological effectiveness—with respect to LEM, especially in the entrance channel, both models are consistent with the data.

Overall, the good agreement with the CHO data suggests that the method described in this section can be applied, in principle, to any cell line. This would allow performing ad hoc RBE predictions specific for the cell type of interest.