Local weighting of nanometric track structure properties in macroscopic voxel geometries for particle beam treatment planning

The research project BioQuaRT within the European Metrology Research Programme aimed at correlating ion track structure characteristics with the biological effects of radiation and developed measurement and simulation techniques for determining ion track structure on different length scales from about 2 nm to about 10 μm. Within this framework, we investigated methods to translate track-structure quantities derived on a nanometre scale to macroscopic dimensions. Here we make use of parameterizations that link the energy of the projectile to the ionization pattern of the track using nanodosimetric ionization cluster size distributions. They were defined with data generated by simulations of ion tracks in liquid water using the Geant4 Monte Carlo toolkit with the Geant4-DNA processes. For the clinical situation with a mixed radiation field, where particles of various energies hit a cell from several directions, we have to find macroscopic relevant mean values. They can be determined by appropriate local weighting functions for the identified parameterization. We show that a stopping power weighted mean value of the mentioned track structure properties can describe the overall track structure in a cell exposed to a mixed radiation field. The parameterization, together with the presented stopping power weighting approach, show how nanometric track structure properties could be integrated into treatment planning systems without the need to perform time consuming simulations on the nanometer level for each individual patient.

without the need to perform time consuming simulations on the nanometer level for each individual patient.
Keywords: Geant4-DNA, nanodosimetry, proton therapy (Some figures may appear in colour only in the online journal)

