Comprehensive track-structure based evaluation of DNA damage by light ions from radiotherapy-relevant energies down to stopping

Track structures and resulting DNA damage in human cells have been simulated for hydrogen, helium, carbon, nitrogen, oxygen and neon ions with 0.25–256 MeV/u energy. The needed ion interaction cross sections have been scaled from those of hydrogen; Barkas scaling formula has been refined, extending its applicability down to about 10 keV/u, and validated against established stopping power data. Linear energy transfer (LET) has been scored from energy deposits in a cell nucleus; for very low-energy ions, it has been defined locally within thin slabs. The simulations show that protons and helium ions induce more DNA damage than heavier ions do at the same LET. With increasing LET, less DNA strand breaks are formed per unit dose, but due to their clustering the yields of double-strand breaks (DSB) increase, up to saturation around 300 keV/μm. Also individual DSB tend to cluster; DSB clusters peak around 500 keV/μm, while DSB multiplicities per cluster steadily increase with LET. Remarkably similar to patterns known from cell survival studies, LET-dependencies with pronounced maxima around 100–200 keV/μm occur on nanometre scale for sites that contain one or more DSB, and on micrometre scale for megabasepair-sized DNA fragments.

In the determination of cross sections for particles with atomic number Z > 2, the Barkas scaling term 1 for the effective charge has been replaced by its relation to the Barkas term for Z=1: σZ = σH Z²(1 -exp(-125(v/c) Z -2/3 ))² / (1 -exp(-125v/c))². (S1) The cross section for hydrogen σ H can be separated into a fraction due to proton interactions, σ 1 , and a fraction due to neutral hydrogen interactions, σ 0 : where f + is the energy-dependent fraction of protons during the slowing-down process.
In a first approximation the energy-dependent fraction of proton states during slowing-down was estimated by the effective charge term for Z=1 according to the Barkas
This estimate resulted in energy dependent stopping powers and ranges for carbon which have been found in reasonable agreement with ICRU data 2 , SRIM calculations 3 and recent measurements 4 . However, a comparison of this approach with a usual Monte Carlo slowing-down calculation for hydrogen has shown deviations in the energy-dependent charged fractions, indicating that the underlying approach was not fully self-consistent. Indeed, the fractions of charged states f + and neutral states f 0 can be derived 5,6 due to their dynamic equilibrium from the cross sections for electron stripping of neutral hydrogen σ 01 and electron pick-up by protons σ 10 f+ = σ01 /(σ01 + σ10), (S4) f0 = 1 -f+ = σ10 / (σ01 + σ10).
For application in the track structure code, however, the primary quantity is the mean free path λ Z , i.e., the inverse of the cross section multiplied by the given target density N. The scaling law for λ Z is λZ = (1 -exp(-125v/c))² λH / Zeff².
Thus, the change from cross sections to their inverse quantities, the mean free paths, entails a transition from fractions of charge states, Eqs. (S4-S5), to fractions of charge-changing processes, Eqs. (S8-S9). The algorithm according to Eqs. (S6-S9) has been implemented in the PARTRAC code for the transport calculation of particles with Z > 2; differential cross sections are taken from the proton or neutral hydrogen data sets with probabilities according to Eqs. (S8-S9).
In Figure S1, various data on fractions of charged states are presented for hydrogen particles as a function of energy, including the fractions f + (Eq. (S4)), p + (Eq. (S8)) and the earlier Barkas approximation Z eff,H (Eq. (S3)). At low energies the fraction of charged states is rather low (e.g. less than 10% below 5 keV). However, even at such low energies the fraction of stripping processes among interactions of neutral hydrogen states is about 25%. Therefore, on average after 3 inelastic interactions of neutral hydrogen the electron is stripped, and subsequently a proton interaction will happen; with high probability (below 5 keV above 80%) this will be an electron pick-up process yielding immediate re-formation of the neutral state after just a short transport distance. Thus, the fraction p + of proton interactions regarding the sequence (or probability) of processes is about 20% at the low energy limit.
Characteristic for charge changing processes during hydrogen slowing-down is that the number of electron pick-up processes is always equal to or exceeds by one the number of electron stripping processes. In order to safeguard this balance between both processes in a Monte Carlo procedure, in this work they have been combined into a single process including an ionization of the liquid water medium (due to pick-up) and an emission of a secondary electron with the speed of the energetic particle in forward direction (due to stripping). The combined charge-changing processes have been associated to interactions of neutral hydrogen states. However, pick-up processes by protons must not be excluded in the particle transport. To obtain the same transport behaviour in the approximate track structure simulation, and consequently the correct stopping power, electron pickup processes have been included as transport events without any interaction.

