Development and implementation of a non Gaussian model for the lateral dose prediction in a proton therapy treatment planning system

Challenging issues in treatment planning system for hadrontherapy are the accurate and fast calculation of dose distribution, the reduction in memory space required to store the dose kernel of individual pencil beams and the shortening of computation time for dose optimization and calculation. In this framework, the prediction of lateral dose distributions is a topic of great interest because currently the double gaussian parametrization is typically used as approximation although other parameterizations are also available. The best accuracy for this kind of calculations can be obtained by Monte Carlo methods, at the expense of a long computing time. This work aims to present a flexible computational model for the calculation of the lateral profile of a pencil proton beam and the results of its implementation in a treatment planning system. The model calculation are compared with the currently used double gaussian approximation and the Monte Carlo calculations, and the tests are performed in water and in presence of inhomogeneities.


Introduction
Particle radiotherapy treatments constitute about the 1% of the total number of treatments for patients receiving radiotherapy worldwide [1]. Why particle therapy? The answer is based on two main characteristic of ions and protons [2]. The physical advantage respect to conventional radiotherapy, is due to the characteristic depthdose profile, analytically based on the Bethe-Bloch formula described by the Bragg curve, characterized by a small entrance dose and a distinct maximum near the end of range with a sharp fall-off at the distal edge, called Bragg peak, that results in a customizable and precise dose deposition. From a radio-biological point of view protons and heavier ions have different Radio-Biological Effectiveness (RBE) [3]; protons and photons have a similar RBE in the plateau region, and higher RBE on the peak region while heavier ions show an even higher RBE the Bragg peak, compared to photons. This is due to the high ionization density produced by the ions as they traverse a cell, resulting in more disruptive damage to the DNA double helix. The higher RBE combined with the lower Oxygen-Enhancement Ratio (OER) [3], especially in and around the Bragg peak, can be a considerable advantage in tumour cell killing, particularly for resistant tumour cells [4]. Analytically, the Bragg curve is the results of three interactions processes: stopping, scattering and nuclear interactions. Focusing on the lateral beam shape (on the plane perpendicular to the beam axis), these effects result in a broader shape, compared with the one expected from a purely electromagnetic calculations. The angular deflection is mainly due to the Multiple Coulomb Scattering (MCS). In a first approximation, the MCS angular distribution is nearly gaussian, and for this reason, is approximated with a double Gaussian function [5,6].
To predict the exact form of the MCS angular distribution and its characteristic width as a function of proton energy as well as scattering material and thickness, we propose a flexible computational model based on Molière's theory [7,8], that requires as input the particle energy (or the beam energy distribution as a phase-space) and the desired material depths.

The model
The model is composed by a purely analytical electromagnetic core based on Molière's theory and a nuclear tails part. It has been developed for a proton pencil beam in water but it can be easily extended to more realistic cases to reproduce mathematically the lateral beam scan in a given region [9].

Electromagnetic core
The Molière theory is mathematically based on the standard transport equation, whit the Bessel transforms and the small angle approximation i.e. sin(θ) ≈ θ. Physically, the theory only accounts for electromagnetic interactions and assumes the Rutherford form for the single scattering cross section: where ρ is the number of atoms per cm 2 , and the integration over φ is performed. χ c and χ α are the two crucial parameters of the model. The first one is connected to the rms scattering angle: where z is the atomic number of the incident particle, Z and A the atomic and the mass number of the target, x is the thickness of the target (g/cm 2 ), p is the momentum in GeV/c and β = v/c. The second one is the parameter connected to the electron screening of the Coulomb potential: When µ = 0 (and consequentially χ 2 α = 0) Eq. (1) is the Rutherford cross section of a point charge. The standard form for the Molière distribution is given by [10] and the projected angle θ x (that has the same trend of θ y ) follows the distribution described in [10]. In the small angle approximation, using simple geometry consideration, the connection between the mean square of the projected angle in the transverse plane and the spatial angle is obtained: Therefore, the connection between the rms of the space and projected angles is: 2.1.1. Mixture and compound The generalization of this theory to a compound or mixture is quite easy, by the modification of Eq. (2), (3) and (5) [9] and considering a proper summation over the constituents as where Z i and A i are respectively the atomic number and atomic weights of the molecule or mixture of n i constituents.

