RBE Model-Based Biological Dose Optimization for Proton Radiobiology Studies

Purpose: The purpose of the current study was (1) to develop a straightforward and rapid method to incorporate a dose-averaged linear energy transfer (LET d )–based biological effect model into a dose optimization algorithm for scanned protons; and (2) to apply a novel beam delivery strategy with increased LET d within the target, thereby enhancing the biological effect predicted using the selected relative biological effectiveness (RBE) model. Materials and Methods: We first generated pristine dose Bragg curves in water and their corresponding LET d distributions for 94 groups of proton beams, using experimentally validated Geant4 Monte Carlo simulations. Next, we developed 1-dimensional dose optimization algorithms by using the Python programming language. To calculate the RBE of protons for biological dose optimization, we invoked the McNamara RBE model and applied the radiobiological parameters of the lung cancer H460 cell line with 137 Cs reference photons. Results: High-accuracy optimization results were obtained. The relative difference between the delivered dose and the prescribed dose was approximately within 6 1.0% in the target. In addition, we obtained the RBE enhancement within the target by applying the LET-painting technique. For example, considering a simple case in which 2 opposed downslope dose fields were superimposed to form a uniform dose in the 5- to 10-cm target region, the center RBE was 1.23 6 0.01, which was greater than the center RBE of 1.16 6 0.01 found when using the traditional method of delivering 2 opposed flat dose fields. Conclusion: We have successfully developed an easy-to-implement method to perform the biological dose optimization procedure by invoking the McNamara RBE model in the iteration process using the Python programming language. According to the RBE model predictions, we conclude that the increased target LET d enhances the RBE. The accuracy of the RBE model predictions needs to be validated in radiobiological experiments.


Introduction
The relative biological effectiveness (RBE) of protons to reference photons currently used in the clinic is assumed to be a constant 1.1 [1][2][3], regardless of proton beam characteristics, depth in tissue, dose per fraction, dose rate, target cell types, and biological endpoints. This assumption was derived from previous experimental data mainly obtained at the middle of a spread-out Bragg peak (SOBP) region, using passive-scattering proton therapy (PSPT) [1,2]. With the development of improved beam delivery techniques, active-scanning proton therapy is becoming more popular and will be predominant in the near future [4]. Using a beam-scanning technique facilitated by orthogonal sweeping magnets allows intensitymodulated proton therapy (IMPT) to be feasibly implemented in the clinic. Moreover, the biological advantage of proton therapy may be enhanced in IMPT compared with traditional PSPT [5]. For example, if the detailed biological characteristics of each beamlet are integrated within the optimization process, biologically optimized IMPT plans can be achieved with an isometric therapeutic effect on the target tumors. However, such biologically optimized IMPT processes require appropriate biophysical models. Many biophysical models have been developed for proton therapy to predict the RBE [6][7][8][9][10][11][12][13][14][15]. Nevertheless, none of these models has ever been applied in any clinical proton therapy treatment planning system to date owing to the discrepancies among model predictions and between models and experimental data [16,17]. This notwithstanding, nearly all of these models and most of the recent experimental data do not support the assumption of a constant proton RBE, but instead support the conclusion that proton RBE varies spatially with respect to the variation of dose and linear energy transfer (LET) [18,19].
Before clinics can use true biologically optimized IMPT plans, more research is required to overcome the many difficulties that arise in the construction of accurate RBE models. First, most of the existing biophysical models were developed by using experimental data generated with PSPT [2], which may not accurately predict the biological consequences for scanned beams owing to the significant differences in beam characteristics between the two [18]. Second, accurate biological effect data for pristine proton beams are limited and difficult to obtain [18]. Moreover, most of the experimental validations of RBE models stay at the stage of irradiating biological samples across the flat physical dose 1 SOBP; few studies have been performed under a uniform biological dose 2 SOBP, the results of which provide much stronger evidence for proving the accuracy and prediction ability of a selected RBE model. Therefore, in the current study, we sought to develop a straightforward and fast method to generate a uniform biological dose SOBP that can be used to validate the accuracy of a selected RBE model. We invoked the McNamara [14] RBE calculation model in the biological dose optimization process to demonstrate the feasibility of our method. The involvement of other models in the biological dose optimization process will be reported in future work.
As in most other RBE models [6,8], the McNamara model predicts that the proton RBE increases with increases in doseaveraged LET (LET d ) [14]. This prediction motivated us to investigate the feasibility of enhancing the target RBE by increasing the LET d therein. With IMPT, the weights of the beamlets can be freely adjusted to form any reasonably desired target dose distribution. Because in general low-energy protons have a relatively high LET, increasing the weights of low-energy beams inside the target can be an effective way to increase the LET d therein. We used 1 simple case to prove this concept: the superposition of 2 opposed downslope dose fields to form a flat target dose distribution. We investigated how these opposed ramped fields affect the target RBE according to the McNamara model predictions.
The current study thus provides the fundamental design methods for performing proton biology experiments. The biological effects from the RBE model predictions reported herein will be validated by using cell irradiation experiments in future work.

