Optimized I-values for use with the Bragg additivity rule and their impact on proton stopping power and range uncertainty

Novel imaging modalities can improve the estimation of patient elemental compositions for particle treatment planning. The mean excitation energy (I-value) is a main contributor to the proton range uncertainty. To minimize their impact on beam range errors and quantify their uncertainties, the currently used I-values proposed in 1982 are revisited. The study aims at proposing a new set of optimized elemental I-values for use with the Bragg additivity rule (BAR) and establishing uncertainties on the optimized I-values and the BAR. We optimize elemental I-values for the use in compounds based on measured material I-values. We gain a new set of elemental I-values and corresponding uncertainties, based on the experimental uncertainties and our uncertainty model. We evaluate uncertainties on I-values and relative stopping powers (RSP) of 70 human tissues, taking into account statistical correlations between tissues and water. The effect of new I-values on proton beam ranges is quantified using Monte Carlo simulations. Our elemental I-values describe measured material I-values with higher accuracy than ICRU-recommended I-values (RMSE: 6.17% (ICRU), 5.19% (this work)). Our uncertainty model estimates an uncertainty component from the BAR to 4.42%. Using our elemental I-values, we calculate the I-value of water as 78.73  ±  2.89 eV, being consistent with ICRU 90 (78  ±  2 eV). We observe uncertainties on tissue I-values between 1.82-3.38 eV, and RSP uncertainties between 0.002%–0.44%. With transport simulations of a proton beam in human tissues, we observe range uncertainties between 0.31% and 0.47%, as compared to current estimates of 1.5%. We propose a set of elemental I-values well suited for human tissues in combination with the BAR. Our model establishes uncertainties on elemental I-values and the BAR, enabling to quantify uncertainties on tissue I-values, RSP as well as particle range. This work is particularly relevant for Monte Carlo simulations where the interaction probabilities are reconstructed from elemental compositions.

