Biophysical modeling and experimental validation of relative biological effectiveness (RBE) for 4He ion beam therapy

Background Helium (4He) ion beam therapy provides favorable biophysical characteristics compared to currently administered particle therapies, i.e., reduced lateral scattering and enhanced biological damage to deep-seated tumors like heavier ions, while simultaneously lessened particle fragmentation in distal healthy tissues as observed with lighter protons. Despite these biophysical advantages, raster-scanning 4He ion therapy remains poorly explored e.g., clinical translational is hampered by the lack of reliable and robust estimation of physical and radiobiological uncertainties. Therefore, prior to the upcoming 4He ion therapy program at the Heidelberg Ion-beam Therapy Center (HIT), we aimed to characterize the biophysical phenomena of 4He ion beams and various aspects of the associated models for clinical integration. Methods Characterization of biological effect for 4He ion beams was performed in both homogenous and patient-like treatment scenarios using innovative models for estimation of relative biological effectiveness (RBE) in silico and their experimental validation using clonogenic cell survival as the gold-standard surrogate. Towards translation of RBE models in patients, the first GPU-based treatment planning system (non-commercial) for raster-scanning 4He ion beams was devised in-house (FRoG). Results Our data indicate clinically relevant uncertainty of ±5–10% across different model simulations, highlighting their distinct biological and computational methodologies. The in vitro surrogate for highly radio-resistant tissues presented large RBE variability and uncertainty within the clinical dose range. Conclusions Existing phenomenological and mechanistic/biophysical models were successfully integrated and validated in both Monte Carlo and GPU-accelerated analytical platforms against in vitro experiments, and tested using pristine peaks and clinical fields in highly radio-resistant tissues where models exhibit the greatest RBE uncertainty. Together, these efforts mark an important step towards clinical translation of raster-scanning 4He ion beam therapy to the clinic. Electronic supplementary material The online version of this article (10.1186/s13014-019-1295-z) contains supplementary material, which is available to authorized users.


Background
With nearly 150,000 patients treated globally to date, particle therapy has revolutionized cancer therapy by offering enhanced precision and radiobiological properties over the conventional photons [1]. At the Heidelberg Ion-Beam Therapy Center (HIT), proton ( 1 H) and carbon ( 12 C) ion beams, the leading modalities in hadrontherapy, are applied clinically, with two additional particle species available for experimentation: oxygen ( 16 O) and helium ( 4 He) ion beams. Interest in medical applications using helium ions began during the clinical trials at Lawrence Berkeley Laboratory (LBL) between the years of 1977 and 1993, with over 2000 patients successfully treated [2]. Since the program's end, 4 He ion beams remain clinically unexploited.
It is well known that, experimentally, heavier ions exhibit greater biological damage and consequently, the biophysical properties of 4 He are intermediate of the two clinically administered particle beams. That being said, application of helium ions provides a distinct clinical advantage, i.e. favorable dose distributions with attributes such as a sharper Bragg peak and lateral penumbra (reduced range straggling and scattering) compared to protons, and similar potential for tumor control with a substantially reduced fragmentation tail compared to carbons ions [3,4]. With these characteristics, helium ions have been proposed as an ideal treatment option for radio-resistant diseases and delicate patient cases e.g. meningioma and pediatrics [5,6].
Next year, HIT will launch the first European clinical program using therapeutic 4 He ion beams, which marks the world's first clinical application of raster-scanning 4 He ion therapy. Over the past decade, substantial efforts have been made at HIT to characterize 4 He ion beams via measurement and FLUKA Monte Carlo (MC) simulation [7,8] both dosimetrically, i.e. in terms of depth and lateral dose distributions with single pencil beam (PB) and spread-out Bragg peak (SOBP) plans, as well as nuclear fragmentation [9][10][11][12]. In addition, classification of the beam's biological effects is in progress, studying both in silico [5] and clonogenic cell survival in clinically-relevant conditions [13][14][15]. Presently, there is no commercial treatment planning system (TPS) available for 4 He ion beams; however, research-based tools were recently introduced or updated to allow planning with 4 He ion beams [10,14,16].
Relative to the clinical standard photons and protons, 4 He ion beams exhibit, in certain cases, more advantageous biological dose distributions with a higher linear energy transfer (LET) [17] in the tumor, resulting in superior relative biological effectiveness (RBE) in the target compared to the entrance channel, a valuable attribute for treatment of deep-seated radio-resistant tumors. To anticipate variability of tissue-specific radio-sensitivity in the clinic, the TPS predictions of physical dose will be coupled with a biophysical (RBE) model for calculation of an effective dose.
In contrast with proton RBE with nearly 300 experimental in vitro measurements, data for helium is relatively scarce (~1/3 as large), leading to larger uncertainties in helium RBE. As for in vivo investigation of 4 He ion beams, few publications examine evidence of enhanced tumor control compared to conventional techniques, most of which originate from the LBL trials from prior decades, yet only a fraction of these works relate findings to RBE [18,19]. In preparation for the first patient treatment with 4 He ion beams at HIT, we compared the predictions of three existing RBE models to biological measurements in vitro with monoenergetic beams and in clinically-relevant scenarios, as well as highlighting the inter-and intra-model variations as a function of tissue type, dose level, LET d , depth and beam configuration in silico. For the in vitro study, a cell line exhibiting substantial radio-resistance was selected for irradiation with both pristine beams and clinical-like fields. These more radio-resistant tissues (α/ β < 4Gy) are of particular interest considering they make up only~5% of the available experimental data in the literature for 4 He ion beams. In addition to in vitro study, patient treatment plans were calculated and compared, applying the various 4 He RBE schemes in place of a constant RBE [20]. The three published models for RBE prediction with 4 He ion beams investigated in this study are as follows: a data-driven phenomenological model (DDM) [13,14] and two biophysical models featuring the Local Effect Model (LEM, version IV) [21] and the modified Microdosimetric Kinetic Model (MKM) [22,23]. With a long-term outlook in mind for 4 He RBE study and clinical integration, this work can serve as a foundation for clinical decision-making regarding effective dose calculation, in preparation for the first 4 He ion beam therapy patient treatments in Europe.