Basic Settings in Monte Carlo Simulations
The general-purpose Monte Carlo toolkit Geant4 [20,21] (version 10.3.p03) was used to perform the particle tracking to generate the depth dose curves and the depth LET d curves for our proton beams. These raw data were used as the input to the physical and biological dose optimization process. Many validated reference physics lists that contain both of the electromagnetic and hadronic physics processes can be applied in proton and particle therapy studies. We compared the simulation results by using 3 typical physics lists, namely, ''QGSP_BERT,'' ''QGSP_BIC,'' and ''FTFP_BERT,'' and found the dose difference was below 1% for protons within the therapeutic energy ranges. In this study, we selected the FTFP_BERT physics list as a representative.
The 94 groups of scanned proton beams used at The University of Texas MD Anderson Proton Therapy Center were selected for the calculations. The energy in these groups varies from 72.5 to 221.8 MeV with a range (depth with 90% of the peak dose in the distal falloff) of 4.0 to 30.6 cm in water.
The spot-scanning nozzle and an 80 3 80 3 40 cm 3 water phantom were modeled by using Geant4. We defined 2 cylindrical scorer sizes in the water phantom for different purposes. A scorer with a small radius of 4.08 cm and thickness of 0.01 cm was used to score the dose that can be validated with the measurements using a Bragg peak chamber (PTW 34070, SN 0016; PTW-Freiburg, Freiburg, Germany) with a radius of 4.08 cm. A scorer with a larger radius of 40 cm and thickness of 0.01 cm was built to cover the large size of the beam spot and to maximally capture the laterally scattered particles to subsequently score the so-called integral depth dose (IDD). The production cut of secondary particles (photons, electrons, positrons, and protons, including all recoil ions) was set to 0.1 mm to match the smallest length (thickness in this case) of the virtual detector recommended by the Geant4 collaboration [20,21]. In addition, the LET d was scored by using the larger scorer. The LET d calculation method for each single beamlet using Geant4 was described in our previous study [22]. Only primary protons were considered in calculating LET d to make it consistent with the definition of LET by the International Commission on Radiation Units and Measurements in which LET only considers the electronic interactions of primary charged particles [23]. The secondary protons are not included in the LET calculation in our studies because they are generated mostly from nonelastic nuclear reactions [24], rather than from electronic interactions.
The number of primary source particles was set as 10 7 for each beamlet to make the calculated dose meet the statistical uncertainty requirement (relative error of the mean value ,1%) when the dose is larger than 5% of the peak dose. All associated simulation data were then written to ROOT histograms [25].