Introduction
An increasing percentage of cancer patients is treated with high energy proton or ion beams. Compared to conventional photon therapy, the tumor can be targeted with higher precision and a dose escalation can be reached in the tumor volume without having to accept a higher dose in organs at risk (Schulz-Ertner and Tsujii 2007). However the enhanced biological effectiveness of ion beams and its underlying physical processes are not yet fully understood. The BioQuaRT 4 project (Rabus et al 2014, Palmans et al 2015 aimed at creating a new measurable dosimetric quantity to define the local radiation quality based on track structure characteristics. This quantity will have a strong influence on the biological response and a future goal is to integrate this quantity into clinical treatment planning systems which simulate and optimise the treatment for each individual patient. To integrate a quantity defined on the nanoand micrometer scale into clinical treatment planning, one could in principle perform Monte Carlo simulations on the nanometer level individually for each patient and beam configuration (i.e. treatment plans). However, this would be extremely time consuming, and we propose an approach to derive nanodosimetric quantities in patient geometries more efficiently. First, this requires a parameterization (depending on particle type and energy) which enables a fast calculation of the relevant track structure properties for individual particles. Then, local weighting functions for the track structure parameters can be defined, which reflect the situation when cells are exposed to a mixed radiation field. This mixed radiation field is composed of particles of various energies and particularly in the case of ion beam therapy-of nuclear fragments as well. We present a stopping power weighting approach to calculate the relevant characteristics of the overall track structure in cells very efficiently, if the energy spectra of all ion species at the point of interest are known. This leads to macroscopic values in millimetresized voxel geometries, which are essential in the clinical situation for treatment planning. With the energy-dependent parameterization of track structure properties and the presented weighting method nanometric track structure characteristics can be projected quickly onto the macroscopic scale. Together with a biological model added on top this method could then be integrated into treatment planning systems to enable an optimisation of treatment plans based on nanodosimetric description of local radiation quality.

Data simulation with Geant4-DNA
Data for the analysis have been produced using the Geant4 Monte Carlo toolkit (Agostinelli et al 2003) version 9.6 and, in particular, the Geant4-DNA processes (Incerti et al 2010a(Incerti et al , 2010b that are included in the low energy package. The simulation with the Geant4-DNA processes is performed on step-by-step basis (i.e. avoiding the condensed history approximation). Due to the low energy cut for the electron transport (7.4 eV in this work) it is possible to simulate the energy depositions with nanometric resolution, i.e. to perform a detailed track structure simulation. For more details, the reader is referred to the Geant4-DNA documentation 5 . We simulated tracks of protons within a homogeneous liquid water cube of 20 μm side length. The primary particles were started in the centre of the water cube at point (0, 0, 0). Traveling in the z-direction this resulted in tracks with a length of 10 μm. In a second step we cut a track of 1 μm length out of the full track starting 1 μm behind the point where the simulation was started. We found that this procedure ensures that the considered track segment also contains energy depositions by secondary electrons created further upstream. The position and the amount of deposited energy was recorded for each individual ionization of a water molecule. The simulations were performed for protons in a clinical relevant energy range between 1 MeV and 100 MeV (1 MeV, 3 MeV, 5 MeV, 10 MeV, 25 MeV, 50 MeV, 100 MeV). As the energy loss of the primary particles along the first 2 μm of the track is lower than 6% for the lowest proton energy simulated, the projectile energy can be considered as constant. The number of simulated tracks was between 100 and 500, leading to at least 13 000 ionization events for high proton energies such as 50 MeV and 100 MeV. The number of ionization events for lower proton energies was even higher.

Calculation of ionization cluster size distributions
It is generally accepted that damage to the DNA in the cell nucleus is the primary reason for the occurrence of biological effects after irradiation (Goodhead 1994, Olive 1998. The field of nanodosimetry aims at describing the stochastic nature of particle track structure on the nanometer scale in order to explain radiobiological phenomena with the help of the concept of local radiation quality (Grosswendt 2004, 2005, Grosswendt 2006, 2010, Rabus and Nettelbeck 2011, Bug et al 2014. Therefore we calculated ionization cluster size distributions (ICSDs) which describe the probabilities for ionization clusters of a particular size to occur in nanometric volumes. The ionization cluster size is defined as the number of ionizations produced by a particular passing particle track in a considered target volume. The size of these volumes was chosen to be cylindrical with dimensions comparable to a DNA segment of ten base pairs (one convolution of the DNA double helix). The higher the probability for large cluster sizes-which may produce irreparable double strand breaks-the higher is the probability for the occurrence of radiobiological effects. Statistical parameters of ionization cluster size distributions (de Nardo et al 2002) are promising regarding their potential to link the ionization pattern of the particle tracks to the radiobiological outcome.
In order to obtain data for parameterizing the track structure properties, we calculated ICSDs for all simulated proton tracks for all considered energies. The cluster size was determined by scoring the number of ionizations in cylindrical target volumes. As usual in nanodosimetry we chose cylinders which resemble the geometry of a DNA segment of 10 base pairs with a height of 3.4 nm and a diameter of 2.3 nm. For each track, a cuboid was placed around the primary particle trajectory, the size of which was determined by the outermost ionizations and a margin of 10 nm that was added to avoid effects at the edge of this volume. Then 10 million cylinders were placed in the cuboid at randomly distributed locations. The orientation of the cylinders was chosen randomly as well, and cylinders could overlap. Conditional ICSDs were then determined by counting the ionizations in each cylinder where at least one ionization occurred. The distribution was then normalized so that the probability to obtain at least one ionization became a hundred percent. With the resulting relative frequency distribution P Q ( ) ν| , where Q is the radiation quality and ν the ionization cluster size, the conditional 5 www.geant4-dna.org.
probability to obtain a particular cluster size can be derived. Figure 1 shows the resulting relative frequency of cluster sizes between 1 and 10. For 1 MeV protons the relative frequencies for large cluster sizes are significantly higher and decrease with increasing energy of the incident proton.
As mentioned above statistical moments of the ICSDs may have the potential to explain radiobiological effects. In this analysis we make use of the conditional mean ionization cluster size M 1 (Q) which is calculated by summing up the products of the cluster size ν and the corresponding conditional probability P: Other moments of the ICSD used in this analysis are the complementary cumulative probabilities to obtain a cluster size larger than two or three, respectively: These moments of the distribution are potentially correlated to the probability to obtain simple or complex double strand breaks of the DNA molecule (Grosswendt 2005, 2010 and hence candidates as input quantities for a treatment planning system based on measurable track structure properties.

Energy-dependent track structure parameterizations
In order to characterize the ICSDs we calculated the statistical moments presented in (1)-(3) of each distribution. These quantities can then be plotted against the energy of the incident projectile to obtain an energy-dependent parameterization (Alexander et al 2015). Figure 2 shows the result for M 1 calculated for all simulated projectile energies E and then plotted as a function of E. The data set was then fitted with a (modified) power law consisting of two terms: The best fit parameters are a 0.67 0.04 1 = ± , b 0.58 0.08 1 = − ± and c 1.32 0.04 1 = ± with 95% confidence intervals. The fit curve is shown as blue line in figure 2. Figure 3 shows F 2 versus the projectile energy. The relationship between the energy E and F 2 can be described by a power law as well: with a 0.28 0.03 2 = ± , b 0.33 0.07 2 = − ± and c 0.17 0.03 2 = ± . Similarly the integrated probability to obtain a cluster size larger than three (figure 4), can be described with the following power law: with a 0.19 0.01 3 = ± , b 0.53 0.05 3 = − ± and c 0.06 0.01 3 = ± for protons. With these parameterizations the statistical moments of the ionization cluster size distribution can be calculated very efficiently for any desired proton energy (within the shown range), without having to perform Monte Carlo track structure simulations during treatment planning.

Weighting functions for a nucleus in a mixed radiation field
In the clinical situation, a cell is typically exposed to a mixed radiation field consisting of primary particles of different energy, high energetic electrons and nuclear fragments. Hence a single cell is hit by different particle species of different energy (i.e. multiple tracks). Therefore The blue line shows the best fit curve for a power function with two terms: the aim in this work is to find local weighting functions for the track structure properties of all tracks hitting this cell in order to enable an efficient calculation of relevant track structure properties in a mixed radiation field. Here we concentrate on protons only, but the approach could also be extended to other particle species.

Nucleus in a mixed radiation field
First, we will analyze the overall track structure properties in a nucleus which is traversed by two different particles. We assume that the nucleus is the primary target and the DNA molecule is homogeneously distributed in the nucleus. The sensitive volume (e.g. the nucleus or a part of it) is modeled by a spherical volume with a diameter of 1 μm (green dots in figure 5). The direction of the tracks traversing this nucleus was chosen randomly. Hence the length of the tracks in the nuclei is different for each track. Blue and red dots in figure 5 show the  Figure 6 shows the individual ICSDs of two single proton tracks with energies of 1 MeV (blue dashed line) and 10 MeV (red dashed line). The relative frequency to obtain a large cluster of ionizations is larger for the 1 MeV proton compared to the 10 MeV proton. In the next step we analyzed the overall track structure in the nucleus, with both tracks present. By scoring the number of ionizations in 100 million randomly distributed cylinders, the situation of a mixed field was simulated. The result shows that the overall ICSD is dominated by the ICSD of the 1 MeV proton. In order to calculate the overall ICSD from  the ICSDs of the single tracks, the trivial approach would be a simple track averaging (i.e. fluence weighting) of the ICSDs of the single tracks. The overall ICSD for a set of tracks t t t , ,... N 1 2 is then given by:

Dose-Weighting in a single nucleus
This trivial approach (orange circles in figure 6) underestimates the influence of the 1 MeV proton. A more sophisticated approach is a dose averaging (i.e. specific energy weighting) according to where d i is the dose (or specific energy) contribution of track i in the investigated volume and D d i i = ∑ . Figure 6 shows that the dose averaged ICSD (turquoise triangles in figure 6) agrees very well with the ICSD resulting from the overall track structure (green triangles). Hence the overall ICSD in a single nucleus can be calculated with the knowledge of the energy of the tracks traversing this nucleus and the dose deposited by these tracks.

Comparison of stopping power weighted values with Monte Carlo results
As shown in section 4.2 the ICSD for a nucleus hit by two particles can be calculated by performing a specific energy weighting of the ICSDs of the individual particles. In the clinical situation our knowledge is restricted to the type of particles hitting the voxel and their energy spectra. Therefore we have to show that a stopping power weighting reproduces averaged values which are in agreement with the results of the overall track structure in the nuclei which were hit. The difference to the situation in section 4.2 is that-in the clinical situation-we do not know the exact, individual ICSDs of the single particles hitting the nuclei (nor their length respectively dose contribution). Hence we will use a modified version of the specific energy weighting, where the contributions of the individual tracks are weighted by the respective stopping power (since the average specific energy will be the stopping power times the average track length in the nucleus).
In a future planning system, the calculations can be performed with the statistical moments of the ionization cluster size distribution as the parameterizations presented in section 3 enable a fast calculation. Thus the following analysis will be based on the statistical moments already parameterized in section 3. Figure 7 shows the result for different combinations of proton energies. The data points connected with solid lines show the results of the analysis of the Monte Carlo generated data set. Each data set corresponding to a data point in the plot consists of one hundred nuclei containing two randomly orientated tracks each. The data points marked with asterisks in figure 7 correspond to the 1 MeV series, i.e. the average of hundred nuclei that are hit by a 1 MeV proton and by a proton of the energy indicated on the x-axis. Data points marked with squares corresponds to the 5 MeV series and circles to 25 MeV. The overall track structure was determined in each nucleus, a mean ICSD was determined for all nuclei and the statistical moments presented in section 2.2 were calculated. The results of these simulations were then compared to the track structure properties calculated with the modified stopping power weighting method: where S is the sum of stopping powers of all tracks, S i is the stopping power of track i and P(t i ) is the track structure property of track i which was obtained from the parameterizations of the original simulation data presented in section 3.

Discussion and conclusion
As it is not possible to perform full Monte Carlo track structure simulations with nanometer resolution for all particle species and energies which are present in the clinical situation, a fast assessment of relevant track structure parameters is crucial. For this purpose we calculated ICSDs for protons of various energies and calculated the according statistical moments of the ICSDs in order to characterize them. The presented parameterizations (4)-(6) in section 3 enable a fast assessment of the statistical moments of the ICSDs of protons in a clinically relevant energy range between 1 MeV and 100 MeV. The upper limit of 100 MeV was given by the current capabilities of Geant4-DNA. For higher energies, variations in radiation quality are small and the data in figures 2-4 could easily be extrapolated to higher proton energies. The preeminent advantage of this nanometric characterization of particle track structure is that ICSDs and the associated statistical parameters are measurable quantities (Garty et al 2002a, 2002b, Bantsar et al 2004, Conte et al 2012. Hence, measurements of these quantities can be used to characterize track structure and radiation quality of the treatment beam. Yet, further effort is necessary to develop portable instruments that can be used in the clinical routine.
When it comes to treatment planning, we have to deal with a mixed radiation field where each cell is hit by particles of different energies from several directions. A full Monte Carlo track structure simulation for each possible treatment plan involved in the optimization process would be too time consuming. Therefore we presented a stopping power weighting approach in section 4.3 that provides an efficient calculation of relevant statistical properties of the overall ICSD in the target (i.e. the homogeneously distributed DNA molecule in the nucleus). When comparing the calculated results for the mean ionization cluster size M 1 shown in figure 7(top) the mean deviation between the calculated and simulated result for M 1 is 3.4%. For F 2 (center) it is 4.7% and for F 3 (bottom) 8.9%. The highest deviations for F 2 and F 3 both occur for data points where 25 MeV proton tracks were mixed with proton tracks of another energy. These deviations are therefore most likely due to the fact that 25 MeV protons ionize very sparsely compared to 1 MeV or 5 MeV protons. Hence the individuality of the tracks influences the result and causes these deviations. Nevertheless in general the calculated result agrees very well with the simulation result. In conclusion, the presented stopping power weighting approach can be used as a tool to obtain statistical moments of the overall track structure in a mixed radiation field. The required input are the energy spectra of the particles hitting the target and the track structure properties of the single tracks, that can be extracted from the parameterizations presented above. The weighting is then performed by using the respective stopping power. In combination with an appropriate biological model added on top to predict the radiobiological outcome, the presented method can be useful to achieve a higher accuracy for biological optimization of treatment plans. One example for such a model is the RMF model which links double-strand break induction to reproductive cell death (Carlson et al 2008, Frese et al 2012.
Overall, the presented parameterizations and the weighting method can quickly provide the crucial features of the overall track structure in a mixed radiation field without having to perform Monte Carlo simulations with nanometer resolution. As an efficient way to calculate these features, they can be the basis for future treatment planning systems with biological optimization based on nanometric track structure properties. Incerti