The FLUKA Monte Carlo code coupled with an OER model for biologically weighted dose calculations in proton therapy of hypoxic tumors

cells compared to well-oxygenated cells is quantified by the oxygen enhancement ratio (OER). In this study we created a FLUKA Monte Carlo based tool for inclusion of both OER and relative biological effectiveness (RBE) in biologically weighted dose (ROWD) calculations in proton therapy and applied this to explore the impact of hypoxia. Methods: The RBE-weighted dose was adapted for hypoxia by making RBE model parameters dependent on the OER, in addition to the linear energy transfer (LET). The OER depends on the partial oxygen pressure (pO 2 ) and LET. To demonstrate model performance, calculations were done with spread-out Bragg peaks (SOBP) in water phantoms with pO 2 ranging from strongly hypoxic to normoxic (0.01–30 mmHg) and with a head and neck cancer proton plan optimized with an RBE of 1.1 and pO 2 estimated voxel-by-voxel using [ 18 F]-EF5 PET. An RBE of 1.1 and the Rørvik RBE model were used for the ROWD calculations. Results: The SOBP in water had decreasing ROWD with decreasing pO 2 . In the plans accounting for oxygenation, the median target doses were approximately a factor 1.1 lower than the corresponding plans which did not consider the OER. Hypoxia adapted target ROWDs were considerably more heterogeneous than the RBE 1.1 - weighted doses. Conclusion: We realized a Monte Carlo based tool for calculating the ROWD. Read-in of patient pO 2 and estimation of ROWD with flexibility in choice of RBE model was achieved, giving a tool that may be useful in future clinical applications of hypoxia-guided particle therapy.

applied to transfer clinical treatment protocols from photon therapy to proton therapy. This is currently implemented by delivering a homogeneous RBE-weighted dose to the target, with the RBE set to 1.1 [3]. However, the RBE is known to vary and several RBE models at normoxic conditions have been proposed [4][5][6][7]. To account for the varying oxygen levels and the variable RBE in particle therapy, a common framework for combing the effects in RBE and OER weighted dose (ROWD) calculations could improve the patient treatment outcome.
A common method of measuring the pO 2 in vivo is by invasive oxygen-sensitive electrodes. However, such methods are limited to accessible tissues and does not give a full three-dimensional representation of the pO 2 levels in the patient [8]. With functional imaging such as positron emission tomography (PET) or magnetic resonance imaging, information about the oxygen level can be obtained on a voxel-by-voxel basis [9]. Several methods for depicting hypoxia with functional imaging have been proposed, and PET is currently the preferred modality when studying hypoxia in clinical settings [10]. Among the tracers in clinical use, two hypoxia radiotracers for PET are 18 F-fluoromisonisazole ([ 18 F]-FMISO) and 18 F-labeled 2-(2-nitro-1 H-imidazol-1yl)-N-(2,2,3,3,3-pentafluoropropyl)-acetamide ([ 18 F]-EF5). While [ 18 F]-FMISO is the most actively investigated hypoxia PET tracer, [ 18 F]-EF5 has the advantage that the non-labeled form of this compound has been widely used to study hypoxia in clinical and preclinical settings. Moreover, the hypoxia-avidity of this tracer has been demonstrated convincingly [11][12][13]. Also, EF5 binding in hypoxic head-and-neck squamous cell cancer is rather well established and has been shown to relate with poor outcome in photon radiochemotherapy [14,15].
Monte Carlo (MC) simulations are frequently used to provide estimates of the RBE-weighted dose in proton and carbon ion therapy, and several RBE models which do not account for hypoxia have already been implemented in the FLUKA MC code [16,17]. This includes the local effect model (LEM) [18,19] and the microdosimetric kinetic model (MKM) [20,21], whereas several methods for accounting for both RBE and hypoxia in particle therapy treatment planning have been suggested [22][23][24][25]. Scifoni et al. [23] and Tinganelli et al. [25] described a ROWD optimization algorithm incorporating OER (through a semiempirical model) and RBE (using LEM IV), and implemented this into the analytical treatment planning system (TPS) TRiP98 [26,27]. They applied an OER model which was intended for particle therapy with heavier ions like carbon and oxygen. For protons, the same ROWD estimation method can be applied, however, an OER model developed specifically for protons would be more suitable, as the universal OER models may overestimate the OER for protons [2]. Wenzl and Wilkens [28] developed an OER model with parameters which depend on the LET and the pO 2 , which was based on experimental data from several types of particles, including, but not limited to, protons and carbon ions. To our knowledge, none of the currently available OER models have parameters based on proton data only, although Strigari et al. [24] included a particle dependency in their calculation of the OER. In this study we developed a method for recalculating treatment plans with a ROWD, with flexibility regarding the choice of RBE model, and implemented the method in the FLUKA MC code. This will indicate potential areas of lower tumor dose effect in a volume with varying degree of hypoxia. Including the method in the dose optimization process would then enable delivering a homogeneous ROWD, i.e. equal biological effect across regions with different oxygen levels.