The McNamara Model for Predicting the RBE of Protons
The phenomenologic McNamara model [14], which was recently developed on the basis of all in vitro proton datasets published through 2015, was used to calculate the RBE for protons in the current study. The RBE calculation formula in the McNamara model [14] is straightforward and requires only the a and b values from reference photons (a x and b x ) for the specified cell line, as well as the dose and LET d from selected proton beams. The McNamara RBE calculation formula is expressed in Equation (1).
Here the fit coefficients are p 0 ¼ 0.99064 (standard error: 0.014125), p 1 ¼ 0.35605 (standard error: 0.015038), p 2 ¼ 1.1012 (standard error: 0.0059972), and p 3 ¼ À0.0038703 (standard error: 0.00091303). The proton dose D p and LET d can be obtained from Monte Carlo simulations for each specified spatial location. We apply the photon coefficients obtained from our previous study to the RBE calculations. For the non-small cell lung cancer cell line H460 with 137 Cs photon irradiation, a x ¼ 0.29 6 0.016 Gy À1 , b x ¼ 0.083 6 0.00395 Gy À2 , and (a/b) x ¼ 3.49 6 0.255 Gy [18].
Based on the error propagation rule, the error of the calculated RBE using Equation (1) depends on the partial derivatives and errors of the following 7 parameters: D p , LET d , D p , p0, p1, p2, and p3. We assumed the statistical errors of D p and LET d are negligible, and then incorporated partial derivatives and errors of other 5 parameters into the Python code to calculate the error of RBE.
According to the definition of RBE for cell survival (RBE ¼ D x /D p ), the RBE-weighted dose from particles, RBE 3 D p , is equal to the dose from the reference photons, D x , at the same cell-surviving fraction. Therefore, we can calculate the cellsurviving fraction S for particles at dose D p by using the linear-quadratic model for photons at dose D x , expressed in Equation (2), without knowing the a and b values for the particles [3,26].
The calculated cell-surviving fraction can be compared with the experimental data from clonogenic assays to further validate the accuracy of the applied RBE calculation model.

Physical Dose and RBE-Weighted Dose Optimization Algorithms
We used the Python programming language (version 3.4.3) and its NumPy and SciPy libraries to perform the physical dose and biological dose optimization procedures. The basic principle of an optimization algorithm is to find the optimal solutions for the weight of each beamlet to form the desired dose distribution within the target, using an iterative scheme.
The 1-dimensional physical dose optimization is relatively simple. Initially, the IDD data generated from Monte Carlo simulations are read into the Python program and each IDD is then normalized by its peak dose. Therefore, without loss of generality we assume the peak dose is 1.0 Gy after normalization if we assign the units of Gy to dose because such assignation is arbitrary. This normalized IDD dataset is then used to solve for the beam weights. We specify the boundaries of the target region and the corresponding target dose distribution, keeping only the beamlets with penetration ranges that can cover the target region when finding the beam weights. For outlying energy groups, the weight is directly set to a very small value, that is, 1.0E-10.
The purpose of the optimization process is to find the optimal beam weights to minimize our objective function. In our optimization algorithm, the Python function curve_fit() from SciPy was selected as the optimizer to solve the beam weights, using the nonlinear least-squares minimization scheme.
The objective function for physical dose optimization is expressed in Equation (3).
Here i is the index of the beamlets in the available energy groups; w i is the weight of the ith beamlet, and it must be a nonnegative value; and D i (z) is the normalized dose at depth z in the ith beamlet's IDD curve. The spatial location z is limited within the target region from z min to z max . We set a relative error of 0.1% for the prescribed physical dose D t (z) within the target to increase the importance of these data points during the data fitting. We can use a simple example to interpret what energies are needed for dose optimization. For a target region with z from 5 to 10 cm, we need only the beams from 81.4 MeV (energy ID ¼ 11 and range ¼ 5.0 cm) to 118.6 MeV (energy ID ¼ 46 and range ¼ 10.1 cm), which means the index i ranges from 11 to 46.
After optimization, the delivered physical dose D r (z) at a specified location z can be calculated by using Equation (4).
The relative difference between the delivered physical dose and the prescribed physical dose for each specified location z within the target region can be measured by using Equation (5).
The RBE-weighted dose optimization process is more complex because the RBE is also a function of the physical dose and LET d , as indicated in Equation (1). Our aim is to find beam weights that make the product of the physical dose and RBE approach the prescribed biological dose value. The objective function of biological dose optimization can be expressed as shown in Equation (6).
Here Bio_D t (z) stands for the prescribed biological dose at location z, with units of Gy (RBE). Similarly, the relative error of 0.1% is also set for Bio_D t (z) within the target. According to Equation (1), the RBE value at location z, calculated by using the McNamara model, can be expressed as RBE . Here the delivered LET d after optimization can be calculated by using Equation (7).
Here LET d,r (z) is the delivered LET d at depth z; LET d,i (z) is the LET d at depth z from the ith beamlet in the raw pristine LET d dataset; w i is the weight of the ith beamlet; D i (z) is the normalized dose at depth z in the ith beamlet's IDD curve; w i 3 D i (z) is the weighted dose from the ith beamlet; and R i w * i D i ðzÞ; LET d ;r ðzÞ is the delivered physical dose at depth z after optimization. Therefore, the biological dose optimization involves iterative calculations using Equations (4), (7), and (1) to minimize the objective function expressed in Equation (6). The relative difference can also be calculated by using Equation (5), but the physical dose should be replaced with the corresponding biological dose. Relative difference within 62.0% is acceptable in the current study.

