Quantitative estimation of track segment yields of water radiolysis species under heavy ions around Bragg peak energies using Geant4-DNA

We evaluate the track segment yield G′ of typical water radiolysis products (eaq−, ·OH and H2O2) under heavy ions (He, C and Fe ions) using a Monte Carlo simulation code in the Geant4-DNA. Furthermore, we reproduce experimental results of ·OH of He and C ions around the Bragg peak energies (< 6 MeV/u). In the relatively high energy region (e.g., > 10 MeV/u), the simulation results using Geant4-DNA have agreed with experimental results. However, the G-values of water radiolysis species have not been properly evaluated around the Bragg peak energies, at which high ionizing density can be expected. Around the Bragg peak energy, dense continuous secondary products are generated, so that it is necessary to simulate the radical–radical reaction more accurately. To do so, we added the role of secondary products formed by irradiation. Consequently, our simulation results are in good agreement with experimental results and previous simulations not only in the high-energy region but also around the Bragg peak. Several future issues are also discussed regarding the roles of fragmentation and multi-ionization to realize more realistic simulations.

(from 10 -15 to 10 -12 s) and the "chemical phase" (from 10 -12 to 10 -6 s). For example, in the PTSim, only energy deposition is taken into account, meaning that DNA damage induced by water radiolysis species in the chemical phase is not considered. Water radiolysis species are abundantly formed along ion tracks. To consider roles played by water radiolysis species, we have to follow the reaction occurring with each water radiolysis product. However, it takes long time for computation. In the "physicochemical phase" for water radiolysis, initial products are reorganized into modified species (e.g., stable molecules and water-free radicals) and thermal equilibrium in the bulk media is established. Reaction processes are then extended up to neutralization of the media in the "chemical phase". During the "physicochemical phase" and the "chemical phase", water radiolysis products can damage DNA molecules. The correlation between yields of water radiolysis products and those of DNA molecules damaged should therefore be assessed to elucidate the contribution of indirect action.
Many researchers have experimentally evaluated radiation chemical yields (G-values), the number of entities formed or destroyed by unit energy (conventionally 100 eV), of · OH, hydrated electrons (e − aq ) and hydrogen peroxide (H 2 O 2 ), all of them are typical water radiolysis products under ionizing radiation 6-10 , using radical scavengers (e.g., coumarin-3-carboxylic acid (C3CA) and phenol). Furthermore, the correlation between generated yields of · OH, which efficiently react with DNA molecules, and those of radiolysis product of an amino acid disassembled in solution is being quantitatively clarified 11 . As is well known, the effectiveness of indirect action increases with the increasing of the energy of incoming ions, while the linear energy transfer (LET) decreases 10,12 . Around the Bragg peak of C ions, the contribution of indirect action is up to 50% for cell killing under humidified air with 5% CO 2 at 37 ℃ 13 . Chemical reactions by water radiolysis products should therefore be taken into account to properly understand the biological effectiveness of cancer treatment with radiation 14 .
In parallel with experiments, the LET dependence of G-values of typical water radiolysis products in the high energy region (> 100 MeV/u) was evaluated using Geant4-DNA 15,16 The simulation results were in agreement with experimental results. However, G-values around Bragg peak energies have not been successfully simulated. Even at the Bragg peak energy region of C ions, the contribution of indirect action for biological effects is still large (~ 50%). We are now trying to establish a platform that can estimate the therapeutic effects and evaluate biological effects by simulations with reference to the amount of · OH obtained experimentally. In this work, we simulate G-values of · OH, e aq − and H 2 O 2 under heavy ions (He, C and Fe ions) up to 10 -6 s using Geant4-DNA. The obtained simulation results are compared to the experimental results around Bragg peak energies and previously obtained ones in the high energy region.