Experimental investigations
Cell culture and clonogenic assay Murine renal adenocarcinoma cells (Renca ATCC® CRL-2947™) were cultured in RPMI-1640 Medium (Gibco, Germany) supplemented with 10% heat-inactivated Fetal Bovine Serum (FBS, Millipore, Germany) and 1% Penicillin/Streptomycin (Gibco, Germany) at 37°C and 5% CO 2 atmosphere. Clonogenic cell survival assay, i.e. seeding, irradiation, incubation and read-out, was performed as previously described using 96-well plates [24]. Image acquisition took place with the IncuCyte® System (Essen BioScience, UK) for colony counting. A baseline characterization of the cell line was performed separately prior to experiment A (pristine peaks) and experiment B (SOBPs), which involved photon irradiation delivery (LINAC, 6 MV, Artist Siemens) with dose levels of 1, 2, 4, and 8 Gy for determination of the LQ parameters (α x and β x ).

Irradiation with Monoenergetic beams
To most closely resemble track segment conditions, cell were irradiated with monoenergetic 4 He beams (E 4-He = 56.66 MeV/u, d BP = 25.9 mm) in experiment A. Two sets of biological measurement points were taken at 6 mm and 12 mm water-equivalent depth (WED). Cell-kill measurements were collected for the pristine beams at dose levels of approximately 0.25, 0.5, 0.75, 1.0, 1.5, 2.0, 2.5, and 3.0 Gy. Dosimetric measurements were performed using a Farmer ionization chamber (TM30010, PTW, Freiburg) for validation of FLUKA MC predictions.

Irradiation with SOBPs
For investigating clinical-like conditions, the same plate configuration was used as in the base-line photon irradiations. Experiment A and B involved 96-well plates positioned against various thicknesses of PMMA, such that each plate corresponded to a specific depth (and hence, LET d ) in the SOBP irradiation [24], with positions of 3.0 cm (p1), 5.98 cm (p2), 7.61 cm (p3), and 8.35 cm (p4) in PMMA. WED values were calculated using a multiplicative factor of 1.165 and are highlighted in Fig. 1 (right panel). SOBP plans were physically optimized in water for the following doses in the 12 cm × 8 cm × 4 cm target region centered at 8 cm depth: 0.5, 1.0, 2.0, 3.0, 4.0 and 6.0 Gy. The 96-well plate geometry with corresponding material composition was integrated into the FLUKA MC simulation.

Models and MC simulation
Modeling the relative biological effectiveness of 4

He ion beams
Biological dose prediction begins with modeling cell survival (S), traditionally described as a linear-quadratic (LQ) trend, with α and β representing the linear and quadratic coefficients, respectively, as a function of physical dose (D). The ratio of the linear and quadratic coefficients, (α/ β) x , is often referred to as a description for the sensitivity of the cell line when exposed to photon radiation (x). The RBE is a multifunctional quantity defined as the isoeffective dose ratio between a reference radiation (D x ) and a particle radiation (D p ), traditionally modeled as a function of three parameters: (α/β) x , LET and D x . Biological (or effective) dose (D RBE ) is defined as the product of the RBE and the physical dose.
Within the LQ framework, we can determine a dependency of RBE on (α/β) x , the helium absorbed dose, RBE α and R β [13,14]: In the next sections, the expressions for RBE α and R β per the three models will be introduced. In the case of the LEM, the LQ approximation for the photon response is valid up to threshold dose D t , which marks the transition dose at which the survival curve for photon irradiation is assumed to have an exponential shape with the maximum slope S max = α x + 2β x D t [25]. In this work, the dose levels have been chosen within the range of LQ applicability, i.e. < D t .
The predictions of the three RBE models have been assessed by comparing RBE α and R β as a function of LET, and the RBE values as a function of LET and dose for two tissue types irradiated with 4 He ion beams. Parameters characterizing the hypothetical tissues considered for this study are reported in Table 1 and labeled water case. The (α/β) x values were selected similar to recent works [26] to represent late-responding tissues (low (α/β) x from 2 to 3 Gy), and early-responding normal tissues and most common tumors (high (α/β) x from around 10 Gy).

