An empirical model of proton RBE based on the linear correlation between x‐ray and proton radiosensitivity

Abstract Background Proton relative biological effectiveness (RBE) is known to depend on physical factors of the proton beam, such as its linear energy transfer (LET), as well as on cell‐line specific biological factors, such as their ability to repair DNA damage. However, in a clinical setting, proton RBE is still considered to have a fixed value of 1.1 despite the existence of several empirical models that can predict proton RBE based on how a cell's survival curve (linear‐quadratic model [LQM]) parameters α and β vary with the LET of the proton beam. Part of the hesitation to incorporate variable RBE models in the clinic is due to the great noise in the biological datasets on which these models are trained, often making it unclear which model, if any, provides sufficiently accurate RBE predictions to warrant a departure from RBE = 1.1. Purpose Here, we introduce a novel model of proton RBE based on how a cell's intrinsic radiosensitivity varies with LET, rather than its LQM parameters. Methods and materials We performed clonogenic cell survival assays for eight cell lines exposed to 6 MV x‐rays and 1.2, 2.6, or 9.9 keV/µm protons, and combined our measurements with published survival data (n = 397 total cell line/LET combinations). We characterized how radiosensitivity metrics of the form DSF%, (the dose required to achieve survival fraction [SF], e.g., D10%) varied with proton LET, and calculated the Bayesian information criteria associated with different LET‐dependent functions to determine which functions best described the underlying trends. This allowed us to construct a six‐parameter model that predicts cells’ proton survival curves based on the LET dependence of their radiosensitivity, rather than the LET dependence of the LQM parameters themselves. We compared the accuracy of our model to previously established empirical proton RBE models, and implemented our model within a clinical treatment plan evaluation workflow to demonstrate its feasibility in a clinical setting. Results Our analyses of the trends in the data show that DSF% is linearly correlated between x‐rays and protons, regardless of the choice of the survival level (e.g., D10%, D37%, or D50% are similarly correlated), and that the slope and intercept of these correlations vary with proton LET. The model we constructed based on these trends predicts proton RBE within 15%–30% at the 68.3% confidence level and offers a more accurate general description of the experimental data than previously published empirical models. In the context of a clinical treatment plan, our model generally predicted higher RBE‐weighted doses than the other empirical models, with RBE‐weighted doses in the distal portion of the field being up to 50.7% higher than the planned RBE‐weighted doses (RBE = 1.1) to the tumor. Conclusions We established a new empirical proton RBE model that is more accurate than previous empirical models, and that predicts much higher RBE values in the distal edge of clinical proton beams.


Supplemental Note 2: Details of fitting
All of the model fitting was performed in MATLAB 2020, using the functions lsqnonlin, fminunc, fmincon, nlinfit and lsqcurvefit. We used lsqnonlin, fminunc and fmincon to determine the free parameter values when minimizing the L2 norm of the data since these functions allow for objective functions to be implemented with more sophistication than simply the sum-of-square distance between the data and the predictions. nlinfit and lsqcurvefit were used when optimizing the models to predict the RBE at the 2 Gy dose level when the residual sum-of-squares was to be minimized.
As our model (as well as the Mairani model) contain a large number of free parameters, the process of fitting such multivariate functions to our noisy training data often depended strongly on the initial conditions selected, since it is possible for the fitting algorithms to converge on local minimums within the parameter space rather than the global minimum. To work around this potential issue, for each fit, we performed 10000 iterations of each fit, using the previously-fit model parameters as initial conditions, but with a randomlygenerated amount of noise added to them (±10% of the initial conditions). This allowed for the fitting to be performed more broadly across the parameter space in an attempt to converge to the true global minimum distance between the predictions and the measured data. The set of model parameters resulting in the smallest L2 norm were then selected as the final parameter values.

Supplemental Note 3: Choice of functions to describe the LET dependence of the slope and intercept of the linear correlations between proton and x-ray radiosensitivity
When selecting functions to model for the slope and intercept's LET dependence, we sought functions using a minimum number of free number of parameters to produce the trends we observed and subject to a limited number of constraints. For the slope, we wished the function that described its LET dependence to tend towards non-negative values at high LET values. This is because a negative slope between, for example D10%,ion and D10%,x-ray implies that modulating a cell's radiosensitivity, for example via DNA-repair inhibition, would render cells more sensitive to x-rays but more resistant to protons which is nonsensical. We also wished that the slope would tend to positive values for LET values approaching zero for the same reason. For the intercept, we subjected it to the constraint that the intercept must be non-negative for all LET values. This is because for cells that are arbitrarily radiosensitive to x-rays, a negative intercept would imply that they have negative radiosensitivities to protons, which is non-sensical. These candidate functions are summarized in Table 1.
For every possible combination of slope and intercept function, we fit the set of free parameters to the training dataset's D5%, D10%, D20%, D37% and D50% values, and SF2Gy values using Matlab's lsqnonlin function, minimizing the relative distance between the predicted endpoint and the measured endpoint. We calculated the Bayesian Information Criteria 24 associated with each endpoint's fit to determine which parameterizations best reproduced the underlying trends in the data. Table 2 and Supplemental Table 2 -Supplemental Table 6 show these values. Generally, the simple exponential function was the best function describing the slope's LET dependence (smallest BIC values), but from these analyses alone it is not clear whether the intercept is best modeled by a constant, a linear function or an exponential function. Further comparisons combining multiple endpoints into a single predictive model suggest that the linear function with no offset offers a suitable description of the data when predicting α and β using multiple endpoints (Supplemental Note 5). Supplemental

