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

Abstract

Background

Helium (⁴He) 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 ⁴He 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 ⁴He ion therapy program at the Heidelberg Ion-beam Therapy Center (HIT), we aimed to characterize the biophysical phenomena of ⁴He ion beams and various aspects of the associated models for clinical integration.

Methods

Characterization of biological effect for ⁴He 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 ⁴He 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 ⁴He ion beam therapy to the clinic.

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 (¹H) and carbon (¹²C) ion beams, the leading modalities in hadrontherapy, are applied clinically, with two additional particle species available for experimentation: oxygen (¹⁶O) and helium (⁴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, ⁴He ion beams remain clinically unexploited.

It is well known that, experimentally, heavier ions exhibit greater biological damage and consequently, the biophysical properties of ⁴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 ⁴He ion beams, which marks the world’s first clinical application of raster-scanning ⁴He ion therapy. Over the past decade, substantial efforts have been made at HIT to characterize ⁴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 ⁴He ion beams; however, research-based tools were recently introduced or updated to allow planning with ⁴He ion beams [10, 14, 16].

Relative to the clinical standard photons and protons, ⁴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 ⁴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 ⁴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 ⁴He ion beams. In addition to in vitro study, patient treatment plans were calculated and compared, applying the various ⁴He RBE schemes in place of a constant RBE [20]. The three published models for RBE prediction with ⁴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 ⁴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 ⁴He ion beam therapy patient treatments in Europe.

Methods and materials

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₂ 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 ⁴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.

Fig. 1

Left: cross-section of schematic for the 96-well plate geometry and composition in FLUKA MC for experiment A and B. Right: central line profile through physically optimized SOBP plan for experiment B, displaying both physical dose and LET_d. The biological measurement positions are designated by the highlighted regions (p1, p2, p3, p4)

Table 1 Photon parameters applied during the in silico investigations. The D_t parameter is required for LEM calculations only

Models and MC simulation

Modeling the relative biological effectiveness of ⁴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]:

$$ \mathrm{RBE}\left({\left(\frac{\upalpha}{\upbeta}\right)}_{\mathrm{x}},\mathrm{D},{\mathrm{R}\mathrm{BE}}_{\upalpha},{\mathrm{R}}_{\upbeta}\right)=-\frac{1}{2\mathrm{D}}{\left(\frac{\upalpha}{\upbeta}\right)}_{\mathrm{x}}+\frac{1}{\mathrm{D}}\sqrt{\frac{1}{4}{\left(\frac{\upalpha}{\upbeta}\right)}_{\mathrm{x}}^2+{\mathrm{R}\mathrm{BE}}_{\upalpha}{\left(\frac{\upalpha}{\upbeta}\right)}_{\mathrm{x}}\mathrm{D}+{\mathrm{R}}_{\upbeta}{\mathrm{D}}^2} $$

(2)

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β_xD_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 ⁴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).

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

Data-driven LET-based model

A phenomenological model for RBE with ⁴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:

$$ {\mathrm{RBE}}_{\upalpha}=1+\left[{\mathrm{k}}_0+{\left(\frac{\upalpha}{\upbeta}\right)}_{\mathrm{x}}^{-1}\right]{\mathrm{k}}_1{\mathrm{L}}^{\ast }{\mathrm{e}}^{-{\mathrm{k}}_2{\mathrm{L}}^{\ast 2}} $$

(3)

where L* represents the rescaled ⁴He LET [13]:

$$ {\mathrm{L}}^{\ast }=\mathrm{LET}-{\mathrm{L}\mathrm{ET}}_{\mathrm{x}}+{\mathrm{L}\mathrm{ET}}_{{}{}^{60}\mathrm{Co}} $$

(4)

LET_x and ${\mathrm{LET}}_{{}{}^{60}\mathrm{Co}}$ are, respectively, the LET of photon under study and of the reference ⁶⁰Co. The parameters used in eq. 3 are as follows [12]: k₀ = 8.924 × 10^− 2 Gy^− 1 and k₁ = 3.368 × 10^− 1 μm·keV⁻¹, and k₂ = 2.858 × 10^− 5 μm²·keV^− 2. For R_β, we have introduced an LET-dependent parameterization fitting the running averages of R_β as function of LET:

$$ {\mathrm{R}}_{\upbeta}={\mathrm{b}}_0{\mathrm{e}}^{\left[-{\left(\frac{{\mathrm{L}}^{\ast }-{\mathrm{b}}_1}{{\mathrm{b}}_2}\right)}^2\right]} $$

(5)

The coefficients for the R_β parameterization are b₀ = 2.66, b₁ = 62.61 keV μm⁻¹ and b₂ = 48.12 keV μm⁻¹ .

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 ${\mathrm{z}}_{1\mathrm{D}}^{\ast }$ divided by the (α/β)_x ratio:

$$ {\mathrm{RBE}}_{\upalpha}=1+{\left(\frac{\upalpha}{\upbeta}\right)}_{\mathrm{x}}^{-1}\bullet {\mathrm{z}}_{1\mathrm{D}}^{\ast } $$

(6)

${\mathrm{z}}_{1\mathrm{D}}^{\ast }$ 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 initial RBE. The resulting best fit values of R_d = 0.3 μm and R_n = 3.6 μm obtained in Mairani et al. 2017 have been used in this work without further adjustments. For the R_β term, it is assumed [28]:

$$ {\mathrm{R}}_{\upbeta}=1. $$

(7)

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:

$$ {\mathrm{RBE}}_{\upalpha}=\frac{1-{\mathrm{e}}^{-{\upalpha}_z{\mathrm{d}}_1}}{\upalpha_x{\mathrm{d}}_1} $$

(8)

with R_β as:

$$ {\mathrm{R}}_{\upbeta}={\left(\frac{\upalpha}{\upalpha_z}\right)}^2\left(\frac{\upbeta_z}{\upbeta_x}\right) $$

(9)

d₁ 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 ⁴He ion beams and is coupled with both biophysical and phenomenological RBE models for ⁴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 ⁴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]:

$$ \mathrm{EUD}=-\frac{1}{2}{\left(\frac{\upalpha}{\upbeta}\right)}_{\mathrm{x}}+\sqrt{\frac{1}{4}{\left(\frac{\upalpha}{\upbeta}\right)}_{\mathrm{x}}^2-\frac{\ln \left(\overline{\mathrm{S}}\right)}{\upbeta_{\mathrm{x}}}} $$

(10)

where $\overline{\mathrm{S}}$ is the mean survival in the target. For the patient cases, we have also analyzed the D_RBE volume histograms (D_RBEVH).

Following the previous works [14], tissue parameters for the prostate case were set to α_x = 0.036 Gy⁻¹ and β_x = 0.024 Gy⁻² for (α/β)_x = 1.5 Gy [39, 40]. For the surrounding healthy tissues, (α/β)_x = 3.1 Gy with α_x = 0.089 Gy⁻¹ and β_x = 0.0287 Gy⁻² was applied [41]. For the head patient case, for the planning target volume (PTV) assuming a glioma tumor, we have used α_x = 0.077 Gy⁻¹ and β_x = 0.009 Gy⁻², yielding (α/β)_x = 8.6 Gy [42] while for the rest of the brain, we have assumed α_x = 0.0499 Gy⁻¹ and β_x = 0.0238 Gy⁻², yielding (α/β)_x = 2.1 Gy [43]. Further details regarding these values are provided in Table 1.

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 (${\mathrm{z}}_{1\mathrm{D}}^{\ast }$) 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.

Results

Investigating model dependencies in silico: SOBPs and patient cases

Clinically relevant scenarios were used to further characterize model variations. Figure 2 presents RBE-weighted 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.

Fig. 2

Biologically optimized SOBPs at 2 Gy (RBE) using the data-driven model (DDM) assuming two distinctive tissue types with (α/β)_x = 2 Gy (left) and (α/β)_x = 10 Gy (right) are displayed as function of the depth in water, plotted against forward calculations with the two biophysical models applied, as well physical dose and LET_d distributions (calculated via MC simulation). RBE variation as a function of depth with percent difference D_RBE (%∆_D-RBE) between planned and forward calculation D_RBE for MKM and LEM are displayed in the middle and bottom panel, respectively. The top axis of each figure segments regions of normal tissue (NT) and tumor tissue (T) for a representative clinical treatment

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_RBEVH) 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_RBEVH, 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, and ~ 10% for both LEM and the reference (D_RBE calculated with DDM). The D_5% for the rectum was 50.2 Gy (RBE) for MKM, 46.0 Gy (RBE) for LEM and 48.2 Gy (RBE) for DDM.