Materials and methods
Experiments. Previously, yields of water radiolysis products ( · OH, e aq − and H 2 O 2 ) under ion irradiation were evaluated. In this section, we reproduce the derivation process of the yields reported in the previous study 11 .
where E 0 is the initial energy of incoming ions. When the experiments are done at several energies, N can be plotted as a function of the energy of the ions, and a resulting fitting curve can be derived. Then, we can determine the track segment yield G′, that is the yield for an infinitesimal energy loss, defined for one specific energy of ion: The so-called "scavenging capacity" of a scavenging reaction is expressed by the product of the rate constant of the reaction and the concentration of the radical scavenger. Its inverse represents the average scavenging time scale. So, by adjusting the concentration of probe, the time dependence of the G-value of · OH can be reconstructed 6 . The detailed procedures for experiments to evaluate G-values of · OH are described elsewhere 10,11,17 .
The aqueous solution containing 1 mM sodium nitrate (NaNO 3 ) and 5 mM di-sodium hydrogen phosphonate (Na 2 HPO 3 ) are used for the measurement of e aq − . e aq − is scavenged by nitrate anion to produce NO 3 2−· , which immediately reacts with one of the surrounding water molecules, resulting in the production of nitrogen dioxide (NO 2 · ). NO 2 · then reacts with hydrogen phosphate anion (HPO 3 2− ) to produce an anion (NO 2 − ), which can easily be determined by applying Saltzman technique 18 with molar extinction coefficient of 42,3000 M −1 cm −1 at 540 nm. Furthermore, Ghormley technique 19 has been adopted to evaluate the yield of H 2 O 2 . H 2 O 2 is very stable compared to other water radiolysis products but is destroyed by e aq − due to intra-track reactions. To minimize the decomposition by intra-track reactions, 2.5 mM NaNO 3 aqueous solution was used for measurements. The detailed procedures for experiments are described elsewhere 8 .
In this work, we re-plotted these experimental results. All experiments were done under neutral pH conditions at 25 ℃ 7,8,10 .

Simulations.
A Monte Carlo simulation with Geant4-DNA was performed to calculate G-values of water radiolysis products [20][21][22][23][24][25][26] . We used the G4EmDNAPhysics_option8 physics constructor with additional subexcitation processes of vibrational excitation and molecular attachment for electrons installed in Geant4-DNA ver. 10.05.p01. The G4EmDNAPhysics_option8 physics constructor covers the energy range from 0.5 MeV/u to 10 6 MeV/u for heavy ions (He, C and Fe ions). Furthermore, G4EmDNAChemistry was used for the simulation in the chemical stage. The simulation geometry is a 10 × 10 × 10 mm 3 air-free water cube. This is the size determined in consideration of the measurement system. When primary particles have deposited 10 keV of their energies into the water, the charged particle tracking simulation of the physical process was ceased and shifted to the chemical process. We aborted the event when total energy deposition of each event exceeded 10.1 keV in order to simulate large enough numbers of tracks in reasonable calculation times. That means that total energy deposition of each event was always between 10 and 10.1 keV. In Geant4-DNA, we can follow a particle-based approach for the simulation of water radiolysis 27 , where molecular species are modeled as point-like objects diffusing in a continuous liquid water medium. In this work, we follow the radical annihilation process from 1 ps to 1 μs. When energies of incoming ions were high enough (e.g., 100 MeV, 25 MeV/u and 400 MeV/u for protons, He and C ions, respectively), they could easily go through a water cube of 10 mm thickness. If not, incoming ions are completely stopped in the water cube. In such cases, energies of incoming ions were rapidly lost, meaning that high ionizing density can be expected around the Bragg peak. This implied that water radiolysis products could be abundantly generated along ion tracks. So, considering the issue of calculation time, it is not easy to follow the chemical process around the Bragg peak energy using a CPU-based simulator. In this study, the issue was overcome by a split simulation. The details of the split simulation are discussed in the following section.
Substantial yields of secondary products generated by water radiolysis (e.g., O ·− , O 2 , O 2 ·− , HO 2 · , HO 2 − ) have been previously evaluated [28][29][30] . In this study, we defined additional molecular species shown in Table 1, which are not implemented in Geant4-DNA ver. 10.05.p01, considering the number of atoms, the number of occupied electrons, the electron occupancy (depending on the molecular composition of each molecule), the Van der Waals radius, mass, charge and diffusion coefficient produced in air-free water due to irradiation. The diffusion coefficients related to the newly defined products are listed in the right column of Table 1 31 . In the physical stage, the water molecules are excited or ionized, and initial radiolysis products are generated by dissociation in the physical-chemical stage 32,33 . The dissociation/association scheme leads to the formation of further radiolysis products such as H 2 , H · and · OH. The branching ratios used in this study are listed in Table 2. This table is recalled from a previous study 15 . In the chemical stage in Geant4-DNA, starting from 1 ps, the molecular species diffuse and can react together based on the diffusion coefficients and the chemical reaction rates. The list of reactions and reaction rate constants are shown in Table 3. The reaction rate constants were referred from a previous study 34,35 . The present simulation was done using G4EmDNAChemistry with the role of the newly defined molecules.