Data-driven LET-based model
A phenomenological model for RBE with 4 He ion beams was developed by fitting in vitro experimental data available in the literature in Mairani et al. 2016a [11] and refined in Mairani et al. 2016b [12]. For RBE α , the following parameterization has been introduced: where L* represents the rescaled 4 He LET [13]: LET x and LET 60Co are, respectively, the LET of photon under study and of the reference 60 Co. The parameters used in eq. 3 are as follows [12]: k 0 = 8.924 × 10 − 2 Gy − 1 and k 1 = 3.368 × 10 − 1 μm·keV −1 , and k 2 = 2.858 × 10 − 5 μm 2 ·keV − 2 . For R β , we have introduced an LET-dependent parameterization fitting the running averages of R β as function of LET: The coefficients for the R β parameterization are b 0 = 2.66, b 1 = 62.61 keV μm −1 and b 2 = 48.12 keV μm −1 .
For comparison in track-segment conditions, we have assumed L * = LET while for the clinically-relevant scenarios and in vitro studies, we used 6 MV photon beams as a reference radiation for calculating rescaled L * values.

Modified Microdosimetric kinetic model (MKM)
In the modified MKM [22,23], for any radiation quality, RBE α is expressed as a function of the saturation-corrected dose-mean specific energy of the domain delivered in a single event z Ã 1D divided by the (α/β) x ratio: z Ã 1D depends on z, the specific energy, and z sat , the saturation-corrected specific energy which accounts for the decrease of RBE due to the overkilling effect for high specific energy values [27]. z depends on the radius of the domain (R d ) while z sat depends R d and the radius of the cell nucleus (R n ) [22]. MKM input parameters (R d and R n ) have been tuned in a previous work [22] to reproduce an in vitro experimental biological database of Table 1 Photon parameters applied during the in silico investigations. The D t parameter is required for LEM calculations only  Table 2 Clonogenic cell survival LQ fit parameters for photon (α x and β x ) and helium ion beam (α and β) irradiation using the Renca cells in vitro with corresponding LET d derived from MC simulation. Data for both experiment A (pristine peaks) and experiment B (SOBPs) are provided Local effect model (LEM) The LEM-version IV developed by the GSI Helmholtz Centre for Heavy Ion Research (Darmstadt, Germany) [21] relates the biological response directly to the double-strand breaks pattern and has been benchmarked by its developers in various publications [10,21]. The LEM intrinsic α z tables are obtained using the PT RBE Generator software by Siemens which is available at HIT, while for β z , we have used the approximation β z = ( s max − α z )/(2D t ), with negative values found at high LET forced to zero [25]. The LQ parameters are calculated at different energies applying the low dose approximation, which describes how to link the input LEM-calculated intrinsic microscopic parameters, α z and β z , to the macroscopic values, α and β. The initial RBE can be written as: with R β as: d 1 is the dose deposited by a single particle traversal [29,30].

MC simulation of the in vitro study
For both experiment A and B, the target (96-well plate irradiation system) was incorporated into FLUKA MC, including a detailed geometry of the HIT beam-line [31], for validating the biological dose models against experimental measurements. Once biological measurements were acquired, simulations were executed to score physical dose and LET d , as well as the various biological parameters necessary for D RBE using the DDM, MKM and LEM. With a detailed geometry of the 96-well plate target, parameters were scored on a per well basis to reduce physical and biological uncertainties during evaluation of measurement and simulation outcomes, as shown in Fig. 1. Cell survival and, in turn, RBE results were compared to MC prediction to validate enhanced cell-kill with increased LET d for helium ions and to evaluate model performance.