Fig. 3

D_RBE comparison illustrated in a clinical case (prostate cancer). a MC-optimized D_RBE distribution applying the DDM for the prostate case with resulting (b) LET_d distribution for Z = 2 particles. c Differences between the reference D_RBE (DDM) and LEM and MKM D_RBE are displayed in panels (c) and (d), respectively. Contours for the PTV, femurs and rectum are displayed in blue, green/yellow and red, respectively. DVH and D_RBEVH for the three biological models are depicted for the PTV and the rectum in the bottom left and right panel, respectively. Note that the critical organs at risk (e.g. anterior rectum) are susceptible to large variations in small volumes (< 5% of the relative total volume per organ) due to overlap with the tumor structure delineation. The asymmetric LET_d gradient (indicated by the yellow arrow) in panel (b) is indicative of the applied beam optimization procedure to meet dose constraints in the rectum while maintaining target coverage

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 LEM- and 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.

Fig. 4

Validation of the FRoG dose engine for helium ion beam therapy dose calculation with a glioma patient case. D_RBE applying (a) DDM and (b) LET_d is displayed, along with dose difference between D_RBE applying DDM and (c) LEM and (d) MKM. DVH and D_RBEVH for the three biological models are depicted for the PTV (right panel) for FRoG versus FLUKA. The yellow arrow directs attention to the LET_d gradient at the distal edge of the tumor, which could lead to larger uncertainty in RBE prediction for both the tumor and neighboring heathy issues beyond the target

As shown in Fig. 4, DVH and D_RBEVH 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_RBEVH 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.

Fig. 5

Clonogenic survival (top) of the Renca cells when exposed to various doses of a monoenergetic ⁴He ion beam at two measurements points. MC-estimated LET_d values are ~ 6 keV·μm⁻¹ at 6 mm depth (upper left) and ~ 15 keV·μm⁻¹ at 12 mm depth (upper right) using a ⁴He beam energy E = 56.65 MeV/u with a BP position (d_BP) of 25.9 mm. FLUKA-coupled biophysical and phenomenological models predicted cell survival and corresponding RBE (bottom) 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 (γ) with (α/β)_x = 1.79 Gy and ⁴He irradiation, respectively

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 to ~ 1.3, ~ 1.5, ~ 1.8 for the low, mid and high LET_d conditions.

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⁻¹, ~ 10 keV·μm⁻¹, ~ 15 keV·μm⁻¹, ~ 27 keV·μm⁻¹. 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 ⁴He irradiation, respectively. LQ-fit parameters for the four LET_d conditions are listed in the Table 2

Discussion

RBE model assessment

To best interpret the biological models for ⁴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 track-segment conditions, one can clearly discern the basis of intra- and intermodal variation as a function of dose, LET and tissue type.

Figure 7.a shows the comparison of RBE_α (top) and R_β (bottom), for mono-energetic ⁴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 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.

Fig. 7

Comparison between the three model predictions. (a) RBE_α (top) and R_β (bottom) as function of LET for (α/β)_x = 2 Gy (left) and 10 Gy (right). (b) RBE as a function of LET for (α/β)_x = 2 Gy (left panel) and 10 Gy (right panel) at 2 Gy reference photon dose. (c) RBE as a function of ⁴He ion beam physical dose for (α/β)_x = 2 Gy (top) and 10 Gy (bottom) at 4 keV·μm⁻¹ and 15 keV·μm⁻¹ as shown in the left and right column, respectively

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⁻¹ 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⁻¹ 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 ⁴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⁻¹ and ~ 1.45 for 15 keV·μm⁻¹ in low (α/β)_x tissues, and ~ 1.1 for 4.0 keV·μm⁻¹ and ~ 1.35 for 15 keV·μm⁻¹ 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⁻¹, R_β converges to ~ 0.6, while for 15 keV·μm⁻¹ 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⁻¹) and higher LET_d (≥10 keV·μm⁻¹) 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% and ~ 8%, 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 clinically-relevant 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 ⁴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 ⁴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 ⁴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 ⁴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 ⁴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 ⁴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/mechanistic- or phenomenological-based approach are essential prior to clinical translation of helium ion beam therapy.