Results and discussion
Time dependence of G-value of water radiolysis products. Figure 1 illustrates the schematic view of diffusion of reactive species produced by water radiolysis from 1 ps to 1 μs after the irradiation of 400 MeV/u C ions. Water radiolysis products along the C ion path are denser than those along secondary electron trajectories. The high density of reactive species can be seen around the C ion trajectory. Reactive species then diffuse to distant locations from the ion path.  ) are shown with a solid blue line and those without the reactions are shown with a solid red line (G4EmDNAChemistry). The previously obtained reference data are also plotted [36][37][38][39][40][41][42][43] . The simulation results are in good agreement with experimental values. The G-values of · OH calculated using G4EmDNAChemistry are slightly higher than our simulation result, and closer to experimental data in several cases. This implies that the added reactions act significantly. However, there are variations in experimental results, so that it is difficult to argue if the present simulation with reactions of water radiolysis product with newly defined molecules is appropriate. In comparison with · OH, no difference is observed between the two simulations for the G-values of e aq − . Furthermore, the G-value of H 2 O 2 with the new reactions added is higher than that without the reactions. Indeed, H 2 O 2 is not directly produced by water radiolysis, meaning that H 2 O 2 is made by radical-radical reaction (e.g., · OH + · OH → H 2 O 2 ). Since the G-value of H 2 O 2 with the new reactions is higher than without them, the role played by secondary products seem to be important. Similar results are observed with 400 MeV/u C ions as shown in Fig. 2d-f, with no dramatic changes in G-values.
When incoming ions (or electrons) have high enough energy, no drastic changes in G-value is visible (Fig. 2). But, to properly simulate the experimental results around Bragg peak energies 17 , at which incoming ions are completely stopped in the experimental set-up, we should consider the role played by secondary products added. Around the Bragg peak, many water radiolysis products would be generated along the ion trajectory due to high ionization density. Because of the high ionization density, it is necessary that radical-radical reactions are properly taken into account. Figure 3 shows the G-values of · OH for 0.75 MeV/u He ions and 0.83 MeV/u C ions without the reaction of secondary products added (red: G4EmDNAChemistry). The experimental results are also plotted in black symbols 17 . The G-values of · OH without the reaction of added secondary products are about twice higher than experimental results. Around the Bragg peak energy (< 6.0 MeV/u), high ionization density is expected compared to the high-energy region and incoming ions would completely stop in the irradiation cell. It is thus very tough to follow the reactions induced by incoming ions by a conventional CPU simulation. To more accurately simulate the G-value around the Bragg peak energy, we performed a split calculation. The number of · OH produced was computed for each energy step with an increment of 0.25 MeV/u, and then, integrated from 0 to any energies as indicated in Fig. 4a,b. Finally, G-values are calculated from Eq. (2). All results shown in Fig. 4 have been calculated with the new reactions added. The time was chosen close to the scavenging ones of experimental data. When the role of reactions of added secondary products (Table 1) are taken into account, the discrepancy between simulation results (blue) and experimental results becomes significantly smaller than that without the reactions represented in Fig. 3. Figure 5 shows the number of reactions related to · OH as a function of LET of C ions, ranging from 10 eV/nm to 700 eV/nm. The number of reactions is cumulative one up to 1 μs after the irradiation. Among newly added reactions, · OH + OH − → O ・− + H 2 O and O ・− + H 2 O → · OH + OH − contribute significantly to the reaction with · OH. This means that we have to consider roles of newly added reactions to accurately simulate the G-values of water radiolysis products because these reactions contribute the annihilation process of · OH. Figure 6a-c show the LET dependence of G-values of · OH, e aq − and H 2 O 2 at 100 ns after irradiation, simulated in this work and compared to experimental data. Overall, G-values of examined reactive species are in good agreement with experimental results (plotted points) in a wide LET range including the Bragg peak energies. In the cases of · OH and e aq − , the G-value decreases monotonically with increasing LET as shown in Fig. 6a,b. In comparison to this, the G-value of H 2 O 2 increases monotonically www.nature.com/scientificreports/ www.nature.com/scientificreports/ with increasing LET, with higher values for lighter ions at the same LET as shown in Fig. 6c. These results remind us that LET is not a universal parameter for describing the G-value. For scaling G-values, Z eff /β (or (Z eff /β) 2 ), where Z eff is the effective charge of incoming ions and β is the velocity of incoming ions normalized by the speed of light in vacuum, has been applied. Indeed, the description of the dependence of G-value is improved compared to that with LET, but Z eff /β has not been recognized as a universal parameter. To universally express  www.nature.com/scientificreports/ G-values of water radiolysis products in a wide LET range, the electron interaction concept, which the number of interaction induced by secondary electrons governs the radiation induced yields could work well 7,44,45 . Furthermore, G-values of · OH for Fe ions by the present simulation are significantly lower than the experimental results. In the experiments, a water equivalent moderator was used for adjusting the incident energy of Fe ions. Under this process, lighter ions, whose LET is smaller than that of Fe ions, could be produced by fragmentation. Considering the influence of lighter ions produced by fragmentation, higher G-values are expected, so it is reasonable that the present simulation result is lower than those of experiments. In the present version of Geant4-DNA, the influence of fragmentation of the ions is not considered. The implementation of the influence of the lighter ions produced by fragmentation would be one of the future issues to properly simulate the contribution of indirect action in the human body 46 . Moreover, discrepancies between simulations and experiments are confirmed in G-values of H 2 O 2. The reason of this discrepancy is the implementation of multi-ionizations.