Depth Dose and LET d Distributions Generated from Validated Monte Carlo Simulations
Because the Geant4-based Monte Carlo system is used to generate physical data that will be used to perform radiobiological experimental designs, it is necessary to validate the calculated dose distributions with physical measurements. The spot-scanning nozzle at The University of Texas MD Anderson Cancer Center Proton Therapy Center in Houston has been commissioned and is thus appropriate for validating dose distributions [27]. Here we compare only the depth dose distributions between Monte Carlo simulations and measurements. To avoid clutter in the data presentation, Figure 1A shows only the dose data for 14 selected beam energies, although the data from all 94 groups of beam energies are available. The data with marks are from the measurements using a large-radius (r ¼ 4.08 cm) plane-parallel ion chamber Bragg Peak Chamber (PTW 34070, SN 0016) in a PTW MP3 water tank (PTW-Freiburg). The data with lines are from Monte Carlo simulations using a scorer with the same radius as the detector. The dose is expressed as Gy per monitor unit (MU). Ten data points at different locations along the raw measured and calculated Bragg curves (r ¼ 4.08 cm) of each beam energy were sampled to calculate the conversion factor between number of primary source particles (used in Monte Carlo simulations) and MU (used in beam delivery). The proton source parameters were tuned to make the calculated dose curves match the measured curves. After multiple adjustments, the standard error of the conversion factor (number of protons/MU) was reduced to between 0.1% and 0.4% for all of the 94 beam energies, which indicates the good agreement between Monte Carlo calculations and measurements.
Because the spot size of some low-energy beamlets is larger than the detector size (radius ¼ 4.08 cm), the measured IDD must be corrected to consider the dose contributed by particles outside of the detector radius, including both the primary protons and laterally scattered particles. Figure 1B shows the corrected IDDs of 14 selected beam energies from Monte Carlo simulations using a virtual detector with a radius of 40 cm to maximally capture all particles radially. The product of the dose and the area of the front surface of the virtual detector (radius ¼ 40 cm) is expressed as GyÁcm 2 /MU. Figure 1C shows the LET d data of 14 selected beam energies. The IDD and LET d datasets for all 94 beam energies will be used as the basic input data for biological dose optimizations.