Energy loss and lateral displacement
Since in realistic hadrontherapy cases the paths of protons in water is very long, the basic equations have to be modified to take into account the energy loss process. Considering that the parameters χ c and χ α depend on p, β, the energy loss problem can be solved if one finds the dependence of these quantities on the water thickness x traversed [9]. Including the energy loss effects, the general formula for passing from angular to spatial displacements can be derived. The rms y M of the transverse displacement on a measuring plane at x due to a layer dt at the depth t is given by (x − t)θ xR , where the angle is given by Eq. (7), that presents a logarithmic dependence on χ c , which contains the thickness x (see Eq.(2), (8)). The physical meaning is that two successive layers act in a dependent manner, since the second layer receives trajectories deflected by the first one. The following final formula can be obtained (details in [9]):

Nuclear effects
To complete the description of the overall phenomena, nuclear effects have to be properly taken into account. An additional term is added to describe the interaction of non-primary particles that affects the tails to the distribution. This is particularly relevant when big thicknesses are involved, such as in the case of hadrontherapy (it can be estimated as about 1% per cm of depth).
To account for these effects, a modified Cauchy-Lorentz distribution [11,12] is considered: where the three free parameters are the amplitude A, the half width half maximum (HWHM) b and the variance σ 2 .

Complete model and validation
Now, the two parts are combined by adding a weight factor W p [13], that express the percentage of events without nuclear interactions (primary protons), as a function of the traversed thickness [9]. The total normalized final distribution for the lateral displacement, can be formulated as: where f M (y) is the distribution of Eq. (10) and W p is the weight factor. Both f M (y) and t(y) are normalized to unit area. The free parameters A, b and σ 2 , have been determined by fitting the complete FLUKA MC lateral displacement distributions, and they are parametrised with Chebyshev polynomials. These results have been validated against FLUKA [14,15] and HIT Heidelberg Ion Beam Therapy Center [16]. A selection of results is reported in Fig.1. For all the analysed depths and energies, the model agrees very well with the FLUKA predictions and with the measurement data.
Energy (a.u.) 2 − In the normalization zone the error of the data points is 2%.

Treatment Planning System (TPS) implementation
To test the possible clinical advantages of using our model, these results have been implemented in CERR -A Computational Environment for Radiotherapy Research TPS [17,18] at Ludwig-Maximilians-Universität München. The rational was to compare the dose deposition obtained using our model with the results obtained using the implemented double Gaussian approximation (DGA), taking as gold standard the results of the MC method. Several energies and depths have been considered for single beams and for a full treatment plan in water. Furthermore, plans containing inhomogeneities i.e. bone and fat, have been studied. For a more detailed explanation see [19]. The doses evaluated by the TPS, have been analysed quantitatively using statistical tests, lateral profiles comparison, residual analysis and DVH (Dose Volume Histogram) for the full plans. A selection of results is reported in Fig.2 for single beam in water and in presence of bone and fat inhomogeneities and for a full plan in Fig.3 (Right). Fig.3 (Left) shows the direct comparison between the model DVH and the DG ones for the full plan; the delivered dose to the PTV (Planning Target Volume [3]) is almost equal in the two cases but the dose on surrounding healty tissues ((PTV3D+1)-PTV) is lower with the model. Overall one can conclude that a plan obtained using the model for the lateral evaluation of the dose, is a promising candidate instead the DG ones in terms of treatment quality.

Conclusion
This study shows that the model has better accuracy than DGA method with adequate computation time in water, and it is in a good agreement with the MC for all the analysed  cases. Since the comparison with data and MC shows very good agreement, the model can be proposed to be used instead the currently used DGA. This solution can enhance the accuracy of the lateral profile evaluation in hadrontherapy practice, that can be a critical parameter especially in the case of deep tumours near organs at risk. Since the model is much faster that the MC [20,21] and has same accuracy in water, it can be proposed in future to be used for the online beam monitoring, and also as a powerful and flexible tool for the forward dose evaluation. A future step will be the study of the effects of different materials and different geometries, including interfaces and the use of other particles already clinically used, e.g. carbon ions, α particles.