LET dependence of reactive species yield.
In agreement with a previous study 47 , the influence of multi-ionization of air-free water molecules act effectively. Meesunguonen and Jay-Gerin suggested that the multi-ionizations, although less frequent compared to single ionization, especially contributes to the primary HO 2 · /O 2 ·− yield in the high LET region. In Meesungunoen  www.nature.com/scientificreports/ J. and Jay-Gerin 's code named IONLYS-IRT, cross sections of elastic, phonon and vibrational electron scattering obtained from electron-impact on amorphous ice films have been employed. To take into account the effects of multi-ionization under high-LET heavy-ion irradiation the model has been extended to incorporate double, triple, and quadruple ionization process in single ion-water collisions. In the present simulation using Geant4-DNA, we consider the elastic scattering, ionization, electronic excitation and charge exchange processes for electron, proton and alpha particles. Furthermore, in the cases of heavy ions, only the ionization process is considered. In the absence of multi-ionization, the G-value of H 2 O 2 increases monotonically with increasing LET as shown in Fig. 6c. When multi-ionizations are considered, we anticipate that the G-value of H 2 O 2 would increase with increasing LET up to 100 eV/nm and then drop. Generally speaking, H 2 O 2 is formed by the recombination of two · OHs ( · OH + · OH → H 2 O 2 ). As LET increases, radical-radical reactions occur more efficiently. Above 100 eV/nm, it is known that O( 3 P) is produced by multi-ionizations process. Since · OH reacts not only ), · OHs are consumed by other radical-radical reaction, causing a reduction of yields of H 2 O 2 yield at high LET. However, we do not consider the contribution of the multi-ionization process in our simulation. Therefore, the roles of O( 3 P) are not added. H 2 O 2 is formed not only by the recombination of two · OHs, but also by the reaction between H · and HO 2 · (H · + HO 2 · → H 2 O 2 ), and the reaction between H 3 O + and HO 2 − (H 3 O + + HO 2 − → H 2 O 2 ). The previous simulation performed by Meesunguonen and Jay-Gerin more properly reproduced experimental G-values of H 2 O 2 above 100 eV/nm. The current version of Geant4-DNA does not handle the multi-ionizations process. We are planning to add it into Geant4-DNA in the future. In comparison to e aq − , the influence of multi-ionizations to the G-value of · OH is not significant. One of the main processes consuming the · OH is its interaction with e aq − . This reaction is a competing process with the interaction of e aq − with the O 2 ·− . However, for the high LET region, · OH react more efficiently with e aq −48 . In brief, · OH are diminished by reactions with e aq − in the track core. This explains why the present simulation results of · OH are in good agreement with experimental results as represented in Fig. 6a. In accordance with previous studies, · OH and e aq − are not affected by the multi-ionization 47,49 . Furthermore, · OH react with e aq − in the high LET region, dense track on the particle path 47,49 . Thus, although multi-ionization processes are not considered in the simulation, the present results regarding · OH and e aq − are in good agreement with experimental results. Figure 7 shows a comparison of the simulated G-values of · OH at 1 μs after irradiation to a previous simulation 47 and experimental measurements 50 . The present simulations (open symbols) are in good agreement with the previous simulations (solid lines) and experimental results (solid symbols). Once again, the present simulations successfully reproduced the experimental results in a wide LET range thanks to the consideration of newly added reactions.