Patient studies and validations
Retrospective study: patient treatment planning and forward computation of D RBE In this work, the MC-based treatment planning tool (MCTP) [32,33] is employed to create biologically optimized treatment plans and to perform forward dose calculation for retrospective study. The MCTP relies on FLUKA's capability to describe the interaction and transport of radiation with matter for 4 He ion beams and is coupled with both biophysical and phenomenological RBE models for 4 He. FLUKA has been benchmarked against dosimetric data, demonstrating overall a satisfactory agreement [11].
The MCTP uses dosimetrically commissioned scanned pencil beams as available at HIT [34]. The data-driven RBE model has been used for treatment plan optimization. The MCTP tool relies on externally generated databases for each biological effect model to calculate RBE and D RBE values [37,38]. To properly calculate effective dose for helium ion beams, Z = 2 primary particles and secondary fragments as well as Z = 1 secondary fragments must be scored separately. Hence, both the DDM and a phenomenological model for Z = 1 were used during biological dose weighting of Z = 2 and Z = 1, respectively [35].
MCTP-based plans have been calculated to achieve a homogeneous three-dimensional D RBE of 2.0 Gy (RBE) and 4.0 Gy (RBE) in the target region with a single field and a two opposing fields arrangement in water. Two targets were chosen: rectangular parallelepiped volumes of 6 cm × 6 cm × 6 cm and 3 cm × 3 cm × 3 cm centered at 12.5 cm water-equivalent depth. FLUKA MC scoring for physical and biological quantities was performed in voxels of 2 mm × 2 mm × 2 mm. The lateral PB spacing was 3 mm while the depth separation between Bragg peak positions of two consecutive energy slices was 2 mm. The plans have been calculated assuming two representative tissues with (α/β) x of 2 Gy and 10 Gy as reported in the first two rows of Table 1. MCTP-based plans for two patients (previously treated with protons at HIT) were simulated using one and two 4 He ion beam portals. Beam configurations for the head and prostate case involved a single field (superior-inferior direction) and parallel opposed fields (anterior-posterior / posterior-anterior direction), respectively. To achieve dose homogeneity in the target of the head case, a ripple filter has been used for broadening the beam longitudinally [36]. The lateral PB spacing was 3 mm while depth separation between two consecutive energy slices was 3 mm. FLUKA MC scoring was performed in voxels of 1 mm × 1 mm × 3 mm. The planned doses were 54 Gy (RBE) in 27 fractions and 66 Gy (RBE) in 20 fractions for the head and the prostate cases, respectively [35], applying the clincal fractionation scheme used at HIT with proton beams.
Forward re-computation of the optimized plans have been carried out to investigate the variation of the D RBE as a function of the depth, applying the biophysical models previously described. LET d distributions were additionally scored for dose values larger than 5% of the maximum Z = 2 dose. For dose distribution characterization in the target of the SOBPs, equivalent uniform dose (EUD) was applied [37]. We have calculated EUD as follows [38]: where S is the mean survival in the target. For the patient cases, we have also analyzed the D RBE volume histograms (D RBE VH).

Development and validation of an analytical biological dose calculation engine: FRoG
Once patient case dose calculation was established via biological dose models coupled with FLUKA, validation of the fast (GPU-based) analytical dose engine, FRoG, was performed [44,45]. Physical and biological parameter database generation took place using FLUKA MC simulation. Corresponding biological parameters for DDM (α He and β He ), LEM (α He and β He ), and MKM (z Ã 1D ) were scored as a function of depth, along with the necessary physical parameters (dose and LET d ). The physical and biological tables were incorporated into the FRoG platform, enabling multi-tissue (variable (α/β) x ) dose calculation for the three biological dose models. The glioma patient plan was executed in FRoG for comparison with the gold standard FLUKA MC.
All patients records were anonymized prior to the study, obtained with informed consent and handled following the Helsinki Declaration. All methods were approved by the Heidelberg University Medical Faculty, following applicable guidelines and regulations of the institution.

Investigating model dependencies in silico: SOBPs and patient cases
Clinically relevant scenarios were used to further characterize model variations. Figure 2 presents RBEweighted dose (D RBE ) for the SOBPs, calculated via MC simulation, as a function of depth in water for the three investigated models, as well as physical dose and LET d . RBE variation and %Δ RBE are also visualized in the following middle and lower panels, respectively. The SOBP plan, biologically optimized using the DDM, was applied to reach a biological dose level of 2 Gy (RBE) and for two tissue types, exhibiting (α/β) x of 2 Gy and 10 Gy, experimental surrogates for testing radio-resistant and radio-sensitive tissues, respectively. A similar investigation was executed for an irradiation plan with two opposing fields, as shown in Additional file 1: Figure S1. For quantification of global difference in the target between the various models, EUD calculations for the SOBPs studied in silico are provided in Additional file 1: Table S2 and S3.
Model dependencies in clinically-relevant scenarios: patient cases In Fig. 3, an investigation of RBE model performance of a prostate cancer patient in silico is displayed. The MCTP calculated D RBE distribution for the pelvic case applying the DDM and LET d distribution are shown as well as dose difference (Δ Gy(RBE) ) from the reference when performing forward calculations with LEM and MKM. The physical dose volume histogram (DVH) and biological dose volume histogram (D RBE VH) for the PTV and rectum, chosen as a representative organ at risk (OAR), are displayed in the bottom panel. DVH statistics for the PTV in terms of D 50% , D RBE-50% and the inhomogeneity coefficient I 5% = (D 5% -D 95% ) / D RBE,p have been analyzed. D RBE-50% , D RBE-5% , and D RBE-95% represent the biological dose received by 50%, 5% and 95% of the PTV volume in the cumulative D RBE VH, respectively. D RBE,p is the prescribed biological dose. I 5% evaluates the biological dose gradient introduced in the PTV by performing forward calculation of the patient plans with the various RBE models. LEM resulted in −5.7% lower D 50% , while applying the MKM yielded 8.3% higher D 50% . The I 5% values were, respectively,~12% for MKM, and1 0% for both LEM and the reference (D RBE calculated with DDM). The D 5% for the rectum was 50.