Validation of cross-section scaling
To test the validity of the generalized Barkas procedure, Eq. (1), it has been applied directly to stopping powers. In Fig. S2a-d, the electronic stopping powers for He, Li, C, and Ne ions (solid red lines) obtained by scaling the ICRU data for hydrogen 7 are compared to the standard values (black points) recommended by the ICRU 2 . In Fig. S2, solid green lines depict stopping powers obtained by an analogous scaling method that starts from ICRU data for helium 7 instead of hydrogen; this procedure exactly reproduces the He data ( Fig. S2a), verifying its self-consistence. The hydrogenbased scaling underestimates the maximal stopping power for helium ( Fig. S2a) by about 4% and shifts its location to about 40% lower energies compared to the ICRU database (217 keV/µm at 0.1 MeV/u vs. 226 keV/µm at 0.18 MeV/u); this reflects the very different atomic structure of helium compared to that of hydrogen. Note, however, that the scaled results for helium are shown here for comparison only; they are not used in PARTRAC, since in the code there is a dedicated cross-section database for helium ions 8 that explicitly accounts for its charge states. For light ions with Z > 2 (Figs. S2b-d), the stopping powers from hydrogen-based generalized Barkas scaling (solid red lines) are in good agreement with the ICRU data 2 , but underestimate the maximum stopping powers by a few percent. The generalized scaling works much better than the original one (dashed red lines), especially at energies below ~0.5 MeV/u.
In the literature, also a three-term expression can be found for the effective charge 9 , As illustrated for carbon ions (Fig. S2c, dash-dotted red line), this alternative effective charge formula, Eq. (S10), produces results very close to those from Eq. (2), used throughout this work.
The refined Barkas scaling works surprisingly well for ions as heavy as neon (Fig. S2d). For all studied light ions with Z > 2, scaling based on hydrogen is superior to the one that starts from helium data (solid and dashed green lines for the refined and original procedures, respectively).
The results obtained with the scaling procedure cover only electronic stopping powers. PARTRAC does not include nuclear interactions, and hence it does not assess nuclear stopping powers or fragmentation of the projectile ion in nuclear reactions. ICRU data are available on nuclear stopping powers for helium 7 (Fig. S2a, blue points), and SRIM calculations 3 do reproduce these values and provide estimates of nuclear stopping powers for all the studied ions (thin blue lines in Fig. S2). Above 100 keV/u, nuclear stopping powers are negligible, and hence they have not been plotted here. At energies below about 10 keV/u, however, nuclear processes do contribute significantly to the ions' stopping. At very low energies, nuclear stopping powers exceed the electronic ones; according to SRIM calculations, this occurs below 0.17, 0.69, 1.5, 1.6, 2.0, and 3.0 keV/u for H, He, C, N, O, and Ne ions; the corresponding values of nuclear and equally electronic stopping powers are 6, 18, 76, 91, 104 and 123 keV/µm, respectively. At about twice higher electronic stopping powers, the nuclear ones amount to about a fourth of electronic ones, i.e. nuclear processes make about 20% of the total stopping only. These values can serve as estimates of the low-energy limit of the present PARTRAC simulations.
To further validate the effective charge concept, it has been adopted in Monte Carlo simulations of the full slowing-down of hydrogen. In Fig. S3a, it is shown that this method provides results (solid red line) that are fully in agreement with usual simulations (grey dotted line) that explicitly track the charge states, i.e. distinguish between neutral hydrogen atoms and singly-charged hydrogen ions (protons) and use dedicated cross sections for these two states. PARTRAC results generally agree with the ICRU data 7 (blue squares). Around the maximal stopping power, the results of PARTRAC are by up to 6% higher than the ICRU values, and closely follow the calculations by the SRIM 3 code (black crosses). To provide a more detailed picture of the processes, contributions to the stopping power have been separately scored for ionisations by protons (green line), excitations by protons (orange line), ionisations by neutral hydrogen (yellow line), and by charge changing processes (violet line). , the ICRU database contains data on nuclear and total stopping; these are shown by blue and magenta points. Also shown are electronic, nuclear and total stopping powers from SRIM 3 (black, blue and magenta lines, respectively). Since nuclear stopping is negligible above 0.1 MeV/u, it has not been included in the plots in that region. In panel c, also shown is the scaling based on an alternative formula for the effective charge, Eq. (S10).
A track-structure simulation of full stopping within the effective charge concept has been performed also for carbon ions, i.e. with interaction cross sections obtained from hydrogen data by scaling according to Eq. (1). 100 tracks of 10 MeV/u carbon ions have been simulated, and the primary particle interactions have been scored. The stopping powers from this Monte Carlo procedure (Fig. S3b, red line) are practically identical with those from directly scaling hydrogen stopping powers (Fig. S2c, solid red line). They closely reproduce the ICRU data 2 (Fig. S3b, blue squares) but underestimate the stopping power around its maximum by almost 7%; in comparison, the SRIM 3 code (Fig. S3b, black crosses) predicts the maximal stopping power to be about 10% higher and shifted to about 30% higher energies than the ICRU. In the low-energy range, the present approach perfectly agrees with the ICRU data; for energies < 25 keV/u, the low-energy limit of the ICRU dataset, the present results show trends similar to those predicted by the SRIM code.