Conclusions
In the present study, we simulated the track segment yields (G′) of water radiolysis products ( · OH, e aq − and H 2 O 2 ) under heavy ion irradiation using a Monte Carlo simulation code in the Geant4-DNA in a wide LET range up to 700 eV/nm. To accurately simulate the G′ around the Bragg peak, we addressed two issues. The first issue is roles played by secondary products generated by water radiolysis. The second issue is the long computation time due to the abundantly generated water radiolysis products around the Bragg peak. To save the computation time, we did the split simulation. Consequently, G′ of · OH and e aq − were in good agreement with experimental results in the examined LET range. However, discrepancies of G′ of · OH between simulation and experiment were seen in Fe ions. This discrepancy suggested that the contribution of lighter particles produced by fragmentation should be considered for a more accurate simulation. At the same LET, G' of · OH and e aq − of heavier ions were lower than that of lighter one, especially noticeable from 40 to 100 eV/nm. This finding is consistent with the fact that LET Figure 7. LET dependent G-values of · OH at 1 μs after irradiation of a water target with different radiation types. The solid lines presented by Meesungunoen J. and Jay-Gerin J.P. represent the results of Monte Carlo simulations incorporating multi-ionization of water. The short-dot lines presented by by Meesungunoen J. and Jay-Gerin J.P. correspond to simulated G-values calculated as a function of LET without including the mechanism of multi-ionization of water 47 . Experimental data presented by Burns W.G. and Sims H.E. 50 were used for comparison with simulated LET dependence of G-values. www.nature.com/scientificreports/ is not universal parameter to express the yields of water radiolysis products especially when it come to complex track structure with high energy secondary electron. At this stage Z eff /β could be a better parameter. It will be crucial to discuss a parameter that universally describe the G′ of water radiolysis products. In comparison to · OH and e aq − , G′ of H 2 O 2 did not agree with experimental results above 100 eV/nm even after the consideration of roles of secondary products generated by water radiolysis. This is because the contribution of multi-ionization is not taken into account in the present stage. To have a more realistic simulation, the multi-ionization should be added in the Geant4-DNA in the near future.