A Uniform RBE-Weighted Dose Distribution in the Target
We set the prescribed RBE-weighted dose level to be 3.80 Gy (RBE) within the SOBP (5-10 cm), corresponding to the reference photon dose at which the cell-surviving fraction is 10% for the H460 cell line, as presented in our previous study [18].
Thirty-six groups of beam energies are needed for biological dose optimization ranging from 81.4 MeV (energy ID ¼ 11 and range ¼ 5.0 cm) to 118.6 MeV (energy ID ¼ 46 and range ¼ 10.1 cm). The optimization results, including delivered physical dose, delivered biological dose, RBE predicted by using the McNamara model within the 95% confidence interval (CI) (61.96 r), and surviving fractions of the H460 cell line predicted by using the linear-quadratic model, are shown in Figure 2A and 2B. The predicted surviving fraction within the SOBP was 0.1, as expected. The relative difference between the prescribed and delivered biological dose within the target is shown in Figure 2C. The relative difference within the target ranged from À1.2% to 0.8%, but at most of the locations, the relative difference was within 60.2%. In Figure 2D, the standard error of RBE (%) is below 1.5% proximal to the dose falloff, but thereafter it increases sharply to the maximum value of 5.0%.

Enhancement of Biological Effects Using Opposed Ramp Dose Fields
For multi-field IMPT, it is not necessary to apply a flat dose SOBP for each field as long as the accumulated dose from multiple fields is uniform in the target. However, different beam delivery strategies with different beam weights can alter the LET d distribution within the target. According to Equation (1), when all other parameters are fixed, RBE is a monotonically increasing function of LET d . Therefore, if we can apply a new beam delivery strategy to increase the LET d within the target while maintaining the target dose level, the biological effect within the target should be increased, and thus, this model-predicted RBE enhancement must be experimentally verified.
Bassler et al [28,29] and Fager et al [30] have investigated the therapeutic effect of applying the LET-painting technique, which aims at placing more low-energy, and thus high-LET, particles in the target tumor region to enhance the biological effect therein [28,29] or to reduce radiation dose with equivalent clinical effectiveness [30]. We first report the results from the traditional method of superimposing 2 opposed flat dose fields to form a flat target dose SOBP, called ''case A'' henceforth. Then, we report results from a simple case of LET painting made by superimposing 2 opposed downslope dose fields to form a flat target dose SOBP, called ''case B'' henceforth, which is similar to the gradient dose distribution reported by Krämer and Jäkel [31] for carbon ions.
The single-field delivered dose and LET d distributions in the target of 5 to 10 cm for both case A and case B before the superimposition of 2 opposed fields are compared in Figure 3A and 3B. LET d within the target from a downslope dose field was enhanced as compared with the traditional flat dose field. The delivered dose and LET d distribution after the superimposition of 2 opposed fields are compared in Figure 4A and 4B. In case A, although LET d at the boundaries increased sharply, LET d was much lower at most of the central target region, with the value of 2.63 keV/lm at the middle point of the target. In case B, the LET d distribution was notably different and varied slowly across the target region, with the value of 4.34 keV/lm at the middle point of the target. The RBE variation in Figure 5A and 5B follows the same trend as the LET d variation in Figure 4A and 4B. The RBE at the target center in case A was 1.16 6 0.01, whereas in case B it was 1.23 6 0.01, and the RBE enhancement ratio of case B to case A was 1.23/1.16 ¼ 1.06 at the target center. However, in case B, the out-of-target dose increased and the dose gradient at the boundary of the target was not steeply contoured, which may limit the clinical application of this dose delivery strategy. The 95% CIs (61.96 r) of RBE for case A and B are shown in Figure 5C and 5D. No obvious overlap can be found between these 2 CIs within the central target region. Hence, we assume these 2 setups can result in distinguishable RBE values. We are planning to experimentally verify this enhanced RBE effects in case B compared with case A. The experimental results will be reported in our future work.
As shown in Figure 3A, the dose gradient was (1.9À0.1) Gy/5 cm ¼ 0.36 Gy/cm. We increased this dose gradient by narrowing the width of the target region. The physical and biological dose results for the target from 5 to 7 cm are compared in Figure 6A and 6B. The dose gradient was then (1.9À0.1) Gy/2 cm ¼ 0.9 Gy/cm for case B. Compared with the results shown in Figure 5, the RBE was increased for both cases. For case A, the 5-to 7-cm target region had a center RBE of 1.21 6 0.01, which was greater than the RBE of 1.16 6 0.01 in the 5-to 10-cm target region. For case B, the narrower target had a center RBE of 1.30 6 0.02, which was greater than the RBE of 1.23 6 0.01 in the 5-to 10-cm target region. The 95% CIs (61.96 r) of RBE for case A and B are shown in Figure 6C and 6D. Similar to Figure 5, no obvious overlap can be found between these 2 CIs within the central target region. The relative difference within the target ranges from À1.2% to 0.8%, but at most locations, the relative difference is within 60.2%. (D) The standard error of RBE (%) is below 1.5% proximal to the dose falloff, but thereafter it increases sharply to the maximum value of 5.0%.