Validating RBE models in a clinical platform: FRoG
A glioma patient case is displayed in Fig. 4 for RBE evaluation and validation of a fast analytical dose calculation engine (FRoG). FRoG calculation run-time for the glioma patient (yielding D and D RBE applying DDM, MKM and LEM) was 142 s, a time gain factor of~225 when compared to MC simulation using a 300 node CPU-cluster. The MCTP calculated D RBE distribution for the head case applying the DDM and the resulting LET d distribution are shown as well as dose difference Δ Gy (RBE) from the reference when performing forward calculations with (c) LEM and (d) MKM. For the LEMand MKM-based forward biological dose calculations, D 50% for the PTV is 1.5% higher and −3.7% lower, respectively, than the reference. Larger I 5% values were found for LEM and MKM of~18% and~14%, respectively, relative to the reference of~13%. The greatest variations between the models occur for the normal tissue with (α/β) x = 3.1 Gy, outside of the PTV, especially in the distal region where the highest LET components of the distribution are prevalent. For the glioma patient case, there are no OARs in proximity of the target.
As shown in Fig. 4, DVH and D RBE VH plots between FRoG and FLUKA are in good agreement. The percent absolute deviations in D 50% and D RBE-50% for the PTV between FLUKA and FRoG for physical dose (D phys ) and the three biological doses are as follows: 0.2, 0.4, 0.4, 0.6%, for D phys , D DDM , D LEM and D MKM , respectively. Further details regarding DVH and D RBE VH statistics are provided in Additional file 1: Table S1.

Experimental evaluation of the RBE models
Enhanced cell-killing was observed in the biological measurements of experiment A for higher LET d (~15 keV·μm − 1 ) compared to lower LET d (~6 keV·μm − 1 ). Figure 5 displays both the experimental findings (points with error bars) and FLUKA MC-coupled RBE model predictions for cell survival and RBE, as well as percent difference in RBE (%Δ RBE ) of the three models against experimental data. Linear quadratic (LQ) fitting of the cell survival data from photon irradiations with the 6MV LINAC yielded α x = 0.034 Gy − 1 and β x = 0.018 Gy − 2 , for an (α/β) x of 1.79 Gy. For the lower LET d condition, LEM exhibited the most stable prediction of RBE as a function of dose below 1.5 Gy with %Δ RBE < 5% but consistently underestimates RBE. On the other hand, DDM and MKM yielded better RBE predictions from 1.5 Gy and above. For the higher LET d condition, DDM and MKM predicted with the highest relative accuracy within the studied dose range, with %Δ RBE < 5% up to 2 Gy. LQ-fit parameters for two LET d conditions are listed in the Table 2. Regarding outcome of experiment B, initial investigation of cell-kill response to photon irradiation yielded α x = 0.050 Gy − 1 and β x = 0.023 Gy − 2 , for an (α/β) x of 2.17 Gy, which is on average 0.38 Gy higher than the (α/ β) x found in experiment A. Figure 6.a displays the cell survival versus dose for the four LET d conditions (~5 keV·μm − 1 ,~10 keV·μm − 1 ,~15 keV·μm − 1 ,~27 keV·μm − 1 ) within a clinically relevant dose range (D phys ≲3 Gy). For both model predictions and experimental data, a dose dependence in RBE was observed in all cases. In general, DDM and MKM performed best for both higher and lower LET d conditions in the studied dose range, consistent with findings from the monoenergetic beam experiment. RBE predictions for all three models agreed within ±5% of the experimental data for the two highest LET d conditions (~15 keV·μm − 1 and~27 keV·μm − 1 ), especially DDM and MKM for dose levels > 2 Gy. For 2 Gy, %Δ RBE for the four LET d conditions (in ascending order) were roughly, + 3.7%, − 1.9%, − 1.9%, − 4.4% for DDM, − 1.7%, − 5.3%, − 3.4% and + 0.9% for LEM, and − 4.1%, − 1.1%, − 1.1% and − 4.8% for MKM. For the lower LET condition of~5 keV·μm − 1 (entrance channel measurement), all models produced RBE predictions within ±5 −10%, reaching~1.3 for 0.5Gy,~1.25 for 1 Gy,~1.18 at 2 Gy and stabilizing to~1.1 for the higher doses. As for the LET d conditions found in the target (~10 keV·μm − 1 ,~15 keV·μm − 1 ,~27 keV·μm − 1 ), representing a low, mid and high range LET d for therapeutic helium ion beams, respectively, greater variability was observed as a function of dose, especially for doses < 2 Gy. For 1 Gy, observed RBE values were~1.8,~2.2,~2.8 for the low, mid and high LET d conditions in the target. At 4 Gy, RBE values decreased to1 .3,~1.5,~1.8 for the low, mid and high LET d conditions.

