Quantification of biological range uncertainties in patients treated at the Krakow proton therapy centre

Variable relative biological effectiveness (vRBE) in proton therapy might significantly modify the prediction of RBE-weighted dose delivered to a patient during proton therapy. In this study we will present a method to quantify the biological range extension of the proton beam, which results from the application of vRBE approach in RBE-weighted dose calculation. The treatment plans of 95 patients (brain and skull base patients) were used for RBE-weighted dose calculation with constant and the McNamara RBE model. For this purpose the Monte Carlo tool FRED was used. The RBE-weighted dose distributions were analysed using indices from dose-volume histograms. We used the volumes receiving at least 95% of the prescribed dose (V95) to estimate the biological range extension resulting from vRBE approach. The vRBE model shows higher median value of relative deposited dose and D95 in the planning target volume by around 1% for brain patients and 4% for skull base patients. The maximum doses in organs at risk calculated with vRBE was up to 14 Gy above dose limit. The mean biological range extension was greater than 0.4 cm. Our method of estimation of biological range extension is insensitive for dose inhomogeneities and can be easily used for different proton plans with intensity-modulated proton therapy (IMPT) optimization. Using volumes instead of dose profiles, which is the common method, is more universal. However it was tested only for IMPT plans on fields arranged around the tumor area. Adopting a vRBE model results in an increase in dose and an extension of the beam range, which is especially disadvantageous in cancers close to organs at risk. Our results support the need to re-optimization of proton treatment plans when considering vRBE.

doses to healthy tissues compared to photon radiation therapy [2]. Proton range for dose distribution is usually defined as the depth where the dose drops to 80% or 90% of planned dose according to Paganetti [3]. Proton radiation therapy gives the possibility to efficiently treat malignant tumors, administering lower doses to healthy tissues compared to photon radiation therapy [2]. Protons also exhibit an increased relative biological effectiveness (RBE, defined as the ratio of photon to proton physical dose needed to achieve the same biological effect), which makes them more effective than photons [4]. In clinical practice, safety margins are added to the clinical target volume (CTV), thus generating the planning target volume (PTV) to ensure that the entire potentially malignant tissue is irradiated by taking into account the range uncertainty and tumor motion. The conformity of PT dose distributions, that is essential for patient safety, is obtained by the combination of proton beam properties, application of multiple treatment fields and intensity modulation technique. Nevertheless, robustness of PT treatment plans is challenged by proton beam range uncertainty in the patient [5].
An emerging issue of PT treatment planning is the variation of RBE values across the target volume; currently, it is assumed that RBE has a constant value of 1.1 in clinical practice. In fact, radiobiological experiments show that the RBE values at the end of the proton range (i.e. Bragg peak distal fall-off ) increases up to 2 [6], according to the review in Ilicic et al. [7]. The clinical evidence for an enhanced proton RBE at the distal edges has been reported in the very recent summary report of an EPTN workshop [8]. This variation of RBE along the proton path is an additional source of uncertainty in PT treatment planning. In order to obtain more precise information on dose distribution in patients, advanced dose calculation methods based on analytical or Monte Carlo (MC) algorithms are used for treatment planning (e.g. Eclipse [9], RayStation [10]) or to support treatment planning, quality assurance and research (e.g., TOPAS [11], GATE [12], Fluka [13], FRED [14], gPMC [15], MCsquare [16]). These tools can calculate physical and RBEweighted dose (D RBE ) distributions. The latter is typically obtained, for protons, with a parametric model using the dose-averaged linear energy transfer (LET d ) [4] and the information on radiosensitivity of different tissue types characterized by the α/β ratio [17,18]. Identification of the increased LET d regions (which also increase the RBE values) during treatment planning might be used to minimize the LET d and RBE outside the PTV, especially in organs at risk (OARs) [19,20].
In this paper we investigate, quantify and compare biological uncertainties of clinical PT treatment plans with constant RBE (cRBE) and variable RBE (vRBE), focusing on the range uncertainty and evaluating clinical parameters of dose-volume histograms (DVHs). We performed MC simulations of clinical treatment plans using FRED (Fast paRticle thErapy Dose evaluator) [14]. We used simulation data to propose a method of quantification of biological range extension (or biologically effective range, as defined in Grün et al. [21]). Further, we compared biological uncertainties of the treatment plan with uncertainties resulting from its robustness to patient positioning and computed tomography (CT) calibration.