Bayesian Information Criterion (BIC) value associated with the function predicting D10% Intercept
Slope

Bayesian Information Criterion (BIC) value associated with the function predicting D10% Intercept
Slope  acting on y provides b that minimizes the sum-of-squares distance between y and D b, the use of additional endpoints will tend to provide b that more faithfully reproduces = .
Consequently, given that the choice of endpoints to use in our formalism is arbitrary, we must decide which, and how many endpoints to use in our final parameterization. To do so, we must weigh how much more faithfully the model predicts the data using these additional parameters against the informational cost of including those extra parameters. To do this, we calculated the BIC associated with the different possible parameterizations of our model as follows: i) For each of D5%, D10%, D20%, D37%, D50% and SF2Gy, we calculated the free parameters ci, fi, and pi, associated with those endpoints ii) Then for each possible combination of endpoints, we predicted the αproton and βproton values of the whole database using Eqs. S9 and S10.
iii) We then calculated the relative distance between the predicted curves and the measured curves using the L2 norm normalized by the mean inactivation dose as follows: This metric was chosen since it incorporates information across the whole survival curve, and not being restricted to a particular dose or survival level, gives a description of how the models perform overall.

iv)
We then defined the maximum log-likelihood function, ln(L), in terms of the sum over the L2 norms of all the data (N datapoints) in the training dataset 24 : Then, for each parameterization, we calculated the BIC, where p is the number of parameters and n is number of survival curves in the training dataset, as:

vi)
We then compared the BIC values and selected the parameterization which minimized the BIC values vii) We also repeated this process using different parameterizations of the intercept, allowing its LET dependence to be described either by a constant, a linear function, a linear function with no vertical offset, or an exponential function with no vertical offset to clarify which parameterization, when used across multiple endpoints, yields the best description of αproton and βproton.

Supplemental Note 6: Parameterizations used for Wedenberg, McNamara and Mairani Models
In assessing the goodness of fit of the previously published models (Table 6 and Supplemental Note 5), we retrained them on our training dataset. In many of the original publications detailing these empirical models, a variety of LET-dependent functions were assessed to determine which function best predicted the training data. In some cases, we noted that the best parameterizations based on the goodness-of-fit metrics we calculated were different from those in the original publications or depended on the goodness-of-fit metric used. For example, the function that minimized the χ 2 /ν for the Mairani et al. model 4 was their fLQC (linearquadratic-cubic) function, and the function that minimized the BIC value was their fL (linear) function. This discrepancy arises largely because the BIC includes a larger penalty for the number of parameters used in the fit, thus favoring the less complicated linear function. In cases such as this, where different parameterizations could be justified among the previously published models, to give as much credit to these models as possible, we conservatively chose whichever parameterization minimized the goodness-of-fit metric investigated. Below we summarize the actual parameterizations of those models we used in this work.
For Wedenberg et al model, with free parameter a1, proton dose D, and LET L the function selected was: For McNamara et al model with free parameters a1, a2, a3 and a4, proton dose D and LET L, the function selected was: For Mairani et al model when minimizing the RRS for the predictions of the RBE at 2 Gy, we used the following function with free parameters a1, a2, a3 and a4, proton dose, D and LET, L: Where fLQC is a function of three additional free parameters, a5, a6 and a7, and the LET, L: For Mairani et al model when minimizing the BIC associated with the L2 norm, we used the following function with free parameters a1, a2, a3 and a4, proton dose, D and LET, L: Where fL is a function of one additional free parameter, a5, and the LET, L: The values of the free parameters found for these models are summarized in Supplemental Table 7-Supplemental Table 9 below. For reference, Supplemental imposing the same LET criteria on the database and retraining the models accordingly (Supplemental Table 12 and Supplemental Table 13), as well as without excluding any high LET data (Supplemental Table 11).
Supplemental Table 11: Goodness-of-Fit metrics for different models fit to our training dataset, minimizing either the RBE at a proton dose of 2 Gy or the L2 norm between the predicted and measured survival curves across all LET values up to 88.8 keV/ m. Regardless of how the higher LET data are handled, our model provides the lowest χ 2 /ν when fitting to the RBE at the 2 Gy dose level, and the smallest BIC values when minimizing the L2 norms. However, when the higher LET data are excluded, the Wedenberg model provides a slightly smaller BIC when minimizing the RSS around the RBE2Gy values. It is important to note here though that this that this does not suggest that the Wedenberg model is the most accurate around the 2 Gy dose level -our model is consistently the most accurate -but rather that when the higher LET data are not included, and one only wishes to predict the RBE at a single dose-level, 2 Gy, then the gain in accuracy our model (as well as those by McNamara and Mairani) provide may not be warranted given the use of the additional free parameters. However, given that the BIC associated with the L2 norms are considerably better for our model than any of the other models, regardless of the LET range investigated suggests that our model is the most robust across dose levels. It is also worth noting that even when the very high LET data are included, our model doesn't see such a drastic change in performance as the other models, further implying that our model is much more robust across a wide range of proton LET values.

Model
Supplemental Note 8: Linear correlations between proton and x-ray sensitivity for the parameters D5%, D20% and SF2Gy Fig. 1 shows the linear correlations between proton and x-ray radiosensitivity for the parameters D10%, D37% and D50%. Similar trends can be seen for the parameters D5%, D20% and SF2Gy as shown in Supplemental Figure   1 below.