RBE model assessment
To best interpret the biological models for 4 He ion beams, one must begin with a survey of their dependencies in track-segment conditions, i.e. monoenergetic beam case disregarding contributions from a mixed radiation field. In  Figure 7.a shows the comparison of RBE α (top) and R β (bottom), for mono-energetic 4 He ion beams as a function of LET for two tissues, (α/β) x = 2 Gy (left panels) and 10 Gy (right panels), representing two distinct tissue types with differing responses to radiation. Comparison of these cases shows RBE α and (α/β) x are negatively correlated. As particle LET increases, an a b Fig. 6 Clonogenic assay for clinical-like fields (SOBPs) for the Renca cell line in experiment B. MC simulation estimated LET d values of biological measurement were~5 keV·μm −1 ,~10 keV·μm −1 ,~15 keV·μm −1 ,~27 keV·μm −1 . FLUKA-coupled biophysical and phenomenological models predicted cell survival (a) and corresponding RBE (b) with varying degree of accuracy as a function of dose. The dotted and solid black line represent the LQ-fit of the Renca cells photon irradiation and 4 He irradiation, respectively. LQ-fit parameters for the four LET d conditions are listed in the Table 2 upward trend for RBE α as a function of LET is observed, until a saturation point, where the RBE α plateaus prior to fall-off. In general, this fall-off is more prominent and occurs at a lower LET range in lower (α/β) x tissues. For lower LET, the largest inter-model variation occurs for the (α/β) x = 2 Gy case between LEM and the other two models, while for the higher LET region, all models exhibit a varying response. For (α/β) x = 10 Gy, the models yield similar predictions for LET values lower than about 20 keV·μm − 1 . The location of RBE α maximum changes as a function of the model applied.
Regarding R β , the models assume or predict different behaviors as function of LET. In the MKM [28], R β is assumed to be unity, i.e. β He = β x , while for the single-hit based version of LEM applied in this work [21], R β decreases as LET increases. In the LET-based DDM approach, R β increases with LET until reaching a maximum at~63 keV·μm − 1 and then drops to zero for LET > 100 keV·μm − 1 . For the data-driven approach, R β is independent of (α/β) x , and therefore it's behavior is consistent between tissue types. These differences in expressing R β lead to significant variations among the models which, in part, reflect the large experimental uncertainties of the available experimental in vitro data [13].
RBE versus LET for the two tissues at physical dose levels of 2 Gy (left column) and 4 Gy (right column) are depicted in Fig. 7 b. As expected, the RBE initially increases with LET, reaches a maximum and then decreases. The RBE decreases for increasing dose mainly for low (α/β) x , and increases for decreasing (α/β) x of the tissue. RBE results at lower LET and higher LET are presented as a function of physical dose for the two tissues. The chosen LET values are representative of the LET d values found in the entrance channel and in the middle of an SOBP, respectively, for the two opposing beam fields arrangement depicted in Additional file 1: Figure  S1. For clinical targets like an SOBP, one must consider a mixed radiation field with a complex LET spectrum, rather than a single LET value as in the case of an ion in the track-segment condition.
As expected, an enhanced RBE is observed at lower doses for all models, and this trend is more pronounced for lower (α/β) x tissues. For the low LET condition, LEM predicted a limited RBE variation within the studied dose level, between maximum and minimum values, of about 20% and of about 4% for (α/β) x = 2 Gy and (α/β) x = 10 Gy, respectively. For 15 keV·μm −1 and for (α/β) x = 2 Gy, MKM and the DDM approach resulted in roughly the same predictions, while for (α/β) x = 10 Gy the DDM estimated about 15% higher RBE. In order to reduce model-related uncertainties in the target region, assuming 15 keV·μm −1 is a representative LET d value for Z = 2 in the target, one could use hypo-fractionated treatments (D RBE > 4 Gy (RBE)) where variations in RBE prediction decrease. In addition, hypo-fractionated treatments reduce the impact of precise (α/β) x value assignment for target tissues on RBE determination. On the other hand, hypo-fractionation may diminish the therapeutic window by reducing the ratio of the target RBE compared to the entrance channel (i.e. tumor to normal tissue effective dose ratio). With typical peak-to-plateau dose ratio of~2 for 4 He ion beams and assuming a dose value of 4 Gy in target, RBE predictions (averaging over the three models in this work) are as follows:~1.1 for 4.0 keV·μm −1 and~1.45 for 15 keV·μm −1 in low (α/β) x tissues, and~1.1 for 4.0 keV·μm −1 and~1.35 for 15 keV·μm −1 in high (α/β)x tissues. Conversely, standard fractionation schemes (~2 Gy (RBE) target doses) can enhance the peak-to-plateau ratio.
Close examination of the R β component for the DDM reveals that for LET of~4 keV·μm −1 , R β converges to~0 .6, while for 15 keV·μm −1 R β approaches~1. As described in previous works [13,14], R β parameterization was obtained by a convenient parameterization which fits the running averages of the experimental data, neglecting any (α/β) x dependencies due to the large uncertainties effecting the β term. Recent works develop a phenomenological model for proton beams from in vitro data following a similar approach to R β handling by assuming a negligible (α/β) x dependency [35,46]. With DDM, parameter fittings are merged to a relatively small amount of data using a running average and thus, this work can shed light on RBE model performance in regions where data is sparse and predictions exhibit large uncertainties. Moreover, existing experimental data is especially scarce for low (α/β) x values (< 3 Gy) [14], where the largest RBE values are expected and the highest variations among the models occur. Further data for low (α/β) x tissues and for clinically-relevant dose levels, especially in standard fractionation regimes (D RBE <~3 Gy (RBE)), are essential for benchmarking the predictive power of these RBE models.