Patients database description
In our study we used CT images (Siemens Somatom Definition AS), as well as treatment plans and dose distributions calculated in the treatment planning system (TPS) Eclipse 13.6 for 95 patients with brain and skull base tumors, treated with protons from November 2016 to September 2018 at the Cyclotron Centre Bronowice (CCB) in Krakow. The patients underwent intensitymodulated proton therapy (IMPT) with fraction and total RBE-weighted doses (accounting for RBE = 1.1) ranging from 1.8 to 2.0 Gy(RBE) and 36 to 74 Gy(RBE), respectively. The treatments were prepared with multi-field optimization (MFO) and carried out in 1-3 stages, where stage 2 and 3 were the boost plans. The patients were divided according to the tumor type into two groups, i.e., brain patients and skull base patients, which is presented in Table 1. Information on the PTV volume, prescribed dose and diagnosis for each patient can be found in Table S1 in Additional file 1. Aiming at the coherence of the patient database and most clinically representative results, we excluded 20 patients from the database with tumor localization other than brain or head and neck (H&N), pediatric cases, replanned patients and treatment plans optimized with single field uniform dose (SFUD) approach.

LET, RBE and dose calculations in FRED MC
In this work we performed MC calculations of D RBE and LET d distributions delivered by treatment plans prepared in TPS. The MC calculations were made in a patient anatomy defined by its CT images, resampled from the original grid size of 0.67 × 0.67 × 1.2 mm 3 to 1.5 × 1.5 × 1.5 mm 3 . The calculations were performed with the MC tool FRED. FRED allows for particle tracking on graphical processing units (GPUs), which shortens the calculation time for the whole treatment plan to several minutes, which is about 1000 times faster than general purpose MC tools executed on central processing units (CPUs) [14]. The beam model used for patient treatment at the Krakow PT center was implemented in FRED and validated against commissioning data, measurements in homogeneous and heterogeneous media and standard clinical calculations [22]. FRED includes several radiobiological models, both parametric and table based (e.g. Carabe [23], McNamara [24], LEM1 [25]), allowing storing RBE distributions voxel by voxel.
In this paper we computed D RBE distributions using the parametric vRBE model by McNamara since it is based on the most recent and comprehensive set of radiobiological data and accounts for the LET d distribution in considered volume [24]. The α/β ratios were chosen based on the tumor type and overview presented in van Leeuwen et al. [26]. We used α/β ratios of 6 Gy and 4 Gy for brain patients and skull base patients, respectively, while D RBE to normal tissue was computed using α/β of 2 Gy. RBE values were extracted for each patient individually for the voxels with the RBE-weighted dose values exceeding 5% of the prescribed dose. We studied the differences in RBE values in the PTV and OARs, where elevated RBE can cause post-treatment side effects in patients. The LET d values were calculated in each voxel from all primary and secondary protons as a sum of deposited energy of all events weighted by the particle stopping power as in Polster et al. [27]: where dE-deposited energy, dx-path length, ρ-density of the medium and ΣdE-the sum runs on all events of energy deposition.

Evaluation of RBE-weighted dose
To compare the D RBE distributions computed with cRBE and vRBE, DVHs and the DVH parameters were computed for each patient and each dose distribution. The mean and maximum doses in PTV as well as in brainstem and chiasm were evaluated. The DVHs were generated using the dicompyler-core library [28]. All data analysis was performed using in-house Python scripts that were validated against Eclipse TPS.