Comparison of DNA damage predictions in lymphocytes and fibroblast nuclei
Additional simulations have been performed for modelled nuclei of human fibroblasts irradiated by 1 H, 4 He, 12   Taken together, the setup for fibroblast nuclei is in several respects an alternative approach to the settings in the calculation for lymphocyte nuclei. Due to the flat shape of the nucleus and shorter transport distance, calculations with lower energies could have been performed. However, the DNA structure on scales up to 50 nm, the track structures, and the method of assessing DNA damage by overlapping the DNA model with the tracks are the same in both approaches.
To point out the effectiveness of different energies, DNA damage yields for fibroblasts are shown in Fig. S4 as a function of initial energy (symbols connected by solid lines) compared to the results for lymphocytes (dotted lines). The results for the two nuclei are in good agreement with each other; this is the case also in dependence on LET or as a function of (Z eff /) 2 (results not shown). On the other hand, the per-track yields (not shown) are lower by a factor of ~2 for fibroblasts than for lymphocytes, since the track lengths are shorter in the flat fibroblast nuclei than in the spherical lymphocytes.

Yields of DNA damage per track and lymphocyte cell in dependence on LET
Cellular DNA damage induction per track (Fig. S5) highlights its increase upon slowing down of the primary particle. For all studied endpoints the yields of DNA damage per track in lymphocyte cell in dependence on LET show continuous increase with increasing LET until the particle is stopped inside the nucleus, and the dependence on particle type is rather low as long as the particles possess energy sufficient to pass through the whole nucleus.

Characterizing ion tracks by (Z eff /)²
A commonly used parameter for radiation quality 10 that may provide a lower dependence on the particle type than the LET is the square of the effective charge divided by the velocity (Z eff /) 2 .
Similarly to the LET discussed in the main body of the paper, also (Z eff /) 2 is approximately constant for high-energy ions but for low-energy ones it varies significantly along their traversal through the cell nucleus. To deal with this issue, values of (Z eff /) 2 have been determined in this work for the mean energy inside the nucleus, which in turn has been calculated from the mean energy loss during the passage of the particles through the nucleus. Note that these values of (Z eff /) 2 hardly represent the track-average values for low-energy ions.
For comparison of the LET and (Z eff /) 2 concepts, the values of (Z eff /) 2 and of the LET resulting from dose and fluence, Eq. (4), for these mean energy values are presented in Fig. S6. At low energies, considering the LET values in dependence on these mean energies (Fig. S6a) leads to an improved agreement with ICRU stopping power data 2,7 than if the LET is related to the initial energies (Fig. 2a).
For the determination of Z eff , the Barkas formula in Eq.
(2) has been used for the sake of consistency with published data; the general findings would not change if the denominator in Eq. (1) were included.
The resulting (Z eff /β)² data span a range from ~2.5 for the fastest H up to ~50,000 for the slowest Ne particles (Fig. S6b). Contrary to the (local) LET and stopping power which increase in the course of an ion's slowing down from an initial value to a maximum (single-particle Bragg peak) and then decrease in its distal part, (Z eff /β)² steadily increases for each particle with decreasing energy. While the results plotted in dependence on LET often exhibit hooks (two distinct values at a single LET corresponding to proximal and distal parts of the single-particle Bragg peak, respectively), (Z eff /β)² can be used as a parameter for describing the results as a mathematical (single-valued) function; the results are shown in subsequent sections. In Fig. S7, yields of double-stranded DNA fragments in four size intervals formed by single ion traversals are presented as a function of (Z eff /)². The yields for H are lower than for the other particle types, the yields for He are similar to heavier particles, and the alignment of the heavier particles is at least partly better than when plotted in dependence on LET (Fig. 4). When considered per ion track, the DNA damage yields in lymphocyte cell (Fig. S8) as a function of (Z eff /)² are for H and He significantly lower than for heavier particles; the overall good alignment seen in dependence on LET for all considered particles (Fig. S5) has been lost. Uniqueness for particles heavier than He is still present up to a (Z eff /)²-value of about 10,000, but also reduced compared to the dependence on LET (Fig. S5).

RBE in dependence on LET or (Z eff /) 2
In Fig. S9, yields of DSB, DSB sites, DSB clusters and a comparison to experimental results are presented in terms of RBE. Since no calculation for a reference photon irradiation has been performed in this work, the yields have been related to the results for 128 MeV hydrogen which turned out to yield for the endpoints considered here the lowest values of all studied radiation qualities. The RBE for DSB induction (panel a) peaks at ~3 for He but saturates for heavier particles at ~2.5, which is also the maximum value for H; the MCDS prediction 11 is slightly above the yields for He and H as long as these particles possess energy high enough to fully pass through the nucleus, but do not reflect the saturation trend for the heavier particles beginning above a (Z eff /)²-value of ~1,000.