FLUKA particle therapy tool for Monte Carlo independent calculation of scanned proton and carbon ion beam therapy

While Monte Carlo (MC) codes are considered as the gold standard for dosimetric calculations, the availability of user friendly MC codes suited for particle therapy is limited. Based on the FLUKA MC code and its graphical user interface (GUI) Flair, we developed an easy-to-use tool which enables simple and reliable simulations for particle therapy. In this paper we provide an overview of functionalities of the tool and with the presented clinical, proton and carbon ion therapy examples we demonstrate its reliability and the usability in the clinical environment and show its flexibility for research purposes. The first, easy-to-use FLUKA MC platform for particle therapy with GUI functionalities allows a user with a minimal effort and reduced knowledge about MC details to apply MC at their facility and is expected to enhance the popularity of the MC for both research and clinical quality assurance and commissioning purposes.

ticles are transported and implemented interaction models are limited. This restricts the usage of commercial, simplified MC TPS with respect to novel ion species, scoring of non-standard quantities (i.e. LET D , experimental and new RBE models) and accurate dose estimation in non-biological materials.
General purpose MC Codes have a full flexibility to support the above mentioned problems. They have demonstrated a potential as a tool for independent dose recalculation for patient specific QA, or for benchmarking for dose calculation algorithms of commercial TPS (IAEA and ICRU 2007, Parodi et al 2012, Grevillot et al 2012, Bauer et al 2014. Several institutions succeeded to develop their own customized solutions based on MC Codes, mainly for research purposes (Trento- (Fracchiolla et al 2015) (also clinical), CNAO- , HIT- (Parodi et al 2012)). Future systematic and user-friendly utilization of the general purpose MC tools for quality assurance (QA) and research, requires reliable and at the same time generally applicable interfaces to connect and embed MC tools in the radiation oncology workflow. One of the general purpose MC Code that has been continuously developed with a focus on medical application is FLUKA (Ferrari et al 2005, Böhlen et al 2014, Battistoni et al 2016, Augusto et al 2018.
The purpose of the current manuscript is to present first, an easy-to-use particle therapy platform based on FLUKA and to demonstrate its potential for clinical applications on selected examples for proton and carbon ion beams, i.e. beam modelling or patient specific treatment plan QA. Furthermore, with the new FLUKA release, the platform currently supports biological dose calculations as well as calculations of dose-averaged LET distributions for treatment plans. These new developments are scheduled for release in 2019 and will hence become widely available.

Materials and methods
Recent developments of the FLUKA particle therapy tool aim at supporting current research challenges. The first section provides general information about its functionalities. The subsequent sections present its usage for the following applications: • Patient specific treatment plan QA for proton therapy in clinical setting (section 2.3) • Biological dose scoring for carbon ion therapy (section 2.4) • Dose-averaged LET determination on the basis of treatment plans generated with a commercial TPS for carbon ion therapy (section 2.5)

FLUKA Monte Carlo code and its GUI Flair for particle therapy
The FLUKA particle therapy tool is based on the established general purpose MC code FLUKA, which includes the state-of-the-art physics models (Ferrari et al 2005, Böhlen et al 2014. FLUKA has been thoroughly benchmarked against depth-dose data and lateral profiles from research and clinical particle-beam therapy accelerators, for protons, carbon and various other light ions with therapeutic potential (Battistoni et al 2016, Parodi et al 2012, Mirandola et al 2015, as well as supported a beam characterization for the design of TULIP (TUrning LInac for Protontherapy) (Cuccagna et al 2018). To provide a flexibility in designing simulations for particle therapy as well as to improve the usability, FLUKA was enhanced with several new functionalities. The FLUKA particle therapy tool is integrated with Flair (Vlachoudis 2009, Vlachoudis and Sinuela-Pastor 2014)-a graphical user interface (GUI) and an integrated development environment (IDE). Flair enables creating and editing error-free FLUKA input files, writing and compiling user routines, executing simulations, data processing and plot generation. It is based on Python and thanks to its modular design the RTmodule for particle therapy was implemented to enhance the use of the new FLUKA functionalities for medical applications.
For particle therapy applications of FLUKA described in this study, the predefined simulation default settings-HADROTHErapy is always chosen to guarantee a reasonable simulation time and reliable accuracy. Particle transport thresholds are set down to 100 keV, except low energy neutrons, which are transported down to 10 −5 eV.
Atomic physics models handle continuous energy loss, energy loss straggling, delta-ray production (production cut at 100 keV) and the multiple Coulomb scattering of charged particles. Nuclear interaction models for hadrons, photons, muons and neutrinos are described using PEANUT model . For ions, in the range down to 0.1 GeV n −1 , FLUKA uses a modified version of relativistic quantum molecular dynamics (RQMD-2.4) model (Andersen et al 2004) implementation, below 150 MeV n −1 nuclear reactions are handled by Boltzman master equation (BME) , while smoothly transiting from one model to another in the overlapping range of energies.
Typical input and output standards used in the clinical experience is recognized by the platform, and the exported data readable by other software. The RTmodule uses the DICOM (digital imaging and communications in medicine (DICOM 2018)) standard for the data exchange with imaging tools or TPS. More specifically, it can process computer tomography (CT) files into the voxelized patient geometry and translate and embed regions of interest (ROI) from DICOM RTSTRUCT file into patient-voxel files readable by FLUKA. The module imports information from each pencil beam spot described in the DICOM RTPLAN file. It processes the DICOM RTDOSE file into the FLUKA format and provides a tool for direct comparison with the FLUKA output. The obtained FLUKA data can be then exported to the DICOM RTDOSE format. For treatment plan calculations, the FLUKA particle therapy tool enables handling a multiple-field treatment at once or performing individual treatment field analysis with possibility of merging data from individual fields. Furthermore, for the sake of the compatibility with clinically used and commissioned TPS, on-the-fly conversion of dose-to-medium to dose-to-water can be performed for both physical and biological calculations (Bauer et al 2014).
The standard patient model for stoichiometric calibration of voxelized CT-based geometry is based on a published Hounsfield unit (HU) calibration procedure (Schneider et al 2000). Creation of a voxel geometry requires two files. The individual assignment of HU (or any other arbitrary unit) to tissue material, and the composition of tissue components. It is important to emphasize that if the individual stopping power (dE/dx) determined by FLUKA does not match exactly the stopping power used for clinical TPS (due to minor differences in the defined material composition), there will be a systematic deviation in the particle ranges. Correction factors for stopping powers can be added to overwrite FLUKA-calculated values.
The beam source can be generated directly from clinical treatment plans. The DICOM RTPLAN file contains the essential information of each pencil beam spot (i.e. energy, position, rotation), which are then translated to FLUKA-proprietary functions. Specific beam model information is treated as supplementary data; details for selected key energies needs to be provided in a separate file, and the values between key energies are interpolated using a spline algorithm. This supplementary information encompasses: beam size, angular spread, momentum spread and other data specific to the treatment and gantry/beamline facility, such as the normalization of monitor units (MU) to the number of primary particles. In the other words, each defined pencil beam spot has its assigned individual set of modifiable parameters. This eases and accelerates the facility-specific beam model creation and verification.
While importing DICOM information into the FLUKA input RTmodule, scoring parameters (such as: voxel size, dose binning and patient/image orientation) are automatically determined to simplify subsequent dosimetric assessment. The RTmodule determines the scoring resolution based on the DICOM RTDOSE file, however, on top of the standard dose scoring, it is possible to modify the pre-generated settings or to specify additional scoring estimators.
Finally, for biological dose (D RBE ) computation as presented here, a radiobiological database for the local effect model (LEM) (Scholz et al 1997) is used for FLUKA simulations (Mairani et al 2010). In the respective FLUKA input file it is possible to include external files that describe alpha and beta parameters within the linear quadratic formalism, of different tissue components of the mixed radiation field as a function of energy per nucleon and ion. This allows to interface easily external RBE models with the FLUKA computation. At the end of the simulation, doseweighted averages of α mixed D and β mixed D are obtained for each voxel using the formalism from (Krämer and Scholz 2006), followed by the calculation of the survival curve and dose-weighted RBE (see appendix).

Commissioning of the FLUKA particle therapy tool for independent dose calculations
For testing and illustrating the versatility of the FLUKA particle therapy tool on clinical cases, the tool was first commissioned for two particle therapy facilities, i.e. the synchrotron based facility CNAO for proton and carbon ion therapy and the cyclotron based proton facility at Trento. For these two facilities two different levels of commissioning were performed, which is described briefly in the following. Then the commissioned data were handled via the RTmodule functionalities.
Commissioning for the CNAO facility: The geometry of the CNAO horizontal beam-lines was modelled, including the vacuum window, the beam monitoring system and the air gaps. For each proton and carbon ion beam within the clinically relevant energy range, both the momentum spread and the size of the beam before the nozzle were adjusted to match the commissioned measured values. Furthermore, for carbon ion beams, ripple filters were embedded into the geometrical description as additional 'lattice' layers. A previously determined facility-specific CT calibration curve was provided. FLUKA-based stopping powers for different materials were slightly adapted to match the predetermined TPS-based dE/dx values.
Commissioning for the Trento facility: For the modelling of the Trento proton facility, no detailed information concerning the beam delivery system was available, thus the methodology presented in Fracchiolla et al (2015) was used. Original beam data used for the MC model in Fracchiolla et al (2015) was retrieved and described in the supplementary beam model file. Simulations for 16 key energies (between 70 MeV and 200 MeV in steps of 10 MeV) were performed and compared with commissioning measurements in terms of range and spot size. The range was defined as the position of the 90% of the maximum of the curve in the Bragg peak's distal falloff (R 90 ) for integral depth dose curves measured in water, while spot size was defined as sigma obtained for a single Gaussian fit in transversal planes for measurements in air at five distances from the isocenter (−19.8 cm, −10 cm, 0 cm, +10 cm, +20 cm). For treatment simulations, the beam energy delivered for the Trento cyclotron was assumed to be continuous in the range from 70 MeV to 226 MeV. Values for the energies in between were interpolated using a spline algorithm. For the commissioning of the facility, a specific CT calibration curve and two phantoms were used to adjust the default HU-tissue calibration within the FLUKA particle therapy tool. The first phantom consisted of two material blocks with significantly different densities (two plates substituting tissues-lung 0.3 g cm −3 and bone 1.819 g cm −3 ), while the second one was an anthropomorphic lung phantom.

Application 1: patient specific treatment plan QA in particle therapy
Representative clinical treatment plans were selected at the two different particle therapy treatment centers. These were generated using commercial and fully commissioned TPS (Syngo at CNAO, RayStation version 6i in Trento). At both sites dose distributions resulting from treatment plans are calculated with pencil beam algorithms implemented in the TPS.
Based on this data, FLUKA input files were built using the RTmodule. Pencil beam scanning data (particle species, position, rotation, intensity and spot size at the isocenter) were extracted from treatment plans. For both cases supplementary beam model file with the additional, commissioned beam source information were used. In the case of CNAO simulations, the model of the beam nozzle was designed using FLUKA combinatorial geometry, however the model, once created, can be easily re-adapted for further simulations.
Two proton chordoma patient cases were selected from each of the facilities. The cases are representative for rather homogeneous tissue conditions. Each treatment plan consisted of three fields. As third proton example, a head-and-neck cancer case, representative for more heterogeneous tissue conditions, was selected. In addition, in this setting a 4 cm thick range shifter was used and modeled in FLUKA geometry. The last clinical case was a chordoma patient, treated with carbon ions using a two field arrangement.
Standard dose scoring in terms of dose-to-water was applied, as well as dose-to-medium. For proton simulations 1% of the planned number of particles was used, while for the carbon ion 10% of planned particles was used to obtain error below 1% for values above 25% of the maximum dose. All calculated dose distributions were evaluated in terms of Dose-Volume-Histograms and several DVH parameters were extracted and compared (D 5 and D 95 for PTV, D 5 , D 7.5 and D 10 for OAR). As in clinical practice a fixed proton RBE of 1.1 was used.

Research application 2: biological dose calculation
For biological dose calculation, a carbon ion treatment plan of a chordoma patient case was selected. Biological dose calculations with FLUKA were performed using the LEM I model with the parameter set, given in table 1.

Research application 3: dose-averaged linear energy transfer scoring
The final research application presents the possibility of extracting an unrestricted dose-averaged linear energy transfer (LET d ) distribution from TPS files, which can serve as a simplified indicator for future approaches of treatment plan optimization or radiobiology oriented research.
The scoring of the LET d (ICRU 1970) was performed with external FLUKA user routines, which will be added as a standardized estimator in the future. This external routine was designed to derive the total LET d , which takes into account the contribution from primary and secondary particles from each beam spot, calculated through scoring separately the numerator and denominator of the equation (1) from (Wilkens and Oelfke 2003): where S i is the is the stopping power of the beam spot i with energy E at depth z, and ϕ i is a local particle spectrum of the beam spot i with energy E at depth z.

Commissioning
Since the commissioning of the CNAO facility with FLUKA was previously performed and the related results can be found in Molinelli et al (2013), Mirandola et al (2015), Parodi et al (2013) and Mairani et al (2013), we retrieved the required beam and calibration data. This information were directly used in FLUKA particle therapy tool for the presented CNAO clinical cases. However, the commissioning of the Trento beam model required model validation in FLUKA. Simulations of pencil beams with selected energies showed a very good agreement with integrated depth doses (IDD) measured in water (figure 1) and Bragg peak shape. The spot size differences in air for x and y axis at the isocenter between measurements and simulations are less than 10% (figure 2), while for all five distances, the mean difference does also not exceed 10% for any of the key energies. This difference is considered to be adequate from a clinical perspective because it is comparable with the spot size fluctuations for each energy as a function of spot position and gantry angle (Schwarz et al 2016). Concerning the CT calibration curve, the discrepancies between the dose distributions obtained with FLUKA and the clinically used TPS were negligible, therefore it was decided not to correct and overwrite stopping powers calculated by MC. A more detailed approach would require calibration  with a tissue characterization phantom, composed of several materials of known elemental composition and electron density.
Based on the commissioned data we tested the functionality and verified our automatized approach on the proton and carbon case presented in this article.

Application 1: patient specific treatment plan QA in particle therapy
Chordoma cases (three field arrangements): figure 3 shows the comparison of dose calculation results (at the isocenter plane) obtained with the respective clinical TPS (a), with the FLUKA particle therapy tool (b) for the clinical chordoma case at CNAO (ChordomaCN), as well as the absolute difference (c). As expected, for this relatively homogeneous case, a good consistency was found for dose-to-water calculation. Figure 4 shows the corresponding DVH, and table 2 summarizes respective DVH parameters. For this chordoma case at CNAO the difference in selected DVH parameters for the PTV did not exceed 0.4% and for OAR the maximum difference of DVH parameters was 3%.
Similar results were obtained for the chordoma case at Trento proton facility (ChordomaTR), although these results were found to be slightly less accurate, with a difference in the DVH parameters up to 1.6% for the PTV and 4% for OAR.
Head and neck (2 fields arrangement): figure 5 shows the comparison of dose calculation results obtained with the respective clinical TPS and those obtained with the FLUKA particle therapy tool for the head and neck case at Trento proton facility, in the principal planes ((a)-(c) axial, (d)-(f) coronal, (g)-(l) sagittal) covering the isocenter. The right column figures show the absolute dose difference at the isocenter. The TPS overestimated the dose in empty cavities up to 16% and thus underestimated the particle ranges. The resulting differences in the DVH are shown in figure 6. Table 3 summarized the most essential DVH parameters. D 95% exceed 14% for PTV1 and D 95% 21% for PTV2. For OAR these values were respectively: 7% for D 5% at brainstem, 22% for D 10% at optic nerve, 8% for D 10% at thyroid.   These significant differences were first associated with the range shifter modelling in TPS. To test this hypothesis, patient specific QA measurements were retrieved from IBA MAtriXX and compared with MC simulations performed at three depths in a homogeneous Gammex phantom: at 2 cm, 3 cm and 5 cm. The calculated doses were subsequently compared to measurements and a gamma index analysis was performed with a dose tolerance (DT) of 3%, distance-to-agreement (DTA) of 3 mm and lower dose threshold of 1% of the maximum measured dose. In case the gamma evaluation is performed with MC simulations, potential artifacts, due to statistical noise, need to be considered (Low and Dempsey 2003). The average score for gamma index analysis at all three depths and for both fields was 96%. Consequently the range shifter model used in MC was ranked as acceptable and the dose differences between MC and TPS calculations were associated rather with shortcomings of pencil beam algorithm modelling for the primary beam interaction in the range shifter and lateral inhomogeneities (Saini et al 2017, Widesott et al 2018.

Research application 2: biological dose scoring in ion therapy
The biological dose scoring for the chordoma case treated with carbon ion therapy at CNAO is presented in figure 7. The right hand figure is the physical dose distribution after recalculation with the FLUKA particle therapy tool, while the left hand side figure shows the resulting biological dose. As mentioned previously, for the scope of this comparison the LEM I was used, as it is the standard clinical radiobiological model employed for carbon ion therapy in Europe. The scored RBE-weighted dose was compared with the results obtained in the TPS, and good agreement is observed in 2D cross-section presented in figure 7 as well as for the DVH comparison shown in figure 8. Figure 9 shows the dose-averaged LET distribution (LET D ) in the isocenter plane of the carbon ion case, and the resulting LET D -volume histograms (LVH). When studying the LVH the pronounced variation of the LET D across the target and OAR can be noticed. While 95% of the PTV received values of 55 keV µm −1 , for the 'most exposed 5%' of the PTV volume the LET D almost doubles to 102 keV µm −1 . The highest LET D values above 100 keV µm −1 were found to be located at the PTV edges. For the selected OARs, i.e. brainstem and optic nerve, 10% of the respective volumes received LET D of at least 67 keV µm −1 and 47.5 keV µm −1 .

Discussion
Particle therapy is a growing field within radiation oncology that faces a number of challenges. Starting from the beam delivery, passive scattering is being eclipsed by pencil beam scanning technology, for which much less experience has been gathered so far (Lomax 2008a(Lomax , 2008b. The biological quantification of particle beams and the assessment of treatment plans are re-discussed, this includes the clinically used constant RBE for proton beams, as well as radiobiological models in general. In addition there is and ongoing discussion about novel ions species beyond proton and carbon ions (Mairani et al 2016. This manuscript illustrates the potential of the FLUKA particle therapy tool in this context, ranging from QA purposes in clinical particle therapy environment to clinical radiobiology oriented research.
The FLUKA particle therapy tool has been specifically designed to support independent dose verification for particle therapy, as a recommended part of the QA procedure (IAEA and ICRU 2007). For fluence modulated radiotherapy techniques, for both photon and particle therapy, the current best practice is to perform experimental treatment plan verification as part of the patient specific QA program. This procedure is beam time and workload intensive, and inefficient in detecting errors (Lomax et al 2004, Arjomandy et al 2010, Furukawa et al 2013. Thus, independent dose calculation (supported by log-file recalculations) is an alternative procedure for patient specific QA, but also for systematic investigations of dose calculation uncertainties or uncertainties in the beam model of commercial TPS based on studies using (anthropomorphic) phantom or patient cases. Modelling of the primary beam interaction within the range shifter or the poorly modelled multiple Coulomb scattering in complex geometries were identified as shortcomings of dose calculation algorithms in commercial TPS (Tourovsky et al 2005, Paganetti 2012, Verburg and Seco 2013, Taylor et al 2017, Saini et al 2017, Widesott et al 2018. When aiming to apply general-purpose MC codes for verifying TPS, and its integration into the clinical research, workflow is key-implying smooth and flexible data import from the TPS and imaging systems. With the proton and carbon ion therapy cases exemplified above, the FLUKA particle therapy tool demonstrated such functionality and usability in a clinical setting. The current implementation allows fast export of the TPS radiotherapy plan files, which was tested on two different TPSs. The Flair geometry editor allows adapting the simulations into the specific treatment scenarios, and if needed supports modelling of the range shifter, ripple filter or other devices used during treatment. Finally it is possible to export the results into DICOM standard files, readable by most clinical and research software.  A limitation of some of the above-presented clinical cases, lays mostly in the not fully exploited commissioning measurements and applied simplifications. Obviously commissioning and beam modelling of the FLUKA particle tool itself has an impact on the achievable dose calculation accuracy. For example, the differences in the chordoma case from the Trento proton facility showed slightly higher differences compared to the CNAO case. This was inevitable due to less rigorous beam modelling as well as the HU-tissue calibration, which was on purpose not exploited fully to illustrate that the FLUKA particle tool can be setup easily, with acceptable results, even if the detailed geometry of the beam delivery system or the HU calibration curve is not known. Consequently, this affected insignificantly the ranges of the primary beam spots, resulting in marginally higher differences in the selected DVH for the Trento cases. On the other hand, HU calibration has also an inherent uncertainty around 2%-3% (Schuemann et al 2014). However, even with a rather generic HU calibration curve and no detailed beam line information the FLUKA particle therapy tool can be used for clinically motivated and oriented research.
One of these applications is the systematic exploration of helium or oxygen ions, which FLUKA can support due to its benchmarked interaction models. Especially helium ions have become a research focus during the last years and it is generally agreed that they are a promising alternative to proton therapy due to their favorable physical and biological properties (Fuchs et al 2015, Mairani et al 2016, Knäusl et al 2016, Tessonnier et al 2017a, 2017b. The versatility of FLUKA for helium ion beam research was recently demonstrated for brain and ocular meningiomas (Tessonnier et al 2018). Despite the promising result achieved so far, further studies for different tumor entities are still needed to investigate the clinical potential of helium ions.
Moreover, the radiobiological modelling of helium remains a major challenge in the light of the limited number of experimental data. Based on the pioneering work at the Heidelberg ion therapy facility it can be expected that other European synchrotron based particle therapy facilities will follow and more radiobiological data will become available in the coming years. Thorough benchmarking of the applied FLUKA model for RBE-weighted dose calculations based on the LEM model was already demonstrated for protons and carbon ions (Bauer et al 2014, Mairani et al 2008. The FLUKA interface for biological calculation based on the linear-quadratic (LQ) formalism allows to apply biological calculations of the plan with several different models, such as LEM IV (Elsässer et al 2010) and microdosimetric kinetic model (MKM) used by NIRS (National Institute for Radiological Sciences, Japan) (Magro et al 2017). Thus the FLUKA particle therapy tool can further serve as support for compariso n between different biological models and for exploring novel ions species.
The uncertainties in radiobiological models motivated to explore the dose average LET (LET D ) as additional descriptor for characterizing dose distributions in particle beam therapy. It is well known that high-LET results in higher cell-killing per absorbed dose, and can overcome radio resistance (Schlaff et al 2014) of tumor cells, but it is also well known that high-LET causes higher toxicity when delivering doses with high RBE to OAR. The significant difference between the RBE-weighted DVH and LVH exemplified above in section 3.4 underlines the importance of LET assessment (Cao et al 2018). In the new clinical era of particle beam therapy that has been initiated with the widespread implementation of pencil beam scanning, such additional scoring of treatment plans becomes essential. The FLUKA particle therapy tool can support such clinical research related to LET D or can be used to benchmark or validate LET calculation modules in commercial TPS.

Conclusion
In this study the performance of the FLUKA particle therapy tool in the clinical environment was presented. The application was tested in the scope of two different facilities and against their commissioned TPS for QA purposes, RBE-weighted dose and dose-averaged LET scoring, for proton and ion treatment plans. Presented work is encouraging and it is planned to be applied in the clinical environment in the near future.
To conclude, the FLUKA MC code using well benchmarked physics models provides good accuracy for clinical purposes, and the newly developed RTmodule in GUI Flair supports user during the simulation workflow. FLUKA particle therapy tool is reliable and well integrated with the code itself, it allows to import data, tailor it according to the user's needs and arrange the simulation input file with no detailed knowledge of the MC code. The user can share the results and export the output for further analysis in external software. Due to its flexibility it is easily applicable also for research purposes. Current work is focused on improving the usability of the particle therapy tool and extending the range of applications for the clinical and research environment.