Biological range extension
Most of the patient plans consist of 3-4 treatment fields from different directions, which were arranged around the target, therefore the range can be extended in different directions. For evaluation of biological range extension we used an approximation, by assuming that the volume covered with at least 95% of the prescribed dose (V95) computed with cRBE and vRBE has a spherical shape. We thus defined biological range extension as the difference between radiuses: where R vRBE -radius of sphere of V95 volume computed with vRBE, R cRBE -radius of sphere of V95 volume computed with cRBE. The schematic representation of this method is presented in Fig. 1.

Comparison of D RBE uncertainty and treatment plan robustness
We evaluated physical and biological uncertainty of treatment plans by comparing the DVHs of D RBE distributions for PTV and OAR structures computed with: cRBE versus vRBE, and cRBE versus treatment plan robustness uncertainty bands. We have chosen the most critical structures (i.e. OARs with the highest dose) for tumor types selected for this study. We also analyzed the dose in other OAR including the spinal cord, which in majority of the cases is substantially below the dose constraint for this organ regardless of the RBE modeling approach. The robustness analysis was performed in FRED following the clinical procedure implemented at CCB: each clinical treatment plan was recalculated with cRBE for 12 cases: modifying CT image Hounsfield units (HU) values by ± 3.5% and translating the treatment plan isocenter of ± 2 mm in X, Y and Z directions. The DVHs with error bands were generated for PTV and OARs. Generally, the PTV was created by expanding the CTV with a safety margin of 3 mm, taking into account patient anatomy variations and patient positioning. In addition, CTV was extended to a complementary technical PTV structure (PTV tech ) that includes the CT calibration curve uncertainty and is equivalent to a safety margin of 3.5% of the maximal beam range plus 1 mm for a given treatment field. The treatment dose was eventually prescribed to the PTV tech that is the envelope of the maximum safety margin of PTV or PTV tech . The aim was to perform a comparative study of the dose variations originating from: the vRBE modelling against the CT translation, which mimics the patient mispositioning and the uncertainty of HU values due to CT scanner calibration, both performed with cRBE.
To quantify the differences between DVHs for cRBE and vRBE models the DVH indices such as D max (the maximum local dose in he considered structure), D mean (the mean dose in the considered structure), D05 (the maximum dose covering 5% of the considered structure) and D95 (the maximum dose covering 95% of the considered structure). The D95 and D05 parameters were analyzed for CTV and OARs, respectively, according to the formulas:

Time performance
The time performance of FRED calculations was discussed before [14,22]. In this study we simulated 10 4 primary protons per pencil beam for each treatment plan. The mean(std) simulation time for brain patients was 2.75(1.35) min and for skull base patients 1.74(0.94) min, which makes FRED calculations sufficiently fast for its application in the clinical routine to support treatment planning and quality assurance.

Statistical analysis
No a priori assumption on the distribution of the dose values was done. The 95% confidence intervals (2.5-97.5 percentile, CI 95% ) were calculated with the bootstrap method (100,000 iterations) around the median value for relative D RBE values and biological range extension. Estimation of statistical difference between doses and changes in proton range using cRBE and vRBE was performed using the two-sided Wilcoxon signed-rank test [29].

Validation of LET d and D RBE calculations in FRED
The accuracy of D RBE and LET d distribution calculations in FRED were validated against TOPAS simulations from McNamara et al. [30], which is presented on Fig. 2. The FRED simulation of a dose cube optimized with Eclipse TPS of modulation width from 150 to 250 mm was performed in a virtual water phantom with statistics of 10 4 primary protons per pencil beam.
The maximum difference in D RBE computed with the vRBE model in FRED and TOPAS (adapted from [30]) was 3%. The difference between maximum LET d values was up to 13%, which does not substantially influence the RBE values in patients.

