An approximate analytical solution of the Bethe equation for charged particles in the radiotherapeutic energy range

Charged particles such as protons and carbon ions are an increasingly important tool in radiotherapy. There are however unresolved physics issues impeding optimal implementation, including estimation of dose deposition in non-homogeneous tissue, an essential aspect of treatment optimization. Monte Carlo (MC) methods can be employed to estimate radiation profile, and whilst powerful, these are computationally expensive, limiting practicality. In this work, we start from fundamental physics in the form of the Bethe equation to yield a novel approximate analytical solution for particle range, energy and linear energy transfer (LET). The solution is given in terms of the exponential integral function with relativistic co-ordinate transform, allowing application at radiotherapeutic energy levels (50–350 MeV protons, 100–600 Mev/a.m.u carbon ions). Model results agreed closely for protons and carbon-ions (mean error within ≈1%) of literature values. Agreement was high along particle track, with some discrepancy manifesting at track-end. The model presented has applications within a charged particle radiotherapy optimization framework as a rapid method for dose and LET estimation, capable of accounting for heterogeneity in electron density and ionization potential.

Quite aside from the technical challenges, there are also a number of theoretical physical and modeling problems that demand the attention of physicists and modellers. One of the most significant difficulties is adequately modeling tissue-particle interaction, and the computational power required for this. Modeling of radiation interaction is crucial for many applications, and if charged particles are to be exploited to their full potential it is imperative that energy deposition through different media is well understood. Monte Carlo (MC) techniques are typically employed for this purpose 13 . MC particle transport packages such as GEANT4 14 , MCNPX 15 and FLUKA 16 are extremely powerful and capable of simulating the interaction of particle radiation with matter as it traverses a medium. These approaches can capture system behaviour at all scales, including secondary and higher order interactions, but are computationally expensive with non-trivial evaluation times. Fast, approximate methods for estimating the energy profile through different media are especially useful in radiotherapy treatment planning, where multiple rapid dose calculations are needed to enable interactive optimization. More than this, such methods should ideally work with various particle energies and types, and be applicable not only to uniform media but also to the heterogenous nature of human tissue.
To achieve this, the purpose of this work is to derive an analytical model which captures the general behaviour of charged particles at the energy ranges pertinent to therapy  MeV for protons, 100-600 MeV/a.m.u for carbon ions). Previous authors have deduced useful functions to describe the Bragg curve in radiotherapy, based on empirically-determined range-energy relationships 17,18 which work well for low energy protons, but these are typically restricted to dose. Other formulations include empirically derived approximations of other radio-biologically relevant parameters, such as range or Linear Energy Transfer (LET) for a specific particle radiation type [19][20][21] . In this work, we take a different approach, starting from fundamental physics in the form of the Bethe equation to yield an analytical general solution for velocity, range, energy, and LET. This is broadly applicable for most charged particles used in therapy, with the notable exception of electrons, which violate important assumptions in the Bethe equation and are not considered in this work. The formulation here is suitable for heavier charged particles across the range of therapeutic energies. We compare results from this work with tabulated data, and MC simulations for protons, to ascertain how effective such an approach is at capturing the gross behaviour of such systems, and demonstrate potential application of the method outlined.

Theoretical background
Charged particles undergo electronic interactions as they traverse a medium, losing energy in the process and slowing down. As the particle slows down, the density of ionizations induced in the medium increases, before. Ionization density drops abruptly to zero beyond the Bragg peak, since all of the particle's kinetic energy has been exhausted and it can be considered stationary. Under certain assumptions, the stopping power in a medium is given by the relativistic Bethe equation 22 . Further corrections can be made to account for high energy particles and high atomic number targets, including the Barkas, Bloch and Fermi corrections, comprehensively reviewed by Ziegler 23 . However, at energies and in materials relevant to radiotherapy, these effects are expected to introduce errors less than 1% (less than 1 mm) whereas in radiotherapy, uncertainties greatly exceed this level 13 . Accordingly we neglect these corrections in order to achieve an analytical solution, without loss of accuracy in the regime of interest. The relevant form of the Bethe equation is then where n is the electron density of the material, e the electron charge, m e the electron mass, I the mean excitation potential, z the multiple of electron charge and β = v c , where v is the speed of the particle and c is the speed of light in a vacuum. Particle energy is a function of particle speed, itself a function of distance. The quantity given by equation 1 can be considered the linear electronic stopping power of a medium for a given particle type, and represents the effect of the medium on the particle. The reciprocal effect of the particle on the medium is termed unrestricted linear energy transfer (LET); the two quantities are equal by the definitions in ICRU Report 85 24 . In general LET only considers collisions where the kinetic energy imparted to secondary electrons is below a given threshold (or restriction). This restricts the quantity to shorter-range electrons, which might provide better characterisation of radiation effects on the cellular scale. However, this study intends to characterise radiation on a millimetre scale, and microscopic dosimetry and radiobiological effects are considered outside its scope. It therefore follows many other medical physics investigations [25][26][27] in considering unrestricted LET a relevant quantity for clinical modelling work. This is justifiable since calculations suggest a monotonic relationship between restricted and unrestricted LET 28 , and a number of radiobiological models 29, 30 make use of unrestricted LET.
The range of a charged particle in a medium may be defined in a number of different ways, depending on the theoretical or experimental context. In this work, we focus on the continuous-slowing-down approximation (CSDA) range, obtained by integrating the inverse of Equation 1 between a particle's initial energy and its stationary state. It should be noted that in practice the projected range (ie. the distance travelled by the particle along its initial direction) is slightly reduced by small deflections caused by scattering. However, the impact of this is minimal at clinical energies in low Z materials, with the 'detour factor' of the projected range compared to the CSDA range being greater than 99.87% for protons of energy >100 MeV in water 31 . Our definition of range is therefore, expressed in both energy and velocity terms as

