Rapid implementation of the repair-misrepair-fixation (RMF) model facilitating online adaption of radiosensitivity parameters in ion therapy

Introduction: Treatment planning for ion therapy must account for physical properties of the beam as well as differences in the relative biological effectiveness (RBE) of ions compared to photons. In this work, we present a fast RBE calculation approach, based on the decoupling of physical properties and the αx/βx ratio commonly used to describe the radiosensitivity of irradiated cells or organs. Material and methods: In the framework of the mechanistic repair-misrepair-fixation (RMF) model, the biological modeling can be decoupled from the physical dose. This was implemented into a research treatment planning system for carbon ion therapy. Results: The presented implementation of the RMF model is very fast, allowing online changes of αx/βx. For example, a change of αx/βx including a complete biological modeling and a recalculation of RBE for 2.9⋅105 voxel takes 4 ms on a 4 CPU, 3.2 GHz workstation. Discussion and conclusion: The derived decoupling within the RMF model allows fast changes in αx/βx, facilitating online adaption by the user. This provides new options for radiation oncologists, facilitating online variations of the radiobiological input parameters during the treatment plan evaluation process as well as uncertainty and sensitivity analyses.


Introduction
In carbon ion therapy, biological models are needed to predict the relative biological effectiveness (RBE) of carbon ions and lighter particles generated by nuclear fragmentation events. In the framework of the linear-quadratic (LQ) model (Kellerer and Rossi 1978) the cell survival is described by two radiosensitivity parameters: α (linear) and β (quadratic). Biological models have been developed to estimate the particle radiosensitivity parameters (α p and β p ) from the tissue specific, x-ray reference values (α x and β x ). In the LQ model, RBE can be described by a function of α x , β x , α p , β p and the physical dose d. Several biological models have been proposed and are currently used for the prediction of carbon ion RBE: different versions of the local effect model LEM1-LEM4 (Scholz et al 1997, Grün et al 2012, the microdosimetric kinetic model (MKM) (Hawkins 1994, Kase et al 2008 and the repair-misrepair-fixation (RMF) model (Carlson et al 2008, Frese et al 2012. The LEM1 and the MKM are currently used in clinical practice. Typically, biological models, as well as their main input parameters, α x and β x , suffer from large uncertainties (Weyrather et al 1999, Carlson et al 2008, resulting in large confidence intervals for α p and β p .
At the moment these uncertainties have to be accounted for by the experience of radiation oncologists while three-dimensional evaluation of these uncertainties in treatment planning systems is very challenging. The current practice is a treatment plan evaluation based on model predictions for one specific set of α x and β x . The radiation oncologists have to include uncertainties in α x and β x as well as in the biological modeling itself into the dose prescription and treatment plan evaluation without a possibility to quantify or simulate different biological modeling results. One of the main reasons is that the biological modeling process is computational intensive and takes long. The biological modeling process is generally slow as it has to be simulated over a broad range of ion types and energies, present in ion beams. This drawback is commonly overcome by precalculated and tabulated biological model output. Nevertheless, if e.g. α x or β x should be changed after an optimization, the complete information needs to be generated again using another table. This does not allow frequent recalculations with different sets of radiosensitivity parameters, which would be needed for a biological uncertainty estimation of a carbon ion treatment plan.
In order to overcome this limitation in ion treatment plan evaluation, we present a decoupling approach within the framework of the repair-misrepair-fixation (RMF) model (Carlson et al 2008). The biological modeling is decoupled from the physical dose calculation and can be implemented very efficiently, facilitating online changes of x-ray radiosensitivity parameters in real time.

Treatment planning for carbon ion therapy
We implemented a multifield biological optimization for carbon ion therapy into an extension of the research treatment planning platform CERR (computational environment for radiation research) (Deasy et al 2003, Schell and Wilkens 2010, Brüningk et al 2015. CERR is a Matlab-based research treatment planning software. The presented workflow was implemented and evaluated using Matlab version 2013b. The ion therapy extension includes biological effect optimization, using a spot scanning beam delivery and a pencil beam dose algorithm. The calculations are based on tabulated Monte Carlo generated fragmentation spectra (Parodi et al 2012). These spectra contain data for all relevant ions types with atomic numbers Z from 1 (protons) to 6 (carbon ions). For the calculations we implemented influence matrices which describe the influence of every spot j on every voxel i for dose d ij , particle fluence and biological parameters (such that e.g. the dose in voxel i becomes d i = j d ij w j for spot-weights w j ). The common, standard way to facilitate RBE calculations in a treatment planning system uses precalculated and tabulated biological modeling results (α p and β p ) as function of water-equivalent depth for every available initial carbon ion energy. This means tables of α p and β p values for every initial energy and every considered tissue type (characterized by α x and β x ) are needed and stored with the tabulated depth dose curves. Depending on the chosen α x and β x the right tables are used for the calculations, resulting in potentially inflexible situations (only precalculated α x -β x combinations can be used and once a combination is set it cannot be easily changed without redoing most of the calculation steps again).

The RMF model and the introduced decoupling approach
The biological modeling and the decoupling were achieved with the RMF model. The RMF model was introduced by Carlson et al (2008) and successfully implemented and validated for carbon ion, helium ion and proton treatment planning (Frese et al 2012, Mairani et al 2016. In the presented implementation, the mechanistic RMF model uses Monte Carlo damage simulation (MCDS) estimates of absolute double-strand break (DSB) yields (Semenenko and Stewart 2004, Hsiao and Stewart 2008, Stewart et al 2011. The link between DSB induction and cell death has been discussed in-depth previously by Carlson et al (2008), Stewart et al (2011Stewart et al ( , 2015 and references therein. Carlson et al (2008) and Frese et al (2012) showed how to determine particle radiosensitivity parameters α p and β p from the reference radiosensitivity parameters α x and β x and the MCDS output DSB induction Σ and frequency-mean specific energy z F . The corresponding values for the reference radiation (Σ x and z Fx ) are simulated based on the secondary electron spectra of a Co-60 source (Hsiao and Stewart 2008). MCDS 3.10A and the default cell nucleus diameter d tar = 5 μm are used. Kamp et al (2015) combined the work of Frese et al (2012) with fragmentation spectra of carbon ion beams and implemented the RMF model in the standard, precalculated and tabulated way in a research treatment planning system, in order to show and discuss the influence of fragmentation spectra on RBE predictions by the RMF model.
The presented decoupling is based on the following equations (see also equations (3) and (4) of the work by Kamp et al (2015)). In order to differentiate between the different spots j and voxel i, the dose averaged radiosensitivity values are written as α p,ij and β p,ij . The depth z is included in the voxel i. Φ is the fluence spectrum of all fragmentation particles and SP the stopping power, both are a function of the energy E and the type (atomic number Z) of the ion.
The RMF model with its decoupling parameters c 1 and c 2 (as introduced in equations (2) and (4) has the very convenient property that α x,i and β x,i can always be factored out. This means that in every voxel the biological radiosensitivity parameters can be separated from the physical beam properties (d, Φ, SP,z F , Z ) of the ions and their simulated DSB induction Σ(Z, E). This decoupling can be done for any desired primary ion type, including their different fragment spectra and would equally work for Monte Carlo based dose calculations. c 1 can also be interpreted as the dose-weighted RBE DSB , defined as ratio of the DSB inductions: In equations (2) and (4) a very important relation becomes apparent. The needed values for α p and β p are only dependent on c 1 and c 2 (besides the dependence on α x and β x which are used to describe different tissue types). In the study by Kamp et al (2015) equations (1) and (3) were introduced, not mentioning and not exploiting the key property that α p and β p depend on the same two integrals (abbreviated with c 1 and c 2 ). A direct consequence is that for a given α x and β x , α p and β p can be directly converted to c 1 and c 2 and vice versa. In the recent study, the c 1 -c 2 -formalism is introduced to fully exploit the simplicity which results from this property for RBE and RBE-weighted dose (RWD) calculations in ion therapy. This relation is, too our knowledge, unique to the RMF model and a key feature for several shown and discussed implementations in the scope of this work.
The following equations give an overview of the implemented matrices as well as the physical and biological quantities that can be calculated in every voxel (i) for a given set of spot-weights ω j . The calculation of linear energy transfer (LET), α p,i and β p,i is based on the work of Oelfke (2004, 2006).
The dose-weighted calculation of c 1,i and c 2,i in the ij formalism is very convenient, because it can be implemented in a treatment planning system exactly the same formalism as d i , LET i , α p,i and β p,i . The same geometry, raytracing and water equivalent depths determination can be used for c 1 and c 2 as for e.g. the dose. This also means that the depth dependence along a pencil beam of c 1 and c 2 can also be precalculated and tabulated, identical to d, LET, α p and β p for different initial carbon ion energies. From the memory and time efficiency this means that using the introduced formalism, α p and β p can be substituted by c 1 and c 2 . Any desired α x and β x can be multiplied later on. For other commonly used models, this is not possible because α x and √ β x cannot be always factored out. Equations (9)-(12) are derived in the Appendix. The here presented formalism is based on the property of the RMF model that α x,i and β x,i can always be factored out and, hence, even be changed after an optimization without the need for much additional computation time. Using the RMF model, α p,i and β p,i can be expressed similar to equations (2) and (4).
Solving equations (11) and (12) for c 1,i and c 2,i shows how to determine these values from known α p,i and β p,i distributions. A change of α x and β x can be done very fast. It is based on α p,i and β p,i only, without even the need to know the ij matrices.
These equations have a favorable consequence: it is actually not necessary to implement c 1,ij and c 2,ij explicitly. The standard α p,i and β p,i implementation can be used with an arbitrary value for α x,i and β x,i and converted at any time to c 1,i and c 2,i . This is useful to stay compatible with existing implementations of other biological models in a ion treatment planning system. Note that this calculation of c 1,i and c 2,i is independent of d i , although d i is present in the doseweighting in equations (9) and (10). Similar to equations (13) and (14), c 1,ij and c 2,ij can be calculated from α p,ij , β p,ij , α x,i and β x,i (α x,i and β x,i are independent of the spot j).
The list of quantities that can be calculated in every voxel is completed by RBE and RWD. The dose d i is the dose per fraction. The general calculation of RBE LQ,i with α p,i and β p,i being a function of α x,i , β x,i and parameters of the used biological model can be rewritten using the decoupling of the RMF model by inserting equations (11) and (12) into equation (15). Hence RBE RMF,i can be calculated as In the framework of the RMF model, the calculation of RBE RMF,i is only dependent on α x,i /β x,i and not on α x,i and β x,i individually. This was previously shown by Frese et al (2012) without the decoupling approach. The RWD i is calculated voxel-wise as

Results
We implemented the introduced c 1 and c 2 decoupling formalism and tested its performance. Figure 1 shows an example of a carbon ion treatment plan for an astrocytoma patient, previously treated with photons. The PTV was optimized on RWD = 3 Gy(RBE) using two fields as in Kamp et al (2015). Panels A to D in figure 1 are the commonly shown optimization results RWD i , RBE i , d i and LET i , respectively. The introduced decoupling parameters c 1,i and c 2,i are shown in panels E and F, together with α p,i and β p,i in panels G and H, respectively. A spatially constant α x,i /β x,i = 2 Gy (α x,i = 0.1 Gy −1 and β x,i = 0.05 Gy −2 for all i) is used in this example treatment plan. These values are commonly used for sarcoma of the scull base (Schulz-Ertner et al 2007).
An evaluation possibility for online changes in the biological parameters is shown in figure 2. RBE-weighted dose-volume histograms (RWDVHs) are displayed for two representative structures. The initial RWDVH, optimized using α x,i /β x,i = 2 Gy can be compared to the result of two modified α x,i /β x,i . Note that the presented +20% change in α x /β x can be due to a 20% increase in α x or a 20% decrease in β x . RBE and hence RWD are only dependent on α x /β x and not α x and β x independently (equation (16)).
Considering the whole patient, the values of the decoupling parameter c 1 range from 1.18 to 2.38 (mean 1.40, standard deviation 0.27, median 1.27), the values for c 2 from 0.37 Gy to 25 Gy (mean 2.64 Gy, standard deviation 3.71 Gy, median 0.75 Gy). The histograms and scatter plots in figure 3 show the distribution of resulting c 1 and c 2 values and their dependency on LET. The left part shows the values for all voxel inside the PTV structure (in total Axial CT slice of a treatment plan using the RMF model. The astrocytoma plan with two carbon ion fields was optimized on 3 Gy(RBE) with the ion therapy extension of the CERR treatment planning system using a spatially constant α x,i /β x,i = 2 Gy (α x,i = 0.1 Gy −1 and β x,i = 0.05 Gy −2 ). The PTV is shown in red, along with 3 organs at risk: left optic nerve (green), left eye (orange) and left lens (brown). The panels show (A) RWD i , (B) RBE i , (C) physical dose d i and the dose-weighted LET i in D. The two decoupling variables c 1,i and c 2,i are shown in panels E and F, along with the doseweighted α p,i and β p,i in panels G and H.
3.5 · 10 4 voxel), the right part all voxel in the body structure without the PTV which receive a RWD i > 0.1 Gy(RBE) (8.4 · 10 4 voxel). Both c 1 and c 2 are closely correlated to the LET distribution (Pearson correlation coefficients ρ > 0.97). The slopes of the linear regression lines for c 1 are 0.0087 μm keV −1 and 0.0078 μm keV −1 for PTV and 'body without PTV', respectively (R 2 > 0.96). The corresponding slopes for c 2 are 0.148 Gy μm keV −1 and 0.156 Gy μm keV −1 for PTV and 'body without PTV', respectively (R 2 > 0.98). Hence, the higher the LET value is the greater are c 1 and c 2 . This trend can also be seen in the increase of α p and β p distal to the PTV, which has been discussed comprehensively by Carlson et al (2008) and (Frese et al 2012). Figure 3 shows that the range of LET values and hence the range of c 1 and c 2 values inside the PTV is smaller than for the rest of the body. Considering the values in the PTV, LET ranges from 43 to 124 keV μm −1 (mean 63 keV μm −1 , standard deviation 14.0 keV μm −1 , median 59 keV μm −1 ), c 1 from 1.60 to 2.23 (mean 1.79, standard deviation 0.122, median 1.76) and c 2 from 3.86 to 16.7 Gy (mean 6.71 Gy, standard deviation 2.08 Gy, median 6.07 Gy). The corresponding values for the 'body without PTV' are LET from 16 to 160 keV μm −1 (mean 63 keV μm −1 , standard deviation 31.8 keV μm −1 , median 58 keV μm −1 ), c 1 from 1.28 to 2.38 (mean 1.74, standard deviation 0.254, median 1.72) and c 2 from 0.55 to 23.85 Gy (mean 7.24 Gy, standard deviation 0.500 Gy, median 6.36 Gy), respectively. The observed spread in c 1 and c 2 for constant LET (right panels of figure 3) originates from the fact that two different particle spectra can have the same dose-weighted LET value, but slightly different dose-weighted c 1 and c 2 and hence α p and β p .
After the optimization an online change of α x,i /β x,i can be done in 4 ms (for a stable timing in Matlab the mean value of 1000 changes was taken) if the value is changed throughout every voxel with RWD i > 0 Gy(RBE) (2.9 · 10 5 voxel). A change in only a subset of the structures is even faster, for example it took 1 ms (mean value of 1000 changes) for the PTV (3.5 · 10 4 voxel, 270 cm 3 ). The calculations were performed on a 4 CPU, 3.2 GHz workstation. In the presented decoupling approach an (online) change of α x /β x represents a change in the biological dose response behavior, including a full biological modeling, whereas physical values as for example local particle spectra, LET and physical dose are kept constant.

Discussion
The presented implementation of the RMF model allows very fast changes of α x,i and β x,i and hence online adaption of α x,i /β x,i and RBE RMF,i . The changes can be done for every voxel independently. The shown decoupling is straight forward to implement in treatment planning systems, as α x,i , β x,i , α p,i and β p,i are commonly used. The decoupling, represented by c 1 and c 2 can be calculated before or after the generation of the ij matrices or even after the optimization or after a forward calculation with given ω j . This can also be done for Monte Carlo based treatment planning, because the integration over the particle spectra (compare to equations (1) and (3)) are still independent of α x,i and β x,i . The calculation time is very fast as it essentially consists of the time Matlab needed for two element-wise multiplications and a summation of vectors for α p,i = α x,i · c 1,i + β x,i · c 2,i . For the calculation of β p,i = β x,i · (c 1,i ) 2 , two element-wise multiplications are needed. The length of the vectors is the number of considered voxel. The time needed to do the same changes with the standard implementation (tabulated α p and β p ) is not easily assessable, because the needed generation of new influence matrices and their subsequent handling is influenced by many other parameters. Common values for this lie in the range of several minutes. The presented fast implementation of the RMF model is applicable for all currently discussed ion beams for radiotherapy. Accounting for the fragmentation spectra of the ion beams is straight forward and does not affect the calculation time per voxel for online α x,i /β x,i changes.
There are several possible applications for the presented biological modeling implementation. First of all this approach facilitates online adaption of the reference radiosensitivity parameters including a full biological modeling for RBE calculation. This can be used for online treatment plan evaluation, which can now, for example, include worst case scenarios for the biological response of tumor and normal tissue. The decoupling can be used for a systematical assessment of uncertainties in the RBE or RWD distributions, originating from confidence intervals of α x and β x . This facilitates comprehensive sensitivity analyses in threedimensional, multifield geometries based on patient CT data (Böhlen et al 2012, Kamp et al 2014. In this context it has to be evaluated if it is sufficient to just change α x and β x , or if it is necessary to include further factors, like e.g. uncertainties in the biological modeling process. Fast biological modeling might in addition be useful for robust treatment plan optimization (Pflugfelder et al 2008, Unkelbach et al 2009, Bangert et al 2013 which aims to minimize the impact of uncertainties on a treatment plan.

Conclusion
The presented decoupling approach can be used for fast changes in the reference radiosensitivity parameters and is hence an important technical step towards online adaptation and evaluation of RBE-based treatment plans in ion therapy. This will help to estimate the impact of uncertainties in the biological modeling and the biological input parameters on treatment plans.