Comparison of D RBE in PTV
Two patient groups: brain patients (Patient 1-50) and skull base patients (Patient 51-95) were analysed separately. Patients were ranked according to increasing mean vRBE-weighted dose in the PTV. The comparison of relative mean dose in PTV to the prescribed dose are presented on the top panels of Figs. 3 and 4, with mean, median, minimum and maximum doses over all patients provided in the legend. For brain patients, the CI 95% of D mean /D p median was: The p value from Wilcoxon test performed between D RBE computed in FRED with cRBE and vRBE models was p < 0.05 (p = 1.8 × 10 -7 and p = 5.2 × 10 -9 for brain and skull base patients, respectively), which means that the mean deposited dose in the PTV relative to prescribed doses in PTV (D mean /D p ) are significantly different.
The values of D95 and D05 parameters computed for both groups of patients are presented in Figs. 3

Comparison of D RBE in OARs
The maximum dose (D max ) in brainstem and chiasm was analyzed for both groups of patients and presented on the fourth and the fifth panel of Figs. 3 and 4, respectively. The D max in OARs was compared to the clinical constraints, which is 54 Gy(RBE) for brainstem core, 64 Gy(RBE) for brainstem surface and 54 Gy(RBE) for chiasm. The DVH analysis was performed on the whole brainstem structure. The median dose values computed with TPS and FRED for cRBE model are close to the dose constraints, while for the vRBE model the doses were much higher and often above dose constraints (up to 14 Gy above the limit), especially for skull base patients (73% and 100% of the patients for chiasm and brainstem, respectively).

RBE distribution for vRBE model
In the PTV for brain patients (α/β = 6 Gy), the mean of the RBE values was 1.11, ranging from 1.07 to 1.31, whereas for skull base patients (α/β = 4 Gy) the mean RBE was 1.15 ranging from 1.09 to 1.49. In the brainstem the RBE values reached up to 2.71 and 2.36, while in chiasm the values were up to 2.16 and 2.19 for brain and skull base patients, respectively. The escalation of RBE values after the Bragg peak around 2.5 was also reported recently in Missiaggia et al. [31]. Note that depending on the local dose in a voxel, the high LET d or RBE value may have severe or negligible impact on the OARs side effects. However, we found that for the majority of investigated patient cases, the RBE-weighted dose with vRBE has a median value of over 50 Gy(RBE), so even the slightly increased RBE can potentially impact the clinically relevant biological dose in brainstem or chiasm (with dose constraints for brainstem core and chiasm of 54 Gy(RBE)), and eventually have severe impact on these organs.
The beam arrangement with the D RBE calculated with FRED cRBE model is presented on Fig. 5. Figure 6 shows the D RBE , LET d and RBE distributions with delineated PTV and OARs for an exemplary brain patient (panels (a)-(f ) on Fig. 6) and skull base patient (panels (g)-(l) on Fig. 6). The dose difference between FRED cRBE and vRBE calculations is visible for both patient cases around the PTV contour. The increase of D RBE computed using vRBE is the consequence of an increased RBE, especially at the distal edge of the dose fields.

Mean biological range extension
The V95 volumes for brain patients computed with cRBE ranged from 59.75 to 1066.97 cm 3 Fig. 6. The mean biological range extension R ext (formula(2)) for brain patients was R ext (brain) = 0.44(0.14) cm and for skull base patients R ext (skull base) = 0.45(0.11) cm.
Results for all patients correlated with the volume of the V95 are presented on Figs. 7 and 8. The Pearson correlation coefficient between the V95 (calculated with vRBE) and biological range extension was 0.25 and 0.11 for brain and skull base patients, respectively. The results show that, at least for the patients considered in this study, the biological range extension does not significantly vary with the tumor localization and type.