Model derivation
We begin from the relativistic definition of particle kinetic energy, expressed in terms of the Lorentz factor: it can be shown that If equation 4 and equation 1 are equated directly, a solution for v(x) is analytically intractable. To circumvent this, we introduce a co-ordinate system x′ where the classical expression for particle kinetic energy of = holds. Thus, we can write It is possible to simplify equation 1 by reducing the β 2 term, yielding

This is justified if
. For a high energy proton at 250 MeV, there is only ≈1% difference between both sides, and ≈2% for high energy (400 MeV/a.m.u) carbon ions. This minor difference rapidly dissipates as the particles decelerate, and so the simplification can be employed with only minor discrepancy for particles in the radio-therapeutic energy range. For brevity, we let = This can be solved, and yields a solution in terms of the exponential integral function Ei. From this, it can thus be shown that particle range in the x′ coordinate system is given by If we change the lower limits on equation 7, we can readily show that , and thus the solution for velocity v(x′) in terms of x′ is simply The velocity formula in equation 10 is sufficient when γ  1, but in practice high-energy charged particles used in treatment often are far beyond this, and a better expression is required to account for them. To allow for relativistic effects, we can introduce a transform from x′ to x. We can re-arrange equations 4 and 5 and equate them; it can thereby be shown that dx = γ 3 dx′ and consequently that This transform can be easily evaluated provided once v(x′) is known, and the transform yields the true solution v(x). In the x frame, the true CSDA range, R, is then given by This can also be calculated using a series expansion method, outlined in the appendix (Supplementary data). An illustration of the relativistic transform is depicted in Fig. 1. Once the transformed velocity profile is known, it can readily be applied to yield the mono-energetic energy curve and LET through identities 3 and 6 respectively. The inverse exponential integral function -Ei −1 . The solution derived above requires knowledge of the inverse exponential integral function. This is not a standard function, and to the authors' knowledge no closed form expression exists for it. Values can be readily derived using simple root-finding methods in most computational software packages. It is also possible to employ some simple approximations -for most instances the approximation Ei −1 (x) ≈ ln(xln(xlnx)) can be used without introducing serious error, and Chebyshev polynomials can be employed to higher precision 33 , though knowledge of the limits of validity are important. A series expansion approach is also possible. For this work, we wished to plot the solution at high resolution (steps of 1 μm), and to speed this up we created a high precision look-up table using root-finding methods. The look up table, as well as examples of how to implement the model, are available in the Supplementary data.