Estimation of oxygen levels in a patient
The patient oxygen levels, given as pO 2 , were estimated using PET images of head and neck cancer patients with [ 18 F]-EF5 as hypoxia tracer. The PET images were acquired between 2013 and 2016 using a GE D690 PET/CT scanner (General Electric Medical Systems, Milwaukee, WI, USA) at Turku University Hospital (Finland) [29]. The relationship between PET uptake and pO 2 was estimated applying the following equation: where A, B and C are reaction-specific parameters [30]. We estimated these parameters for [ 18 F]-EF5 PET by non-linear least square curve fit of the uptake data from PET images of six head and neck cancer patients in organs with pO 2 values reported in literature, see Appendix A for details. The PET uptake was normalized to mean muscle uptake, as a tumor-to-muscle (T/M) ratio of 1.5 at 3 h from injection has been found to be an appropriate threshold for determining clinical significant hypoxia in [ 18 F]-EF5 PET images [13].

Calculations of oxygen enhancement ratio
The OER was calculated applying a version of the model by Wenzl and Wilkens [28], modified by obtaining the model coefficients by regression to in vitro proton data only. The Wenzl and Wilkens [28] OER model calculates the survival fraction of cells, S, after ion irradiation according to the linear quadratic (LQ) model [31]: where D is the physical dose. and are the radiosensitivity parameters, defined by Wenzl and Wilkens [28] as follows: where L is the dose-averaged LET (LET d ), p is the pO 2 , and K is set to 3 mmHg [28] . By non-linear least square curve fit of in vitro proton data ( Fig. 1 and Table B where p h and p a are the oxygen pressure at the hypoxic and aerobic (or normoxic) conditions, respectively. In our implementation p h was obtained from the PET image and we set p a equal to 30 mmHg as this is a typical pO 2 value for tissues [28]. This value for p a is lower than what is given in Fig. 1. Fig. 1 applies the aerobic condition in in vitro experiments, which is generally fixed to 160 mmHg and higher than what is normally found in vivo [28]. As the model parameters are based on a survival fraction of 10%, we also set the survival level to 10% in our calculations.

RBE-weighted dose calculations for normoxic conditions
To adapt the RBE-weighted dose for hypoxia, a short description of RBE calculations for normoxic conditions is needed. For a more thorough description, see Rørvik et al. [32]. In general, the RBE is calculated from the LQ model as follows: where D p is the physical proton dose deposited per fraction and α a , β a , α x and β x are the normoxic proton and photon radiosensitivity parameters, respectively. The various RBE models differ in how these radiosensitivity parameters are defined [32]. We applied two methods to calculate the RBE, as described below.
In the first method, the RBE-weighted dose (D RBE1.1 ) was calculated using a constant RBE of 1.1. To achieve this using Eq. (6), which is required to later include the OER in the calculations, we used the same radiosensitivity parameters for photons and protons, i.e. = a x and = a x . This would give the same survival curve for protons and photons, and thus an RBE of 1. The D RBE1.1 is then given by Eq. (6) multiplied by 1.1 and total physical dose, D: where = a x and = a x , and are given as in the previously described OER model (Eqs. (3) and (4)), with a pO 2 value set to 30 mmHg and a constant LET d of 0.3 keV/µm.
In the second method, the RBE-weighted dose (D ROR ) was estimated applying the variable RBE model by Rørvik et al. [5] (see Appendix C for details). We applied an ( / ) x = 2 Gy for all structures. However, published ( / ) x ratios for head and neck cancer vary significantly, and when applying the tool for further testing, a larger range of ( / ) x ratios should be applied [33]. As the goal here was to demonstrate the method, we chose a ( / ) x ratio representing organs at risk (OARs) and late responding target volumes which have low ( / ) x values [32].

RBE and OER weighted dose calculations
To account for hypoxia in the calculations of ROWD, we have, similar to Scifoni et al. [23] and Tinganelli et al. [25], replaced the normoxic proton response parameters in Eq. (6) with hypoxic proton response parameters which vary with the OER: where the OER is the OER at 10% cell survival calculated according to Eq. (5). This ROWD model allows us to apply a , a , x and x as defined in most of the existing RBE models derived from normoxic proton and photon in vitro data and then estimate h and h applying Eqs. (8) and (9). The ROWD accounting for hypoxia is then given as follows: where D is the total deposited physical dose and D p is the physical dose deposited by protons only. While this ROWD model will give an isoeffective dose accounting for hypoxia, we called it D OER,RBE to emphasize that it is not directly the RBE-weighted dose, as only the proton parameters take hypoxia into account and not the photon parameters. The x and x parameters are always the normoxic photon radiosensitivity parameters.
The D RBE1.1 and the D ROR models were adapted to variable oxygen levels by changing their respective a and a parameters to h and h according to Eqs. (8) and (9), and keeping x and x as defined in each model. To distinguish the models adapted for hypoxia and not adapted for hypoxia, the models adapted for hypoxia will hereafter be referred to as D OER,RBE1.1 and D OER,ROR .

Implementation of the method in FLUKA
Our in-house made system [34] for recalculating treatment plans in FLUKA was modified to include the hypoxia inclusive ROWD calculation method. This tool translates DICOM information from a treatment plan generated in an analytical TPS to files readable by FLUKA. This includes information on the radiation beam and scrips translating the CT image into FLUKA geometry. When simulating a treatment plan, we scored LET d , pO 2 , physical dose D, D a and D a in FLUKA. The OER was then calculated in Python according to Eq. (5), using LET d , pO 2 and D from FLUKA, and the ROWD was subsequently calculated using Eq. (10).
LET d , D a , D a and pO 2 were scored using the FLUKA subroutine fluscw (FLUence Scoring Weight). LET d was here estimated applying the same approach as in Fjaera et al. [34], i.e. using the following equation: where E ( ) i is the fluence of particle species i with kinetic energy E and E LET ( ) i is the unrestricted LET for the same kinetic energy. To score pO 2 , we implemented a method for importing estimated pO 2 values in a patient geometry on a voxel-by-voxel basis in FLUKA using fluscw. This subroutine reads the pO 2 values from a table, created from the PET images; the coordinate system of the PET image was converted into FLUKA coordinates, and for each particle position, the corresponding pO 2 value from the PET image was given. Subsequently, pO 2 times dose to water and dose to water were scored in FLUKA. Then, pO 2 was given by the quotient of the two scored quantities. To achieve the needed accuracy of the scored values without increasing the simulation time significantly, we set the fraction of the kinetic energy to be lost in a simulation step to 0.0125, which was 25% of the default value when applying the FLUKA HADROTHErapy defaults.

Treatment plans
Two spread-out Bragg peak (SOBP) scenarios in water and a treatment plan for head and neck cancer were used to test the MC tool and investigate the impact of hypoxia on ROWD estimations. The treatment plans were optimized using an RBE of 1.1 in the Eclipse (Varian Medical Systems, Palo Alto, CA, USA) TPS, not accounting for hypoxia, and then recalculated to account for hypoxia and RBE using our FLUKA based tool. For all simulations, the number of simulated primary particles was chosen to have a voxel mean statistical uncertainty < 2%.
The first SOBP scenario had the same dimensions as the pO 2 table, with scoring voxels of (1.95 × 1.95 × 2.4 mm 3 ). This SOBP scenario was only used to check the FLUKA estimation of pO 2 values. The second SOBP scenario (starting at depth 8 cm and with size 4 × 4 × 4 cm 3 ) was optimized to an D RBE1.1 dose of 2 Gy(RBE) and applied in four simulated water phantoms. Each simulated water phantom had homogeneously distributed pO 2 values, with each of the four phantoms having a given pO 2 : 0.01 mmHg, 0.5 mmHg, 3.0 mmHg or 30.0 mmHg. The simulated phantoms had scoring voxels of same size as in the TPS (2 × 2 × 2 mm 3 ).
The head and neck cancer patient was originally treated at Turku University Hospital with chemoradiotherapy [29], and the PET images of this patient were included in creating the PET uptake to pO 2 conversion curve (Appendix A). For this study, we created a proton treatment plan for the patient, optimized (using an RBE of 1.1) in the Eclipse TPS with prescribed dose of 70 Gy(RBE) in 35 fractions to the planning target volume (PTV). pO 2 values of the patient were estimated from [ 18 F]-EF5 PET/CT images of the patient, and the FLUKA scoring was done using the same grid and voxel size (1.95 × 1.95 × 2.50 mm 3 ) as used during the TPS optimization.

FLUKA readout of pO 2
The generated pO 2 map was successfully read by FLUKA, as shown in Fig. D.1 in Appendix D. The process of importing patient pO 2 data in FLUKA started with a [ 18 F]-EF5 PET image (Fig. D.1a). The conversion from PET uptake to pO 2 was done applying Eq. (1) with reaction-specific parameters found to be A = 2.60, B = 1.98 and C = 2.5 mmHg, see Appendix A for details. The estimated pO 2 table (Fig. D.1b) corresponding to the PET image showed areas of low pO 2 values in the PTV, with minimum, mean and maximum pO 2 in the PTV of 2.2, 16.8 and 60.0 mmHg, respectively. This was also reproduced by the FLUKA imported pO 2 values (Fig. D.1c).

Water phantom
Compared to D RBE1.1 , the ROWD accounting for hypoxia decreased with decreasing pO 2 values, as seen in Fig. 2. The decrease rate typically followed the factors corresponding to the OER values for respective oxygenation levels in Fig. 2. This shows that the proton irradiation is predicted to be less effective at low oxygen levels, as expected from the definition of the OER. At 30 mmHg, which is a normoxic condition, the D OER,RBE1.1 dose was equal to the prescribed dose, as expected. Applying the D OER,ROR model, however, an RBE of 1.09 was observed at the beam entrance, increasing to 1.38 at the distal end of the beam, at normoxic condition. For all oxygen levels, the OER value was seen to be mostly constant along the beam path, with a slight decrease at high LET d values. This decrease was most significant at low oxygen levels. At 0.01 mmHg, which is an extremely hypoxic condition, the OER value started at 2.7 for low LET d values and decreased to 2.3 at the end of the beam range.

Patient case
The pO 2 map of the head and neck cancer patient is shown in Fig. 3a. The treatment plan had LET d values in the range of 1.2-4.3 keV/µm in the PTV (Fig. 3b). The dose volume histograms (DVHs) of the D RBE1.1 , D ROR , D OER,RBE1.1 and D OER,ROR are shown in Fig. 3c. The D ROR model gave the largest dose to the PTV, and the D RBE1.1 the second largest dose, except in 18% of the volume where the D OER,ROR had larger dose, as D OER,ROR gives the same dose as the D ROR model in the normoxic parts of the PTV. The normoxic RBE-weighted doses (D RBE1.1 and D ROR ) had the steepest DVHs, which was also reflected in the dose distributions shown in Fig. 3d-f, where a homogenous D RBE1.1 to the PTV was observed, while the ROWDs accounting for hypoxia had areas with clearly lower dose. The corresponding D ROR distribution is given in Fig. 4, showing also a mainly homogeneous dose distribution to the PTV.  Fig. 3g, showing even clearer that low oxygen levels result in lower biological effect of the radiation, as expected. This is also seen in Fig. 4, where the D ROR -D OER,ROR dose difference is given, showing that the differences D ROR -D OER,ROR and D RBE1.1 -D OER,RBE1.1 are comparable. The similar shapes of the dose difference can be attributed to the approximately homogeneous LET d and physical dose distribution in the target.

Discussion
In this study we developed a tool for calculating the ROWD in hypoxic conditions for protons, by coupling the FLUKA MC code with a ROWD calculation model. The tool was tested with simple treatment plans for homogeneous water phantoms with several oxygen levels as well as on a head and neck cancer patient plan with pO 2 estimated from [ 18 F]-EF5 PET images. Applying a method based on the LQ model has several benefits, one being that the RBE effects are also directly included in the calculations. Such calculations may become useful in future clinical applications of hypoxia-guided particle therapy.
Scifoni et al. [23] and Tinganelli et al. [25] accounted for hypoxia applying the same approach as in the current work, by using the LEM and a dedicated OER model implemented in the TRiP98 TPS, which facilitated RBE and OER driven optimization (kill painting). Wouters and Brown [35], and later Malinen et al. [36], also applied this formalism to estimate hypoxic α and β values from normoxic α and β  parameters. Using this approach, the physical dose is weighted by the inverse of the OER before calculating the RBE. Due to the dose dependency of a variable RBE, the RBE will therefore be slightly dependent on the OER. The estimated ROWD is therefore not equal to the RBE-weighted dose divided by the OER, but the difference is small. We applied our tool using the variable RBE model by Rørvik et al. [5], in addition to a constant RBE of 1.1. The tool could also be applied on most other normoxic proton RBE models, as well as RBE models for heavier ions, which follow the general RBE equation (Equation (6)). The OER model used in this study was based on the model by Wenzl and Wilkens [28], but with model parameters estimated from proton in vitro data only, as the OER has been seen to decrease significantly faster with the LET d when applying low-LET irradiation compared to high-LET irradiation, especially towards the end of the particle range [2]. Strigari et al. [24] observed this same trend, i.e. the lighter the particles the faster the OER decreases with increasing LET, when computing OER vs. LET curves for different particle types. This was also observed in our study, as the OER decreased faster when applying proton parameters than when applying parameters based on the data for several ion types ( Fig. 1). When applying the tool to high-LET radiation, e.g. carbon ions, the OER model parameters must be adjusted to data taken from heavier ions, for instance using the parameters given by Wenzl and Wilkens [28].
The Wenzl and Wilkens [28] OER model has been seen to deviate from other existing OER models, resulting in lower OER for low LET radiation compared to the models by Scifoni et al. [23], Strigari et al. [24] and Antonovic et al. [37]. However, this may be due to the optimization of the model parameters [24]. For instance, using the model parameters in the Wenzl and Wilkens [28] OER model one assumes that all hypoxic in vitro data has a controlled pO 2 of 0.01 mmHg. To tune the existing OER models and judge which one is best, more experimental data is required. To reduce the uncertainty in experimental OER data, it is necessary to estimate LET values accurately and consistently between experiments, preferably using Monte Carlo simulations of the irradiation setup [38,39]. In addition, as the OER is dependent on the tissue type, there may be differences between the clinical OER and the in vitro OER, which may influence the results since the OER model parameters are generally based on in vitro data. However, Wenzl and Wilkens [28]   showed that there was generally good agreement between their OER model and results from preclinical and clinical studies (mainly performed applying photon irradiation). We therefore decided to use this model in our calculations, although uncertainties in the OER model must be acknowledged.
In the current OER model, the parameters were designed for a cell survival of 10%, ignoring a possible dose dependency. Wenzl and Wilkens [40] studied this dependency and found that the OER did not vary significantly with the dose, especially at high LET values and high pO 2 values. At lower LET values there was a clear dose dependency, but this was most pronounced for extremely hypoxic cells with lower pO 2 than what was seen in our patient case. The dose dependency also differed significantly between different cell lines. We therefore chose to maintain the survival level at 10% in our calculations, even though this may be a minor shortcoming. To extend the model to other survival levels, more experimental data is needed to fit the model parameters. Scifoni et al. [23] and Tinganelli et al. [25] did also ignore the dose dependency of the OER, although justified as their method was applied to high LET radiation.
The MC tool presented in this study depends on estimations of the pO 2 levels in a patient on a voxel-by-voxel basis. As our main objective was to demonstrate the tool performance, the pO 2 estimates should be regarded as an approximation in order to provide an indication of how pO 2 values will influence clinical scenarios. Toma-Dasu et al. [30] assumed that PET uptake as a function of pO 2 follows a sigmoid shaped curve, as the inhibition of chemical reactions typically follows this form. They showed their method for [ 18 [30]. This method has also been applied for other PET tracers, including [ 18 F]-HX4 PET [41]. While uncertainties are expected, the pO 2 conversion curve was only applied to the tumor area and its proximity to save simulation time, and the pO 2 of the remaining tissue was set to 60 mmHg. Thus, uncertainties in pO 2 estimation will not have significant impact on the dose to the normal tissue. Also, the lowest pO 2 value in the tumor applied in this study was estimated to be 2.2 mmHg, which corresponds to an OER of approximately 1.6. If the conversion curve (Appendix A, Fig. A.1) was less steep, as the data may suggest, this OER value would be an underestimation. An underestimation may not give high enough dose escalation to the tumor; however, it will be safer for the normal tissue.
The method to account for hypoxia in our study is, as mentioned, similar to the one on which the kill painting algorithm is based [23,25]. Kill painting was first proved experimentally, using cell exposure chambers which allow irradiating of cell cultures under defined oxygenation conditions [23,25]. Later, kill painting was applied in a treatment planning study with TRiP98 and an artificially generated [ 18 F]-FMISO-PET derived hypoxia map of a patient [42]. In addition to this method, hypoxia has earlier also been accounted for in treatment planning studies applying dose and/or LET painting, i.e. where the hypoxic volumes receive either increased dose and/or increased LET values, respectively [43][44][45]. Malinen and Søvik [43] found that hypoxia dose painting appeared to be especially attractive for protons, while both LET and dose painting has the potential to increase the tumor control probability for heavier ions like lithium and carbon. Chang et al. [46] also studied dose painting and showed that dose painting in conventional radiotherapy based on [ 18 F]-FMISO PET is technically feasible and has clear benefits in terms of increasing the tumor control probability. Kjellsson Lindblom et al. [41] showed, in addition, that dose painting in photon therapy based on [ 18 F]-HX4 PET appears feasible in non-small cell lung cancer. Our tool does not yet include optimization of treatment plans based on the oxygenation levels; however, future work includes implementing the method in a FLUKA Monte Carlo based treatment planning system.

Conclusion
The Monte Carlo based tool introduced in this study facilitates the determination of potential regions with underdosage in a tumor due to hypoxia. By further implementing the method in a treatment planning tool, it will be possible to adapt prescription doses based on the local tumor oxygen level, also accounting for RBE effects through coupling with existing RBE models. The resulting biologically weighted dose at a hypoxic condition compared to that in a normoxic condition was lower by a factor similar to the OER. For the head and neck cancer patient, this resulted in large areas of underdosage compared to the prescribed dose. While variable normoxic RBE models generally estimate higher biologically weighted dose than the D RBE1.1 , areas of underdosage were also observed when applying the D ROR model adapted for hypoxia in biologically weighted dose estimations. The results in this study confirm that neglecting the effect of hypoxia in proton therapy potentially could compromise the expected tumor control probability and should therefore be kept in mind in clinical practice.