Comparison of D RBE uncertainty and treatment plan robustness
The results were analyzed and presented as DVHs of CTV (magenta curves), chiasm (green curves) and brainstem (blue curves) volume in Fig. 9. The shaded regions around the DVH computed for the dose calculation performed with cRBE model show the variation of the dose distribution resulting from the robustness analysis. The mean difference for CTV from D95 parameter (calculated by formulas (3) and (4)) were: Gy(RBE) for brain and skull base patients, respectively. The mean differences for OARs form D05 parameter (calculated by formulas (5) and (6)) were: iii. OAR(brainstem for brain and skull base patients, respectively. Results for the selected brain and skull base patients show that incorporating the vRBE model can modify the shape of the DVH (due the dose escalation) compared to the cRBE model (as also reported in Tommasino et al. [32]), to a level comparable or even higher than changes caused by patient misalignment or CT number uncertainty (shaded areas), especially in high dose regions.

Discussion
The brain and skull base tumors are particularly challenging targets, because they are often localized very close to the critical structures (e.g. brainstem or chiasm, see   the application of vRBE model causes an increase of the dose in PTV and OARs (see Fig. 3 and 4), which in turn enlarges the region covered by the high D RBE . For brain and H&N patients Yepes et al. [33] showed RBE values greater than 1.4 and 1.3, respectively using α/β = 2 Gy, which is a similar result to one obtained in this study. We have shown that for glioma tumor tissues with high α/β ratios of 6 Gy, the McNamara RBE model on average does not predict much higher RBE values than 1.1 (see Fig. 3, first panel).
By analyzing volumes receiving at least 95% of the prescribed dose using cRBE and vRBE we estimated the average biological range extension of over 0.4 cm. This result is similar to the predictions of biological range with vRBE obtained by comparing the edges of high dose regions (i.e. range at 80% of the prescribed dose) in Giovannini et al. [34] and Grün et al. [21]. We also showed that incorporating the vRBE model, the range of D RBE distributions in brain and skull base patients was extended on a similar level in both groups and was not correlated with the volume of the target exposed to high dose (more than 95% of prescribed dose). Minimizing the extension of the proton range and not exceeding dose limits in normal tissue is expected to be the factor limiting brain toxicity in H&N patients [35].
Monte Carlo simulation methods provide a reliable prediction of dose deposited in patients [22,36]. Due to the wide availability of inexpensive GPU cards providing high computational performance, the MC calculation algorithms, including implementation of vRBE models can support clinical routine, helping to identify the regions potentially exposed to increased LET d and/or D RBE . As our results show, the MC simulations indicate that real dose distribution can be more inhomogeneous than predicted by the analytical TPS, which results in lower mean dose in PTV using the same RBE model, which was also reported in Tseng et al. [37] and Ytre-Hauge et al. [38]. The RBE value depends on several variables and varies from patient to patient. It was shown that median RBE values may vary between patients by up to 15% [38].
We suggest that MC calculations of LET d and variable RBE should support the treatment plan evaluation. Since including vRBE modeling in clinical treatment planning is still controversial, we believe that treatment planning optimization in the near future should include simultaneous optimization of LET d and biological dose with cRBE in critical structures [39]. This will allow to reduce vRBE in OAR without substantial modification of biologically weighted dose computed with cRBE. Our results show that the region of the high dose computed with vRBE may expand significantly in comparison to cRBE. We should particularly consider these cases, when the PTV area is located near to the critical structures, which, exposed to high doses, can be affected by side effects after proton therapy. The short time needed for dose calculation in FRED enables vRBE evaluation in parallel to the clinical treatment planning workflow, while not affecting the treatment plan preparation time.

Conclusions
As presented in this work, incorporating the vRBE approach would modify the D RBE distribution, leading to an additional source of D RBE and range uncertainty in treatment planning. We suggest that the treatment planning procedure should account for the uncertainty of the D RBE , particularly in high dose regions, predicted with the vRBE models. The biological range extension due to vRBE can be even larger than planned margins in PTV. Our approach of estimating the biological range extension is a suitable method for IMPT plans because, differently to the analysis of transverse profiles, is insensitive to dose inhomogeneities and radiation field directions.