Methods
Comparison with numerical methods. To ascertain how well the solution proposed fits the both forms of the Bethe equation, numerical solutions for equation 1 (fully relativistic) and 6 (simplified form) were found using a Runge-Kutta method. These solutions were contrasted with the analytical method outlined in this work, and the level of agreement between the two approaches was quantified. Physical parameters used in this work are given in Table 1.
Comparison with tabulated data, Monte Carlo methods and experimental data. MC particle transport codes are traditionally used as the primary method for proton dose-depth calculations. We used the MCNPX v2.7.0 package 15 to simulate proton pencil beams with a range of initial energies in water-like tissue and contrasted the predictions with the model outlined here. There has been less work done on carbon ions as these are less common in clinical use, though MC techniques using e.g. GEANT4 are emerging to explore this promising modality 34 . As implementing MC for carbon ions is still an area of active research, we have instead taken the measured range of carbon ions in water from various experimental investigations and contrasted this with the estimated range using the model to gauge the validity of model predictions for heavy ions like carbon.
The proton MC geometry comprised an infinitely narrow pencil beam of monoenergetic particles, incident on a water phantom with length 20% greater than the proton range. Phantoms of two different lateral cross sections were considered: a 'broad' 10 cm × 10 cm case, in which protons undergoing multiple Coulomb scattering (MCS) would be expected to stop, and a 'narrow' 1 mm × 1 mm case, in which the majority of protons undergoing MCS would be expected to escape. Mean energy and LET profiles were recorded at 1 mm intervals using surface tallies (MCNPX type F2, binned into 1000 subdivisions using the E card, modified to tally unrestricted stopping power by including the LET flag on an FT card). Only protons were tracked, from the source energy down to 1 keV. Default physics options were employed, with the addition of light-ion recoil. In addition to this, range estimates and v(x) (note Velocity shown on the vertical axis as a fraction of the speed of light c) for a 250 MeV proton. In this example, the projected range of the proton is 38.15 cm when the relativistic transform is considered. If this were neglected, projected range would be only 24.36 cm. This substantial difference suggests that taking account of the relativistic transform is vital for particles in the radiotherapy energy range. were compared to fitted data from the National Institute of Standards and Technology (NIST) PSTAR database for protons in liquid water 35 and compared to results derived from this method. Carbon ion ranges were compared to tabulated results from ICRU report 73 36 and the MSTAR 37 software empirical estimation.
Model applications. The chief advantage of charged particle therapy is that dose can be accurately delivered to a tumour region whilst largely sparing healthy tissue downstream, in contrast to photons which deposit dose relatively uniformly throughout the beam path. However, this advantage also presents a significant difficultly -dose must be well-aimed and delivered to very high precision 38 . If this is not the case, then the bulk of dose deposition could occur away from the tumour and into healthy organs, negating any benefit from therapy. To complicate matters, for example there is a six-fold difference between the electron density of inflated lung and cortical bone 39 -further examples are given in Table 2. Contemporary fast dose calculation methods employ pencil beam algorithms, which account for heterogeneity by scaling dose-depth curves measured in water; such an approach does not provide estimates of LET, which is an important consideration for the radiobiological evaluation of treatment approaches 40 . The model here could be useful for such applications, as it can be applied to quickly estimate energy distribution through inhomogeneous media, such as tissue and we demonstrate this here by applying the model to inhomogeneous tissue derived from a typical CT scan as proof of principle. Additionally, the model here can be used to explore the effect of the physical parameters on energy deposition profile. Most of the parameters required in this model are accurately known physical constants, as can be seen in Table 1. Of all these parameters however, there has been considerable variation in literature for the mean ionization potential for water. The value used of 75 eV is in line with the International Commission on Radiation Units and Measurement (ICRU) guidelines 36 , with a more recent ICRU report suggesting a lower value of 67 eV 41 . Other authors suggesting a much higher value of 80.2 ± 2 eV 42 based on experimental range measurements. Using the model here, we can readily vary this and see what impact this has on range and energy profiles obtained.

Results
Comparison with numerical methods. Equations     Comparison with tabulated data, Monte Carlo methods and experimental data. Figure 2 depicts the simulated energy profiles from the model established in this work and from traditional MC methods. In general there is good agreement between both approaches, quantified by the data in Table 3. The coefficient of determination column, R 2 compares the fit of the MC data with the model data, and predicted ranges, R T , are shown also. For the model, this was theoretically predicted by equation 12, and for MC simulations this was the depth beyond which no protons penetrated in the 'broad' phantom. While the ranges are in close agreement, there are some interesting differences. Most noticeable is that the MC simulation predicts a slightly lower energy along the curve, as quantified in the table. Secondly, it predicts a small tail that extends beyond the predicted R T of the model. This is considered in the discussion, but most likely arises from the fact that a solution of the Bethe equation doesn't consider energy and path length straggling, arising from stochastic effects and MCS respectively, or secondary charge particle production -all of which are incorporated in MCNPX. Further evidence for this is seen in comparing the curves in Fig. 3 for the 'broad' MC phantom to those for the 'narrow' phantom, where scatter effects are lower and the discrepancies between the model and simulation lesser. Table 3 also shows the percentage error between range estimated by method in this work and PSTAR NIST 35 data, with errors of ≪1% between database and model. The tail discrepancy is also obvious in the LET curves, depicted in Fig. 4. Obeying the Bethe equation, the model predicts a smooth asymptotic curve by definition from equations 1 and 3. The model closely matches the MC LET for most of the particle range, but diverges close to the terminal end of the track, as can be seen in the   figure. This difference may be explained by straggling and secondary interactions, which the Bethe equation does not consider a priori. Table 4 shows tabulated and experimentally determined ranges for carbon ions using different methodologies. For ease of comparison, the experimental range was determined as the depth of the zenith of the Bragg peak. The predicted range from the model agrees within ≈1% of empirical values derived from MSTAR 37 , and also well with the experimentally determined values across all energy ranges between 135 MeV/a.m.u to 380.45 MeV/a.m.u with   relative errors ranging from 0.21 % to 6.14 %. Part of this discrepancy might be apparent in 2(b), in which a small reduction in particle energy at a given depth is observed when using the full Bethe equation, compared to the simplified form. A further part is likely attributed to straggling and secondary interactions. Figure 5 depicts the model prediction for range versus the tabulated values in ICRU 73 36 , showing high agreement for carbon-ions in the radiotherapy range of 100-600 MeV/a.m.u, with a co-efficient of determination between model and tabulated data of 0.9968.