dose calculation techniques (Schneider et al 2000). Since SECT does not provide enough information to predict these parameters accurately, uncertainties are introduced (Schaffner and Pedroni 1998). These uncertainties are taken into account during treatment planning by adding uncertainty margins to the treatment volume.
In recent years, dual-energy CT (DECT) was explored by several groups. DECT provides additional information and therefore has the potential to increase accuracy in the prediction of tissue parameters (Alvarez and Macovski 1976, Torikoshi et al 2003, Williamson et al 2006, Bazalova et al 2008. Different formalisms exist to either predict RSP values (Yang et al 2010, Hünemohr et al 2014a, Bourque et al 2014, Möhler et al 2016, Taasti et al 2016 or elemental compositions (Landry et al 2013, Hünemohr et al 2014b, Lalonde and Bouchard 2016 from a DECT scan. Some of these formalisms have been successfully validated in animal tissue studies (Taasti et al 2017, Bär et al 2018, Möhler et al 2018, Xie et al 2018, and it was shown that uncertainties on the beam range prediction can be reduced by 0.34% using DECT , Bär et al 2018. With DECT being on the edge of clinical implementation for radiotherapy, one major remaining source of uncertainty lies in the determination of the mean excitation energy, or I-value, of patient tissues. The portion of range uncertainties arising from I-values was previously estimated to 1.5% (Paganetti 2012). A recent survey (Taasti et al 2018) shows that particle therapy facilities apply a relative range margin or 3.5%, which was recommended by Paganetti and includes the before mentioned 1.5% portion coming from the I-values.
In practice, the I-values for compounds are calculated from elemental I-values using the Bragg additivity rule (BAR). Current clinically used elemental I-values were estimated by Berger andSeltzer in 1982 (Berger andSeltzer 1982), and in 1984 those values were adapted as recommendation in the ICRU report 37 (ICRU 1984), and later taken over for the use in proton and ion radiotherapy in ICRU report 49 (ICRU 1993). Berger and Seltzer (1982) used a large set of different compound I-values to estimate the elemental I-values for elements in compounds. The compound I-values they used were taken from previous publications using two different methods. Firstly, they took measured stopping power data (Thompson 1952, Tschalär and Bichsel 1968, Nordin and Henkelman 1979, Bichsel and Hilko 1981 from which they derived the compound I-values and their uncertainties according to equations 3.9 and 3.11 in ICRU 37 (ICRU 1984). Secondly, they took calculated compound I-values from dipole oscillator-strength distributions or dielectric-response functions (Bader et al 1956, Zeiss et al 1977, Thomas and Meath 1977, Ashley et al 1978, Painter et al 1980, Jhanwar et al 1981. Details of how those I-values were derived from the oscillator-strength and dielectric data can be found in ICRU 37. Based on those compound data, they estimated two sets of elemental I-values to use in compounds in combination with the BAR, one set for gases and one set for liquids and solids. It is important to note that the I-value of an element or a molecule depends on whether it is unbound or bound, and the type of chemical bond. Hence, elemental I-values for the use in compounds are different from I-values of unbound elements. Since the elemental I-values of Berger and Seltzer were recommended by the ICRU, they are in clinical use to calculate the I-values of compound materials such as human body tissues. The estimation of range uncertainties arising from I-values has always been challenging. Andreo (2009) showed differences in beam ranges of 0.3 g cm −2 for a 122 MeV proton beam when the water I-value varies from 67 eV to 80 eV, covering the variety of values proposed in literature. Differences get larger when considering different tissue types or different particle species. Besemer et al (2013) performed a variation study that uniformly varies the tissue I-values, and evaluated the influence on patient dose distributions. They showed that a 10% variation of I-values influences the R 80 beam range by up to 4.8 mm, and resulting dose distributions by up to 3.5%. Although the uncertainties on I-values can be relatively high, two studies by Yang et al (2012) and De Smet et al (2018) suggest that the resulting RSP values are much lower since correlations between water and medium need to be taken into account. Another study by Doolan et al (2016) investigated the influence of different correction terms to the Bethe formula on the calculated stopping power. They suggest to use the I-value as a free parameter to optimize according to which corrections to the Bethe formula are used, in order to avoid systematic errors on RSP values. Recent work suggests the estimation of a patient specific tissue I-value from MRI imaging (Sudhyadhom 2017).
It was shown in recent studies that accurate knowledge of the I-values allows for reduced uncertainties of the water-to-air stopping power ratio of carbon ion beams, an important quantity in reference dosimetry (Sánchez-Parcerisa et al 2012). Furthermore, I-values are one of the main sources of uncertainties in water equivalent range calculations and beam transport models for protons and heavier ions (Zhang et al 2010, Paul andSánchez-Parcerisa 2013).
The aim of this work is to revise the currently used elemental I-values. We believe that a revision of I-values and an addition to the pioneer work by Berger and Seltzer is necessary for the following reasons: 1. The original publication by Berger and Seltzer does not give details on how exactly the elemental I-values are derived from the given compound measurements, making it difficult to reproduce their values. 2. Since the work of Berger and Seltzer, several new stopping power measurements were performed (Bichsel and Hiraoka 1992, Hiraoka et al 1993, 1994, Bichsel et al 2000, Kumazaki et al 2007 which can be included into the estimation of elemental I-values for the use in compounds.
3. The values by Berger and Seltzer are quoted without uncertainty budget, which makes it difficult to estimate resulting uncertainties on RSPs and ranges.
In order to revise the currently applied I-values, we develop a mathematical model to find, based on old and new measurement data, an optimal set of elemental I-values for the use in compounds. Furthermore, our model establishes an uncertainty budget on our newly found set of elemental I-values as well as on the BAR. Our uncertainty budget allows the propagation of uncertainties from elemental I-values to relative stopping powers of tissues and ultimately to beam ranges, to give a rigorous estimate of range uncertainties arising from I-values.

Optimal elemental I-values to estimate compound I-values
The stopping power of charged particles in a medium m is described by the Bethe formula (Bethe 1930). We can formulate the RSP in terms of relative electron density ρ e and stopping number L as where L med and L w are the stopping numbers of the medium and water, respectively. The stopping number of an arbitrary medium is expressed as with m e the electron mass, c the speed of light, β is the velocity in units of c, and I is the mean excitation energy of the medium. Please note that the Bethe formula is valid for protons and heavier ions, and that the I-value is a tissue-specific parameter and does not depend on the particle type. The herein presented values are hence valid for electrons, protons and heavier ions. Let us define a series of M media indexed by i = 1, . . . , M consisting of N elements and with given elemental weights w med,ij , with j = 1, . . . , N. The Bragg additivity rule (BAR) allows an estimation of the mean excitation energy I med,i of the ith medium using the weighted sum of the logarithmic elemental mean excitation energies: where λ med,ij is the fraction of electrons from the jth element in the ith medium and given by with Z j and A j the atomic number and molar mass of the jth element, and Z A med,i is the number of electrons per unit mass in the medium (in mol/g). Using matrix notation, the BAR can be written as the following estimator where y med is a M × 1-dimensional array containing the logarithm of experimental I-values and ŷ elem is an array of dimension N × 1 containing the optimized logarithm of elemental I-values for use with the BAR defined aŝ The matrix Λ med of dimension M × N contains the fractions of electrons for the respective materials and its elements are written as λ med,ij , corresponding tor the ith medium and jth element. We propose to determine a new set of optimized elemental I-values, i.e. Î j , by finding the weighted least square solution of equation (5). To take measurement and model uncertainties into account, we introduce weighting factors accounting for uncertainties: with u med,i being the relative uncertainty of the I-value experimental measurement of the ith medium. Note that because these uncertainties represent the absolute uncertainty of the natural logarithm of the I-value, they are reported in relative uncertainty on the I-value (i.e. in %). These weighting factors multiply individually both sides of the equation system (5) to account for uncertainties, leading to a new equation system where the elements of ỹ med are ω i y med,i and the elements of Λ med are ω i λ med,ij . We can now find the least square solution to equation (8): with y elem being the estimation of the optimized logarithmic elemental I-values and M a projection matrix (from the measurement to the solution) defined to ease the notation. To find the uncertainties on elemental I-values, we construct the covariance matrix of y elem as follows: with V (ỹ med ) being the covariance matrix on measured material I-values with each line weighted by its corresponding ω i . Note that because the measurements are assumed independent, V (ỹ med ) is diagonal. Equation (10) is defined by combining two terms. The first one is obtained by applying the rule of uncertainty propagation on equation (9), since the Jacobian (∂ y elem /∂ỹ med ) is the projection matrix M . The second term is added to account for the model uncertainty. Indeed, the rule of uncertainty propagation can only yield accurate uncertainty estimations if the model is exact. Because the BAR is not completely accurate, it is judicious to add a model uncertainty component u 2 BAR affecting each optimized values individually and in an independent manner. This way, we divide the uncertainties involved in the estimation of elemental I-values into experimental type A uncertainties and model-related type B uncertainties. The resulting V ( y elem ) is a non-diagonal square matrix of dimensions N × N accounting for statistical correlations in the solution.
The solution expressed in equation (9) using indirectly measured (from stopping power measurements) and calculated (from dipole-oscillator and dielectric data) compound I-value data yields a new set of elemental I-values Î j for use with the BAR. The required I-values were taken from different sources as listed in table 1. They include the data provided in ICRU report 37, table 5.3 of ICRU (1984) and more recent publications. The materials, elemental compositions and corresponding I-values can be found in tables A1 and A2. This formalism starts with compound I-values, as extracted from stopping power measurements, dipole-oscillator and dielectric data. Like Berger and Seltzer, we divide the data into two groups: (1) gases; (2) liquids and solids. For gases, the data used in this work are the same than the ones used by Berger and Seltzer (1982) to determine the recommended elemental I-values in ICRU report 37. For liquids and solids, we added data published in recent literature (see table 1, numbers 11-15). In total, we use 74 liquids and solids for calibration, including six different I-values for water (Thompson 1952, Nordin and Henkelman 1979, Bichsel and Hiraoka 1992, Bichsel et al 2000, Kumazaki et al 2007, Emfietzoglou et al 2009. We refer to this method as a calibration because calculated or measured compound I-values are used to derive model parameters (elemental I-values) that will be used for the prediction of unknown values, such as I-values in patient tissues. The elemental I-values serve as model parameters. We use the set of data given in table 1 to derive those model parameters, however, the initial set could be different or expanded. We obtain two sets of elemental I-values which are optimized for the use in compounds (one set for gases, one for liquids and solids) in combination with the BAR. We use equation (10) to report uncertainty values u elem for the newly determined I elem . The here reported uncertainties are standard uncertainties (68% confidence interval). Some of the data used for our analysis, especially the ones used by Berger and Seltzer, are quoted for a 90% confidence interval, as reported in ICRU report 37, footnote 10. Whenever this was the case, the uncertainties were divided by a factor of 1.6 to convert from the 90% to the 68% confidence interval.
To test the validity of our data, we perform a self-consistency test. For this, we use the optimized elemental I-values to reproduce the calibration data set. The root mean square (RMS) errors between actual and predicted calibration data are compared with predicted data using ICRU-recommended elemental I-values. Again, we separate gases from liquids and solids.

Estimation of u BAR
While u med,i can be derived from experimental uncertainties, u BAR needs to be estimated using a model. To estimate u BAR , we use the calibration data set as listed in table 1, and calculate the residual error r between the estimated and experimentally measured values for y med : with K = ΛM and where the elements of K equal the ones of K divided by the weighting factors, i.e. K ij = 1 wiK ij . We can now estimate the covariance matrix V(r) as To solve for u 2 BAR , we find the value such that the sum of the residuals squared normalized to their variance equals its number of degrees of freedom, that is: Note that the approach is based on the assumption that the experimental data for a particular element follows a Gaussian distribution. The resulting equation (13) follows a chi-square distribution with expectation value equaling its number of degrees of freedom N − M, N being the number of experimental data and M the number of optimized elemental I-values.

Application of optimal elemental I-values to water and reference human tissues
Using the optimized set of elemental I-values, we determine compound I-values for water and a set of 70 human reference tissues White 1986, White et al 1987). For the I-value of water, we need to differentiate between three different values: I w,BAR represents the I-value of water calculated with the elemental I-values quoted in this paper and the BAR, I w,ICRU37,BAR represents the I-value of water calculated with the ICRU 37 recommended elemental I-values and the BAR, and I w,ICRU90 represents the value recently recommended by ICRU 90. The compound I-values are compared to the results obtained with ICRU 37 recommended values. We establish the uncertainties on compound I-values using the covariance matrix of the elemental I-values V ( y elem ):

Uncertainties on relative stopping powers
Once the uncertainties on the mean excitation energies are determined, it is possible to propagate these into stopping power uncertainties. In this way, we quantify the uncertainty on the RSP of medium to water originating from uncertainties on I-values. Since the uncertainties of medium and water can be correlated depending on the water content of the medium (Yang et al 2012, De Smet et al 2018, it is important to consider covariances. Combining equations (2) and (3), the stopping number can be expressed as The derivative of the stopping number with respect to ln I i is then found to be Table 1. Literature used to retrieve the I-values of compounds, by either dipole oscillator-strength distributions, dielectric-response functions, or measurements of the energy loss. Numbers 1-10 were used by Berger and Seltzer (1982)  We can now express the derivative of the RSP with respect to ln I i as The variance on the RSP can now be written as using the following rule

Uncertainties on beam ranges
To quantify the impact on beam ranges, we perform Monte Carlo transport simulations of a pristine proton beam in homogeneous media (volume: 30 × 30 × 30 cm 3 ). We score the energy loss and position of each interaction of the beam with the medium. We choose water and five different human reference tissues (Adipose 3, skeletal muscle 1, brain white matter, femur whole, cortical bone) relevant to proton therapy. For every material, four simulations are performed: (1) using ICRU-recommended I-values; (2) using our suggested I-values; (3) using the upper uncertainty limit and (4) using the lower uncertainty limit. For the simulations, we use the Geant4 code (Version 10.03.p02) with the QBBC physics package (Ivantchenko et al 2012). To calculate the energy loss in a medium, Geant4 uses the restricted Bethe formula with shell, density and higher order correction terms (i.e. the restricted Bethe-Bloch equation). We simulate proton beams of 173 MeV using 10 6 particles per beam and a 1 mm cut-off value for secondary particles.

Optimal elemental I-values to estimate compounds I-values
The proposed approach results in a set of optimized elemental I-values for the use with the BAR with compounds. Using the same measured compound I-values than Berger and Seltzer and more recent literature on measured I-values, we calculate optimized elemental I-values for the use in gases and for the use in liquids and solids separately. Our optimized elemental I-values differ from the ones suggested by Berger and Seltzer, as tabulated in tables 2 and 3. Table 4 shows the correlation coefficients of the optimized elemental I-values for liquids and solids. We use both sets of elemental I-values to perform a self-consistency test on the calibration data. Using the ICRU 37 recommended elemental I-values suggested by Berger and Seltzer, we observe RMS errors of 1.02% (gases) and 6.17% (liquids and solids) when using the BAR to predict the underlying experimental I-values. Using our optimized elemental I-values, this prediction error can be reduced to 0.05% (gases) and 5.19% (liquids and solids). The model uncertainty arising from the BAR , i.e. u BAR , is quantified as 4.42%.

Application of optimal elemental I-values to water and reference human tissues
Using our method, we estimate the compound I-value of water to I w,BAR = 78.73 ± 2.89 eV. This value is in good agreement with the recent recommendation given in ICRU 90 (ICRU 2014), which is based on the value I w,ICRU90 = 78 ± 2 eV given in Andreo et al (2013

Uncertainties on relative stopping powers
We use our optimized I-values to calculate RSP values for 70 human reference tissues. Figure 2 shows the uncertainties on RSP values of 70 human reference tissues, arising from uncertainties on I-values only. We observe uncertainties on RSP values between 0.002% (mammary gland) and 0.44% (adipose tissue 3). The uncertainties observed are the smallest for the soft tissues since their water content is the highest between adipose tissue, soft tissues and bones. Figure 3 shows the percentage depth dose (PDD) curves for water and five human reference tissues, each using four different sets of I-values: (1) the ICRU-recommended values;

Uncertainties on beam ranges
(2) the optimized I-values resulting from this work, (3) the optimized I-values resulting from this work plus 1 standard deviation and (4) the optimized I-values resulting from this work minus 1 standard deviation. We observe differences in the range of the distal 80% of the maximum dose (R 80 ) between ICRU-recommended values and our values of 0.75 mm (adipose tissue 3)-1.10 mm (water using I w,ICRU37,BAR ). We find range uncertainties between 0.31% and 0.47%, with the lowest uncertainty found in femur tissue, while the highest uncertainty is found in water (see table 6).

Discussion
In this work, we investigate RSP and range uncertainties arising from mean excitation energies. We establish a mathematical model to optimize elemental I-values for the use in gases and liquids and solids with the BAR. To calculate our optimized I-values and establish an uncertainty budget, we utilize I-value and stopping power measurements from literature, most of which were used by Berger and Seltzer to establish the ICRU 37 recommended values, however we also include more recent measurements. The set of optimized elemental I-values for the use with the BAR in compounds and the reported uncertainties that can be used to accurately assess the uncertainties on     tissue I-values. Furthermore, we provide an estimation of the uncertainty coming from the BAR itself, and include this uncertainty in our model. Our model allows the propagation of uncertainties to RSP values and beam ranges, providing a better understanding of the resulting uncertainties in proton therapy treatment planning. This study focuses on the effect of the I-value of a medium on the range of a particle beam and associated uncertainties. The magnitude of this effect depends on various factors such as the dose calculation algorithm and the chosen modality to determine the treatment planning inputs. A higher precision in I-values and their uncertainties leads to a higher precision in range prediction for both, analytical and Monte Carlo based dose calculation approaches since in both cases I-values are needed to calculate the stopping power using the Bethe formula.
Elemental I-values of liquids and solids are of special interest for proton therapy treatment planning. We propose a set of values which differs from the values recommended in ICRU 37. This is due to the fact that the underlying model used to find the optimized I-values differs from the methods used by Berger and Seltzer to determine the originally proposed values. For the majority of our proposed values (C, N, O, Al, P, Cl), the ICRU-recommended values are within the uncertainty budget given by our model. However, few exceptions are observed. We find considerably higher elemental I-values for the elements H, F, and Ca than Berger and Seltzer. In turn, our values for Si is lower. While the optimized value for P is close to the ICRU recommended value (195.5 eV versus 199.93 eV), we observe a high uncertainty of 42.45 eV. The high uncertainty can be explained due to the lack of data available for P. From the 74 liquids and solids, only 3 materials contain traces of P. The high observed uncertainty however is of little concern to clinical applications. First of all, P is only abundant as trace element in the human body. Secondly, it is only observed in bones, where we observe a strong statistical anti-correlation with Ca, compensating for the high uncertainty in P. The calibration materials used herein are taken from various sources of literature, most of them were already utilized by Berger and Seltzer to recommend the currently clinically applied elemental I-values. As a self consistency study, we evaluate the accuracy of our optimized elemental I-values in comparison to ICRU-recommended elemental I-values to predict the I-values of the calibration material. Our optimized elemental I-values show a lower RMS prediction error than the ICRU-recommended values, indicating that our values are well suited to predict the I-values of tissues when the BAR is used. When calculating the I-values of 70 human reference tissues, we observe generally larger tissue I-values using the optimized elemental I-values. Using our technique and optim ized elemental I-values, we achieve an I-value for water of I w,BAR = 78.73 ± 2.89 eV, which is in good agreement with the value recommended in the the recently published ICRU 90 report. Using the elemental I-values recommended in ICRU 37, the I-value for water is estimated as I w,ICRU37,BAR = 75.32 eV. Doolan et al (2016) tested different sets of elemental I-values derived with corresponding corrections to the Bethe formula. They concluded that the accuracy of the RSP values with and without corrections is equivalent as long as the correct combinations of I-values and corresponding corrections is used. Here, we used the Bethe form ula without corrections to calculate the RSP uncertainties. If one wishes to include corrections to the formula, this could be easily incorporated into equation (2), but it is not expected to largely influence the herein quoted RSP and range uncertainties.
Our mathematical model incorporates estimates for uncertainty values associated with each optimized elemental I-value, and allows the propagation of uncertainties on tissue I-value uncertainties. We show the importance of taking into account statistical correlations between uncertainties of the different elemental I-values. Our analysis shows that if statistical correlations are ignored, the uncertainties might be overestimated by up to 0.5%.
The method proposed herein allows the estimation of uncertainties on tissue RSPs resulting from I-value uncertainties. We observe the highest RSP uncertainties in adipose tissues and bones, which is expected as those are the tissues with a low water content. Those findings were already discussed in previous studies by Yang et al (2012) and De Smet et al (2018). We observe low RSP uncertainties in soft tissues, which can be attributed to their high water content.
The resulting beam range uncertainties are assessed in water and five selected human reference tissues. We choose the tissues according to their importance to radiotherapy and abundance in the human body. Observed range uncertainties are between 0.31% and 0.47% for a 173 MeV proton beam. Overall, the ranges predicted using our optimized elemental I-values are systematically larger than the ranges calculated based on ICRU-recommended I-values. This is expected since our optimized I-values calculate larger tissue I-values than the ICRU recommendation. In our study, the observed RSP and range uncertainties are much lower than the currently assumed 1.5%.
It should be noted that the systematic shift in depth dose curves observed in the Monte Carlo simulations will not necessarily be observed with more simplistic dose calculation algorithms as implemented in treatment planning systems, since they usually use the stopping power relative to water for dose estimation. However it should be emphasized that our model yields a more realistic I-value for water (and closer to the recent ICRU 90 recommendation) than ICRU 37. Future work will focus on the application of optimized I-values to CT scans of tissues, which can potentially proof the validity and superiority of one set of elemental I-values over the other.

Conclusion
We propose a new set of optimized elemental I-values for the use with the Bragg additivity rule in compounds. Our mathematical model establishes an uncertainty budget on elemental I-values that can be propagated to compound I-values, accounting for experimental uncertainties as well as the uncertainty on the Bragg additivity rule itself. With our model, we provide realistic uncertainties estimations on proton RSP and beam range values in human tissues. The herein presented data show that the currently assumed range uncertainty originating from I-values may be overestimated. Proposed elemental I-values and corresponding uncertainty budgets provide evidence for a reassessment and possible reduction of those range uncertainties in clinical routine. Table 6. Calculated beam ranges in terms of R 80 (in mm) using Monte Carlo proton beam transport simulations. The uncertainties reported are resulting from the uncertainties on our optimized I-values and the differences are taken between ranges simulated with ICRUrecommended I-values and ranges simulated with our optimized I-values.

Acknowledgments
This work is funded by EPSRC UK (case studentship No. 20873) and the National Physical Laboratory (UK).