Experimental benchmarking (in vitro)
RBE model benchmarking through in vitro experimentation with a low (α/β) x cell line was the next logical step to verify the significant RBE enhancement observed in the models for dose levels < 4 Gy, a clinically relevant range bearing in mind the typical fractionation size for proton beams of~2 Gy (RBE). Qualitatively, the study investigated both lower LET d (< 10 keV·μm −1 ) and higher LET d (≥10 keV·μm −1 ) values, pertinent endpoints for both normal tissue complication and tumor control probability (TCP). In addition, critical structures surrounding or distal to the target are also associated with the highest LET d values in the study. It is important to note, however, that the in vitro data available in the literature is solely based on cell-kill of tumor tissues with RBE as the end point. Therefore, the models provide insight into RBE from the perspective of TCP rather than normal tissue response, which requires the immortalization of normal cell lines to investigate relevant end points [47].
For RBE prediction versus measurement in experiment A (Fig. 5), LEM exhibited the highest accuracy for low LET d at dose levels <2Gy, while MKM and DDM performed best for the higher doses. For higher LET d conditions, MKM and DDM both outperformed LEM in predictive power, with local %Δ RBE between~1% and8 %, as the dose increases. Although direct comparison of the track-segment condition in silico study shown in Fig.  7 and the monoenergetic beam in vitro study is incompatible due to the oversimplification of LET d (neglecting mixed field spectra) and the inherently non-linear relationship of RBE and LET, general trends between the models are consistent.
As for investigations in experiment B (clinical-like fields in Fig. 6), interpretation becomes more convoluted when considering the complex mixed radiation field. In general, DDM and MKM demonstrated the lowest local |%Δ RBE | of < 10%, overall. As anticipated, |%Δ RBE | decreased with increasing dose for all three models. Disagreement in the lower LET d condition can be explained by the scarce amount of data for low LET d , especially with cell lines with (α/β) x < 3 Gy, which suggests that further in vitro study and tweaking of the models could yield improved RBE predictions. Nevertheless, 5% to 10% predictive power for RBE in the target region is acceptable considering the uncertainty of the reference photon sensitivity measurement. For the entrance channel condition in Fig. 6, all three models (especially DDM) tend to overestimate RBE for < 1 Gy, a typical fractionation treatment dose range, offering a conservative estimate for normal tissue in the plateau region.
DDM depends only on the (α/β) x ratio while the MKM, instead, depends also on the absolute value of β x , which contributes in the determination of z sat [22]. Low β x values result in a reduced saturation coefficient, leading to RBE enhancement. To further shed light on this point, calculations were performed with the two fields arrangement applying (α/β) x = 2.0 Gy, planned D RBE = 4 Gy (data not shown) and β x = 0.02 Gy − 2 , finding consistently higher D RBE values (about 8%). In contrast, LEM depends on multiple parameters, including α x , β x and D t . By varying α x and β x by 25% but maintaining the same (α/β) x , no measurable dependence of RBE α was found for clinicallyrelevant LET values using carbon ion beams, with a limited effect on the RBE at 10% survival [48].