Discussion
For a given RBE calculation model and associated biological dose optimization results, we can use the biologically optimized beam delivery strategy to verify the accuracy of the predicted biological effects. If the model prediction is sufficiently accurate, the selected cell samples at different locations across the biologically optimized SOBP are expected to have identical cellsurviving fractions within the specified CI, that is, 95%.
Only the non-small cell lung cancer H460 cells were selected in this report to demonstrate the feasibility of our method. More clinically relevant cell lines will be tested in our future studies to validate the applicability of the selected RBE model. We must note that each cell line has its specific radiobiological parameters a and b due to the differences of their radiosensitivity [32]. Therefore, each cell line has its unique biological dose optimization results even if the same surviving fraction is specified.  The optimization time depends on the width of the SOBP, which is directly related to the number of beam energies. In addition, the configuration of the computer can also affect the optimization time. The optimization tasks in this study can be completed rapidly even by using a commonly configured personal computer. The optimization time for the physical dose SOBP in Figure 3A was 7.3 seconds, while for the biological dose SOBP in Figure 2A, the time was increased to 22.4 seconds. The optimization time varies slightly in different runs even for the same optimization task because the time also depends on the current resource allocations of the computer.
The RBE enhancement predicted by using the McNamara model for the target region irradiated by 2 opposed downslope dose fields is caused by the increased weights of the low-energy beams, which have relatively higher LET values. However, the target RBE is enhanced at the cost of the unneeded dose outside of the target as shown in Figures 5 and 6, which can be a limiting factor to apply this beam delivery strategy to clinical treatment plans. Nevertheless, the construction of treatment plans is out of the scope of this article. Although we only focus on the biological response of cancer cells within the flat dose region in this study, we will also investigate the response of normal tissues outside of the flat dose region in future studies. In addition, the increased RBE within a smaller target is not directly caused by the dose or dose gradient characteristics, but rather is caused by the increased LET d within the target, which can be attributed to the increased ratio of low-energy and high-LET protons in the target. Although we have obtained the high-accuracy optimization results, we must point out that one of the limitations of the method developed in this study is that it can only be applied in the uniform medium such as water. We expect to extend the application to complex patient geometries by considering the effect of heterogeneities in dose and LET d calculations.

Conclusion
We have demonstrated the feasibility and flexibility of applying the McNamara proton RBE calculation model into a biological dose optimization algorithm using the Python programming language. High-accuracy optimization results have been obtained. The methodology developed in the current study can be extended to apply to other RBE calculation models for different modalities of charged particle therapy to generate biologically optimized treatment plans. Such work will be reported in a future study.
In addition, we have demonstrated RBE enhancement from using the LET-painting technique. Using 2 opposed downslope dose fields has the potential to increase the LET d within the target so as to enhance the biological effect therein, and this RBE enhancement can be quantified by using appropriate RBE calculation models. Validating the accuracy of these RBE model predictions in cell irradiation experiments will be our next step.