Model application.
A major motivation for the model outlined here is to rapidly simulate likely energy profiles for dose optimization. This model also readily allows for fast simulation of dose through different tissue types. A typical prostate radiotherapy planning scan is depicted in Fig. 6(a), with the planning target volume outlined in blue. This image slice shows at least three distinct regions, and the energy deposition was calculated with this model for hypothetical proton beams assuming varying tissue electron density as given in Table 2. The results of this are shown in Fig. 6(b) and (c) for proton and carbon ion plans respectively. The energy profile for both proton and carbon ions with variation of mean ionization potential of water is shown in Fig. 7 for a 250 MeV proton. The standard ICRU value of 75 eV was used for simulations in this work, but Paul 42 has suggested 80.8 eV better fits some experimental data. The net effect of an increasing estimate of I is an increase in observed range, and converse a decrease in range if lower estimates of I are employed. The physical explanation for this is relatively intuitive, as higher estimates for mean ionization potential translates to fewer ionization events and thus an increase in possible range.

Discussion
The model presented here is an accurate and analytical solution to the Bethe equation - Fig. 1 shows how the continuous solution compares to that derived from conventional Runge-Kutta methods. Results from both are virtually identical, suggesting that the model captures the equation dynamics well. A further advantage of the analytical solution is that it explicitly yields a particle range R T and consequently the Bethe equation-predicted energy, velocity and LET at any point throughout the track. When contrasted to MC simulations, Figs 2 and 3 show that the model captures the gross behaviour of a proton beam well but there are some important caveats to this. Primarily, the model slightly over-predicts particle energy along the track on average, with a mean error ranging from 0.81 MeV for 100 MeV protons up to 5.89 MeV for 250 MeV protons (Table 3), and discrepancies growing at an accelerating rate with particle energy. The most likely explanation for this is due to an inherent limitation of the Bethe equation, which does not intrinsically consider MCS or secondary (or higher order) production, which the MC simulation does. In particular, the magnitude of straggling has been observed to be approximately proportional to total particle range 4 -a trend also displayed in the mean errors reported in Table 3. The result would be additional observed energy loss at a given depth in MC, compared to the Bethe equation. This interpretation is supported by the inclusion of a "narrow" MC simulation also shown in Fig. 3 -particles in this simulation undergo less scattering and the results under this constraint lie closer to the continuous model solution. Model ranges were exceptionally close to tabulated PSTAR estimates, with errors of ≪1% for protons.
This has an interesting consequence around at the end of the particle track (typically defined in medical physics literature as the distal edge), at the particle's maximum range. As can be seen from Fig. 4, the Bethe equation predicts a asymptotic increase in LET at R T . The MC simulation, by contrast, agrees closely with model prediction right up until the Bragg peak, where higher order effects manifest. An important consequence of this is that the model (and indeed, the Bethe equation itself) does not capture the dynamics of particles at the distal edge completely, instead only approximating the general behaviour around this point. For carbon ions, the predicted ranges and experimentally measured Bragged Peaks agreed closely (mean percentage errors 0.21%-6.14%) and again the bulk of this discrepancy might be explained by fundamental limitations of the Bethe equation around the distal edge. Model agreement with MSTAR data was high (≈1% error) even to energies beyond those employed in therapy (600 MeV/a.m.u).
It is crucial to note that the solution outlined here is not intended as a replacement to conventional MC methods, but as a useful means of rapid optimization. As the analytical form presented here can be rapidly implemented with minimal computational costs, it could readily be applied in optimization when multiple iterations of dose calculation are needed. Using the model here, a huge array or potential doses can be quickly simulated, and narrowed down to only the most promising, which then can be more robustly simulated with MC techniques as a final calculation, potentially saving considerable time and expense. The other major advantage of the expression presented here is that it allows rapid simulation of complicated cases involving heterogeneous tissue, based on well-defined material constants such as electron density and mean ionization potential instead of water-equivalent stopping power ratios (SPRs). An example is shown in Fig. 6. This is readily implemented and,  as it is analytical, the approximate dynamics of particles in even complex media can be readily estimated. There is still debate over the ideal value for mean ionization potential for water, and the model was also applied to investigate this. Figure 7 suggests that this difference chiefly manifests at the distal edge.
The method outlined here produces a solution for the Bethe equation, and hence the velocity, energy and LET curves for any charged particle rapidly from first principles. This is very useful, but it could potentially be extended to calculate Bragg curves in 1D by modelling energy variation and range straggling, or even in 3D by detailed consideration of MCS. These effects are the likely explanation for the departures between MC and model data seen at the distal edges in Figs 3 and 4. A number of semi-empirical formulations exist for multiple Coulomb scatter 43 and energy variation 4 , but a detailed comparison of such approaches beyond the scope of this work. Examples of such Bragg curves are illustrated in Fig. 8, which were created using energy profiles from this model, a Gaussian approximation to the Vavilov energy straggling distribution 44 , and step-wise calculation of multiple Coulomb scatter using the Rossi formula 45,46 . Additionally, the Bethe equation considers only stopping due to electrons within a material, but a full consideration of dose deposition mechanisms in tissue would also include nuclear effects.
For majority of the track length, the Monte Carlo LET agrees closely with the model as illustrated in Fig. 3. However, towards the distal edge the model (and Bethe equation itself) displays asymptotic behaviour. We can naively calculate the maximum obtainable LET from the Bethe equation presented in eq. 6 -defining this as S, we can readily calculate dS dv . At some velocity v m the particle is at its minimum speed and thus S tends to a maximum. Setting This analysis yields v m = 0.0141c for charged particles. The maximum possible LET for protons from the Bethe equation is thus 84.35 KeV/μm and 3036.72 KeV/μm for carbon in water. In reality, particles will not reach anything near this, due to deflections and energy losses from MCS, nuclear interactions and secondary production 4 . These are factors not intrinsically considered by Bethe equation or this model, and this in part explains the disparity between the model and Monte Carlo results at the distal edge as illustrated in Fig. 3. Tellingly, LETs are higher in the 'narrow' phantom Monte Carlo, reflecting that fact that maximum LET increases as the potential for scatter decreases.
It is important to note too that the Bethe equation is a continuous function, whereas charged particle interactions are in reality discrete events. Water and other discrete biological targets also have finite extent for energy transfer 47 , and we do not expect LET to reach these maximums. The resolution of the simulations is also worth considering; the model was evaluated continuously and plotted with a fine 1 micron resolution, whereas the MC has a much coarser scale of 1 mm. If a similar length scale is used for the model, the asymptotic behaviour at the distal edge disappears. The model presented is certainly not a substitute for MC calculations of LET, but captures the general behaviour well, though caution should be taken in interpreting results near the distal edge.
The ease and speed with which this formulation of the Bethe equation can be applied is a particular strength, as illustrated by the profiles shown in Fig. 6 for tissue of varying electromagnetic interaction properties. LET is also important in determining oxygen contribution to radiation damage, as molecular oxygen modulates the effectiveness of standard photon radiotherapy by up to a factor of three 48,49 . It is known that oxygen enhancement varies markedly with LET 50 , and in principle accurate determination of LET for charged particles could be applied in conjunction with information on regions of chronic and acute hypoxia 51 , potentially improving treatment further. While this aspect remains speculative, the model presented in this work should lend itself to rapid approximation of particle range, energy and LET even in complex tissue, and ultimately should be of benefit in optimizing treatment and outcome.

Conclusion
The solution presented here to the Bethe equation can be rapidly applied to any charged particle type, and yields accurate estimates of energy profile, velocity and LET through a medium. These are shown to agree well with experimental and Monte Carlo results for both protons and carbon ions, with some discrepancy at the distal edge. The approach outlined also lends itself to the rapid calculation of dose profile through different media, and could serve as a basis for rapid dose optimization for charged particle therapy.