Clinical outlook
Regarding patient dose calculation, LET d prediction for the prostate case was in line with the findings from the SOBP study; however, the head case plan exhibited lower LET d values since the energy spread of the beam is increased by the ripple filter (RiFi) to reduce BP sharpness for clinically acceptable target dose homogeneity. Furthermore, FRoG calculated physical and biological dose distributions were in good agreement with FLUKA MC and well within clinically acceptable tolerances. At HIT, both the MCTP and FRoG dose engine are functional for helium ion beam therapy, enabling future treatment planning comparison and robust RBE optimization studies necessary before and during clinical trials, as performed in previous works for carbon ions [49]. In addition, the FRoG platform will support the development and validation of the first analytical TPS for helium ion beams, providing multiple biological models for clinical research.
As HIT prepares for clinical translation of 4 He, the findings and efforts in this work may serve as a starting point for clinical decision making. Currently, there is no official consensus as to which RBE model for helium ions is best suited for treatment and whether a single tissue approximation for biological dose prediction will be used as done with carbon ions. In light of these issues, the FRoG platform includes all three models presented in this work, as well as tissue-dependent biological dose calculation, providing valuable insight into radiological uncertainty during treatment planning. Regarding optimization of a next generation TPS for particle therapy, advanced optimization strategies are recommended considering the large uncertainties associated with biological modeling and the lack of evidence supporting in vitro model applicability to in vivo settings [50]. With technqiues like RBE/LET gradient minimization in the target, constant over-or under-estimation of D RBE could be detected in an initial dose-escalation phase. At HIT, a systematic clinical investigation with an initial group of patients is anticipated to observe and analyze clinical outcome.
All presented RBE models are based on the same set (or sub-set) of the published biological in vitro data, used repeatedly for model tuning and benchmarking purposes. In vivo data is sparse at best and rarely used to verify the models' predictions [51]. The experimental and intrinsic uncertainties in the data constrain the confidence in these models to a degree which is less than clinically desirable, yielding model fits with significant variation. It is worth noting here that the agreement of the LEM used for this study with respect to the other models might further improve if the same set of in vitro data would have been used for tuning the LEM, as done for the DDM and MKM. These findings suggest that systematics in RBE predictions in the high dose region for clinical 4 He ion treatment fields due to different choices of RBE modelling approaches can be restricted to be mostly within 10% to 15% when tuning the parameters of the RBE models to the same (or a similar) set of the available in vitro cell data for 4 He ions.
In turn, this may imply that systematic uncertainties in the prediction of RBE for helium ions for clinical scenarios are not primarily dominated by the choice of the RBE model, but instead dictated by the choice of the in vitro dataset and methodology used for tuning the RBE model parameters. Similar conclusions might hold true for RBE models of higher Z ion species. Additional systematic RBE uncertainties arise from differences between in vivo and in vitro data; however, due to their scarcity, in vivo and clinical data are hardly used to tune RBE models, but rather for validation of commonly established RBE models [52], exception being the neutron-equivalent scaling point used for carbon ions [53,54]. Previous works also propose application of clinical data for RBE model tuning in addition to in vitro and in vivo measurements [55].
For helium ions, it is certainly challenging to make definitive statements about RBE considering the lacking of experimental data. To reduce RBE model uncertainties for 4 He, collecting additional evidence, especially in vivo, is recommended before clinical application. However, the differences in RBE predictions found in this study for the three presented models are similar to the RBE variation for in vitro data in proton beams, which are typically knowingly accepted when assuming RBE = 1.1 [47]. Ultimately, the choice of model and tissue type for biological dose optimization is a clinical decision to ensure the most safe and effective patient treatment and care possible.

Conclusion
Before the start-up of a 4 He ion beam therapy program, a comprehensive evaluation of the variable RBE and the associated models is critical. The main dependencies of three RBE models for 4 He ion beam therapy were studied in silico and validated against in vitro experimentation with a radio-resistant tumor cell line. Clinically relevant uncertainties were observed, especially for low (α/β) x values where the available literature data are scarce. The observed uncertainties between the models as well as variability of RBE as a function of its dependency (especially for low (α/β) x tissues commonly treated with particle therapy) suggest that the selection, refinement and validation of either a biophysical/mechanisticor phenomenological-based approach are essential prior to clinical translation of helium ion beam therapy.

Additional file
Additional file 1: Supplementary data analysis for biological dose prediction using 4 He ions, including SOBPs for a parallel opposed beam plan (two-field), DVH statistics for FRoG against FLUKA MC for the two patient cases, and EUD calculations comparing the three investigated RBE models. (DOCX 1421 kb)