A method to implement inter-track interactions in Monte Carlo simulations with TOPAS-nBio and their influence on simulated radical yields following water radiolysis

We used TOPAS-nBio version 1.0 (Schuemann et al 2019) with TOPAS version 3.7 (Perl et al 2012, Faddegon et al 2020) to perform Monte Carlo track structure simulations based on the open access Monte Carlo Code GEANT4/GEANT4-DNA version 10.05.p01 (Agostinelli et al 2003, Allison et al 2006, Incerti et al 2010a). Based on GEANT4, TOPAS applies the same physics modules and cross section data as in the natural GEANT4 code. However, in the toolkit TOPAS various simulation configurations are already predefined, which can be accessed via simple text-based commands. This way, the toolkit enables the application of the Monte Carlo Code even for users without advanced programming skills or experience. TOPAS has already been well-validated against experimental data (Perl et al 2012, Testa et al 2013).

While TOPAS was initially designed for simulations at the macroscopic scale, the extension TOPAS-nBio allows simulations on the nanoscopic scale since particles are tracked step-by-step down to kinetic energies of some eV. ⁴ Thus, the radiation effect can be studied at cellular and sub-cellular levels. In this study, the interactions of primary particles as well as chemical radicals and species in the vicinity of the DNA are of particular interest as this is the key component of radiobiological effects. For radiobiological simulations, a variety of biological structures such as different DNA models, cell nuclei and cell models are included in TOPAS-nBio. McNamara et al (2018) have published a detailed description of these geometries. Along with these predefinded geometry components, TOPAS provides several scorers for classifying the DNA damage by e.g. single strand breaks (SSB) and double strand breaks (DSB).

In addition to the extended physics and geometrical features, chemical interactions of the radiolysis products can be simulated (Karamitros et al 2011). Thus, the indirect damaging of the DNA can be investigated. An even more accurate simulation of the expected biological response can be achieved by combining DNA repair models with the simulation. Currently, two different mechanistic repair models can be applied for modeling biological end points like cell survival curves in conjunction with the track structure simulations performed with TOPAS-nBio (McMahon et al 2016, Warmenhoven et al 2020).

In previous studies, validation of TOPAS-nBio was performed and the radiobiology extension of TOPAS has been proven to be an accurate, well performing Monte Carlo tool regarding radiolysis and DNA damage simulations. McNamara et al (2017) showed a good agreement of simulated yields of direct SSB and DSB on simple DNA models with other simulation results and experimental data. Ramos-Méndez et al (2018) performed chemistry simulations and compared the resulting yields of chemical products following water radiolysis with other Monte Carlo simulation results and experimental data. Furthermore, TOPAS-nBio was used in several studies to simulate the DNA damage induced by ionizing radiation (Zhu et al 2020a, 2020b, Ramos-Méndez et al 2021). Another application of the toolkit are investigations of irradiation including gold nanoparticles for dose enhancement (Rudek et al 2019, Hahn and Villate 2021, Klapproth et al 2021). Moreover, the radiobiology effect of gadolinium dose enhancement in neutron capture therapy was studied using TOPAS-nBio by Van Delinder et al (2021).

The purpose of this study was to examine the yield of chemicals, the so called G-value, in dependence of inter-track interaction. According to Karamitros et al (2011), the G-value G is defined as the total number N(t) of chemicals at a given time t produced or consumed per 100 eV of total deposited energy E in an irradiated system:

$$\begin{eqnarray}&&G=\displaystyle \frac{N(t)\cdot 100}{E\,({eV})}.\end{eqnarray} \tag{ 1 }$$

Since, per default, it is not possible to model inter-track interactions TOPAS-nBio, we first developed an approach for modeling inter-track interactions using this Monte Carlo toolkit.

In general, in TOPAS-nBio the radiation effect, i.e. the physical and chemical stage, of one primary particle, or more precisely of one history, is simulated in one event after each other and independently from each other. That means, at the chemical stage, reactions between the chemical species produced by the same primary particle and its secondary particles are taken into account. However, no chemical reactions between chemicals produced in two or more different events (i.e. primary particles) can be considered, so that inter-track interactions are by default not possible in TOPAS-nBio. In order to enable inter-track interactions using TOPAS-nBio, we benefited from the fact that TOPAS-nBio is able to simulate the reactions between chemical species produced by the same primary particle and its secondary particles, i.e. all chemicals with the same event number. That means to allow inter-track interactions of chemicals produced by N primary particles and its secondaries, we simulated N original histories in one event. This way the radiation effect of each history is not simulated one after the other as per default but N histories are simulated simultaneously. By varying the number N of histories simulated simultaneously in one event, we could directly handle a factor proportional to the particle fluence and dose rate. In brief, in order to realize the simulation of N histories in one event in TOPAS, we scored a phase space file of the source and modified. A more detailed description of the implementation of our method in TOPAS is given in appendix A.

All in all, to perform simulations in TOPAS-nBio with inter-track interactions according to our approach, we followed these three steps:

• 1.  
  Scoring a phase space file of the original primary particle source.
• 2.  
  Modification of the phase space file so that N histories are simulated simultaneously in one event.
• 3.  
  Using the modified phase space file in further simulations to investigate water radiolysis and inter-track interactions by examining the yield of chemical radicals.

We applied our approach for realizing inter-track interactions to three different particle sources, which are characterized in the following subsections. First, we investigated the effect of inter-track interactions in detail using a simple source of electrons of low energy (60 eV). With these simulations, we aim to obtain fundamental and basic understanding of the impact of inter-track interactions on the chemical stage of radiobiological simulations. Using this low energy, we compromised between simulation time and detailed track information for analysis. Second, once the basic impacts are identified, we applied our approach on proton beams of two different, higher energies (10 and 100 MeV).

All simulations were performed in G4 _WATER, which has a density of 1 g cm⁻³ and a mean ionization potential of 78 eV. Reference simulations, in which no inter-track interactions take place, were performed without modifying the phase space file.

2.3.1. Particle source a: electron source

In the first simulation scoring the phase space file to enable inter-track interactions, 60 primary electrons with an energy of 60 eV were generated using a volumetric source in TOPAS. This energy was chosen since on the one hand it is the mean energy produced by primary ions of clinically relevant energies (Pimblott and LaVerne 2007) and on the other hand simulation times are acceptable, due to less chemical species produced by only a few inelastic interactions. The last point is particularly crucial, since with increasing N, simulation times grow due to a larger amount of chemical species that interact via diffusion controlled reactions with each other. In figure 1(a), a schematic illustration of the geometry set-up using electrons is depicted. The electrons are emitted in a spherical shell with a thickness of 0.1 nm (reason explained below) and an outer radius r_e. In a distance of 0.1 nm from this spherical shell, the phase space file is scored on the outer curved surface of a sphere with radius r_phsp = r_e + 0.1 nm.

Zoom In Zoom Out Reset image size

In order to investigate the influence of inter-track interactions on both the time component, here realized by different N, and the spatial component, we set r_phsp to the following different values: r_phsp= 1, 10, 20, 30, 40, 50, 75, 100 nm. Since the electron density and hence the density of chemical species gets influenced by the source geometry, we expected a change of the G-value due to different r_phsp. To ensure inter-track interactions, we started with the densest arrangement of electrons and then systematically decreased it. When scoring electrons at r_phsp, we applied filters to the scorer so that only primary electrons with a kinetic energy of 60 eV at the scorer surface were included in the phase space file. This way, it can be avoided that the electron has already undergone interactions like ionization or excitation resulting in chemical species in the following radiolysis. These molecules would not be included in the subsequent simulation calculating the G-value. Thus, the thickness of the spherical shell was set to 0.1 nm, which is the smallest geometrical dimension that can be set in TOPAS-nBio. This way it can be avoided in advance that the particles scoring on the phase space surface already undergone inelastic interactions in this geometry component.

After the first simulation, the phase space file was modified in order to enable inter-track interaction of N = 2, 4, 5, 6, 10, 12, 15, 20, 30, 60 primary particles. These factors were chosen so that the same total number of 60 primary particles could be simulated in each simulation regardless of N. For example, choosing N = 2, chemicals produced by two different primary electrons can interact with each other. For a consistency of presentation, the reference simulations are specified as N = 1.

In the second simulation scoring the G-value, this modified phase space file was set as the particle source. The G-value was scored in a sphere with the radius of r_G = 1 μm for all simulation set-ups of r_phsp and N. The chosen radius ensures that all electrons and their chemicals are stopped inside this volume.

The simulations were repeated 100 times with different random seeds and final G-values were received by averaging the G-values of these 100 independent simulations.

In brief, with this simulation set-up, we studied the effect of inter-track interactions on the G-value by controlling the density of chemicals the following two ways:

• By varying r_phsp, i.e. the geometrically setup, we varied the average distance between the initial electrons.
• By varying N, i.e. varying the setup on the time scale, we varied the number of tracks between which inter-track interactions can occur.

2.3.2. Particle source B + C: proton sources

In the physical stage, interactions of all primary and secondary particles with the surrounding water molecules were simulated step-by-step using the GEANT4 physics constructor G4EmDNAPhysics _option2, which is currently the default setup in TOPAS-nBio. This constructor is in detail described by Incerti et al (2018) with regard to the included physics models and was further investigated in TOPAS-nBio in our previous work (Derksen et al 2021). The following pre-chemical stage, in which the chemical species are generated, was characterized by the branching ratios and dissociation schemes, which are specified by Ramos-Méndez et al (2018) and shown in table 1. Thereby, the placement of all dissociation products is specified for each product and dissociation channel in relation to its inelastic process generating the chemicals (Bernal et al 2015). In the chemical stage, diffusion of the produced chemical species were simulated step-by-step using the Brownian motion. All chemical reactions between the molecules including their reaction rate, which were implemented in GEANT4 by Karamitros et al (2011, 2014) and extended by Ramos-Méndez et al (2018), are shown in table 2. To access these chemistry models, we used the TsEmDNAChemistryExtended module in TOPAS-nBio. In table 2, the included reactions of all eleven chemical species considered in the extended chemical list are shown. Even though oxygen is specified as a molecule in table 2, dissolved oxygen is not simulated in our simulations. The time end for chemical interactions was set to 1 μs which is the maximum value in TOPAS-nBio.

Table 1. Dissociation schemes implemented in TOPAS-nBio. Adapted from Ramos-Méndez et al 2018. © 2018 Institute of Physics and Engineering in Medicine. All rights reserved.

Process			Probability (%)
Ionization state	Dissociative decay	H3O+ + OH•	100
A1B1	Dissociative decay	OH• + H•	65
	Relaxation	H2O + ΔE	35
B1A1	Auto-ionization	${{\rm{H}}}_{3}{{\rm{O}}}^{+}+{}^{\bullet }\mathrm{OH}+{{\rm{e}}}_{\mathrm{aq}}^{-}$	55
	Auto-ionization	•OH+ •OH + H2	15
	Relaxation	H2O + ΔE	30
Rydberg, diffuse bands	Auto-ionization	${\rm{H}}3{{\rm{O}}}^{+}{+\,}^{\bullet }\mathrm{OH}+{{\rm{e}}}_{\mathrm{aq}}^{-}$	50
	Relaxation	H2O + ΔE	50
Dissociative attachment	Dissociative decay	OH−+ •OH + H2	100

Table 2. Chemical reactions and reaction rates k considered in the chemical stage used in the module TsEmDNAChemistryExtended. H₂O molecules are not listed in the reaction formulas and no product means that the reaction product is H₂O. In this context, ${{\rm{e}}}_{\mathrm{aq}}^{-}$ describes an electron solvated in water. Moreover, an electron generated in the physical stage with an energy lower than the threshold of physical interactions (here 7.4 eV) is also considered as solvated. This way, its diffusion and reaction is then simulated in the chemical stage.

No.	Reaction			k (1010/M/s) b
1 a	${{\rm{e}}}_{\mathrm{aq}}^{-}+{{\rm{e}}}_{\mathrm{aq}}^{-}$	⟶	2OH− + H2	0.647
2 a	${{\rm{e}}}_{\mathrm{aq}}^{-}{+}^{\bullet }\mathrm{OH}$	⟶	OH−	2.953
3 a	${{\rm{e}}}_{\mathrm{aq}}^{-}+{{\rm{H}}}^{\bullet }$	⟶	OH− + H2	2.652
4 a	${{\rm{e}}}_{\mathrm{aq}}^{-}+{{\rm{H}}}_{3}{\rm{O}}$	⟶	H•	2.109
5	${{\rm{e}}}_{\mathrm{aq}}^{-}+{{\rm{H}}}_{2}{{\rm{O}}}_{2}$	⟶	OH−+ • OH	1.405
6 a	•OH+•OH	⟶	H2O2	0.475
7 a	•OH + H•	⟶	No product	1.438
8 a	H• + H•	⟶	H2	0.503
9 a	H3O + OH−	⟶	No product	11.031
10 a	H2+•OH	⟶	H•	0.0045
11	•OH + H2O2	⟶	HO2	0.0023
12	•OH + HO2	⟶	O2	1.0
13	${}^{\bullet }\mathrm{OH}+{{\rm{O}}}_{2}^{-}$	⟶	O2 + OH−	0.9
14	•OH + HO−2	⟶	HO2 + OH−	0.9
15	${{\rm{e}}}_{\mathrm{aq}}^{-}+{\mathrm{HO}}_{2}$	⟶	HO−2	2.0
16	${{\rm{e}}}_{\mathrm{aq}}^{-}+{{\rm{O}}}_{2}$	⟶	${{\rm{O}}}_{2}^{-}$	1.9
17	${{\rm{e}}}_{\mathrm{aq}}^{-}+{\rm{O}}{2}^{-}$	⟶	OH− + HO−2	1.3
18	H• + H2O2	⟶	•OH	0.01
19	H• + HO2	⟶	H2O2	2.0
20	H• + O2	⟶	HO2	2.0
21 a	H• + OH−	⟶	${{\rm{e}}}_{\mathrm{aq}}^{-}$	0.002
22	${{\rm{H}}}^{\bullet }+{{\rm{O}}}_{2}^{-}$	⟶	HO−2	2.0
23	${{\rm{H}}}_{3}{\rm{O}}+{{\rm{O}}}_{2}^{-}$	⟶	HO2	3.0
24	H3O + HO−2	⟶	H2O2	2.0
25	HO2 + HO2	⟶	H2O2 + O2	0.000076
26	${\mathrm{HO}}_{2}+{{\rm{O}}}_{2}^{-}$	⟶	O2 + HO−2	0.0085

Notes.

^a These reactions can occur directly after the pre-chemical stage. ^b M = 1 mol dm⁻³.

3.1.1. Radical yield at the end of the chemical phase in dependence of inter-track interactions

In figure 2, the G-value of ^•OH at the end of the chemical stage (t = 1 μs) is shown for all different r_phsp in dependence of the number N of histories simulated simultaneously in one event using electrons with an energy of 60 eV. We chose ^•OH as the molecule of interest for this comparison because it is assumed that indirect DNA damages resulting from interactions of chemical species with the DNA sugar phosphate backbone are mainly caused by ^•OH radicals (Chapman et al 1973, Achey and Duryea 1974, Roots and Okada 1975). For all r_phsp, the G-value of ^•OH reduces with increasing N (see figure 2). However, the magnitude of the decrease depends on r_phsp. While the decrease of the G-value is approximately exponential for small r_phsp, for larger r_phsp (approx. r_phsp ≥ 40) it is rather linear. The reason is, the smaller r_phsp, the higher the density of electrons and the larger the probability of chemical reactions between molecules of different primary particles. For large r_phsp, the spatial distance between primary particles is on average larger. Correspondingly, chemical species have to diffuse farther to undergo inter-track interactions. Hence, the effect of increasing N is less pronounced. Indeed, along with N, the geometry considerably effects the density of primary electrons and, consequently, the concentration of chemical radicals. This way, N and the initial, spatial distribution of the electrons affect the dynamics and complexity of the whole chemical stage significantly. Thus, we also analyzed the G-value of all six most observed molecules (H₃O , H₂O₂ , OH⁻ , ^•OH , e_aq , H^• , H₂ ) in dependence on N in addition to ^•OH . Therefore, in figure 3, the G-value of these chemicals is presented as a function of N for the smallest r_phsp (r_phsp = 1 nm) and largest r_phsp (r_phsp = 100 nm) at the end of the chemical stage. For both radii, a reduction of the chemical species e_aq , H₃O and ^•OH with increasing N can be observed. In turn, the number of the species OH⁻ , H₂ and H₂O₂ increases with increasing N. For example, comparing the reference simulation (N = 1) and N = 60 using r_phsp = 1 nm, the G-value of H₂O₂ is twice as high. Using r_phsp = 100 nm, the amount of H₂O₂ increases only by around 10 % in comparison to the reference. However, for H^• , no significant difference can be observed independent on the geometry and N. Furthermore, the G-value of the molecules H₂O₂ , ^•OH and OH⁻ goes into saturation for high N, which means that further increasing N has no effect on these G-values (see figure 3). Transferring this to FLASH-RT would mean that a maximum dose rate could exits, above which an increase of the dose rate does not have an influence on the creation of chemical species in case of the inter-track interactions anymore.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

However, the total G-value of all molecules decreases with increasing N. This can be explained by a significant increase of chemical reactions with an increase of N which is illustrated in figure 4 (left panel). For example, a reaction that occurs very frequently is reaction no. 6 in which H₂O₂ is formed from two ^•OH molecules (see table 2). This explains why the G-value of ^•OH decreases rapidly with N (see figure 3). Even though the total number of chemical species decreases with increasing N, figure 3 shows that the G-value of some molecules increase and the changes in the G-values depend on the molecule type. The reason for the different modifications of the G-value of different molecules is that some reactions occur more frequently than others with increasing N. This leads to a greater ratio of chemical reactions, that are consuming each molecule type, and reactions, that are producing each molecule type. This variable ratio in turn is responsible that the G-values of some molecules vary stronger with N than the others. The difference of the number of product and educt reactions of each chemical species is illustrated in figure 4 (right panel). Regarding the reactions of ^•OH exemplary, there is a negative and decreasing difference of reactions producing and consuming ^•OH molecules with increasing N. This means that the number of ^•OH molecules that are consumed in chemical reactions increases stronger with N than the number of ^•OH molecules that are produced in chemical reactions ⁷ . Hence, this results in a decrease of the G-value of ^•OH with N at the end of the chemical stage (see figure 3).

Zoom In Zoom Out Reset image size

All in all, the ratio of the chemical reactions consuming or producing the same chemical species is the reason why the G-value of some chemicals increases with N and for other molecules it decreases or even remains rather unchanged.

3.1.2. Time resolved radical yields

In order to validate that the diffusion ranges are responsible for the influence of the spatial separation of the primary particles, we investigated the time-dependent G-values of different r_phsp, which are shown in figure 5 for ^•OH and H₂O₂ using r_phsp = 1 nm and r_phsp = 100 nm. In general, the amount of ^•OH reduces with time for all different set-ups and more H₂O₂ is produced towards the end of the chemical stage. For the reference simulations, this is generally in agreement with the time-evolution of the G-value calculated by Ramos-Méndez et al (2018), ⁸ who validated the chemistry included in TOPAS-nBio against other simulation studies and experimental data.

Zoom In Zoom Out Reset image size

In figure 5, using r_phsp = 1 nm, the G-value of ^•OH decreases for all N with time. Already at the beginning of the chemical stage, the decrease of the G-value is more pronounced for higher N than for smaller N. The reason for this is, that for high N, the initial density of chemicals is very high and the chemicals do not need to diffuse large distances to meet a reactant. Thus, chemical reactions can already occur at the very beginning of the chemical stage (t ≤ 1 ps) explaining the differences of the G-value for different N. To quantify this, we calculated the number of reactions that occur at t ≤ 1 ps in dependence of N, which is in detail described in appendix C. All in all, the number of reactions at t ≤ 1 ps increase with N (see figure C2) explaining the differences of the G-value for different N at t = 1 ps. At the end of the chemical stage, using r_phsp = 1 nm, the G-value of ^•OH saturates for all N. Thus, this indicates that using this setup, the G-value of ^•OH does not change much with time after 1 μs. For r_phsp = 100 nm, the profile of the G-value of ^•OH is unchanged for N = 1 (see figure 5). This is because no inter-track interactions are possible and thus the density of the electrons does not matter. While for r_phsp = 1 nm and N ≥ 1 variations of the G-value are already present at the beginning of the chemical stage, for r_phsp = 100 nm changes due to inter-track interactions only become apparent from about 10⁴ ps onwards. After this time period, the G-value changes for different N, which indicates that the chemicals of different primary particles diffused far enough so that inter-track interactions can occur. This can be further illustrated by calculating the theoretical mean distance between two primary particles which is approximately 1.3 nm for r_phsp = 1 nm and approximately 132.3 nm, 100 times larger, for r_phsp = 100 nm.

Since H₂O₂ is created by a reaction of two ^•OH molecules, the relation of the time-dependent G-values relative to N and to the different electron distributions can be explained similar to the time-dependent G-value of ^•OH . For r_phsp = 1 nm, the G-value of H₂O₂ increases more strongly for higher N than for smaller N (see figure 5). For small N, the G-value is saturated at the chemical time end, whereas for N ≥ 10 the G-value decreases after about 10³ ps to 10⁴ ps. We analyzed this characteristic profile in more detail by comparing the number of reactions in which H₂O₂ is consumed and produced. We could show that the varying ratio between consuming and producing H₂O₂ reactions is responsible for the changes in the profile with increasing N. For more details, see appendix C.

For the less dense electron distribution (r_phsp = 100 nm), the time-dependent G-value of H₂O₂ does not vary much with N similar to the G-value of ^•OH (see figure 5). This can be explained, as mentioned previously for ^•OH, by a low number of inter-track interactions due to large spatial separations of the primary tracks.

All in all, analyzing the time-dependent G-value of different r_phsp, we could indicate that the variations in the G-values for different source geometries is related to the density of chemical species and the complexity of chemical reactions.

Once the fundamental influence of inter-track interactions has been analyzed and explained using a low energy electron source, we further examined the effects of inter-track interactions of high-LET radiation sources. In figure 6, the G-value at the end of the chemical stage at t = 1 μs of different chemical species in dependence of the number N of histories simulated simultaneously in one event is shown for the 10 MeV (top) and 100 MeV (bottom) proton beam. Similar to the electron source, the total number of molecules (dashed line) decreases with increasing N for both proton sources. As before, the G-value of e_aq , H₃O and ^•OH decreases with an increasing N. In turn, the number of the species OH⁻ , H₂ and H₂O₂ increases with increasing N. For 100 MeV protons, the G-value of H^• stays constant as it is the case for the electron source. However, for 10 MeV protons, the G-value of H^• decreases slightly with increasing N. For the proton sources, saturation effects of the G-value can also be observed for high N for some molecules, such as ^•OH and OH⁻ . However, using the high LET (10 MeV) proton source, these effects occur already at lower N compared to the low LET (100 MeV) source.

Zoom In Zoom Out Reset image size

When comparing the G-values of the proton sources without inter-track interactions (N = 1), it is obvious that the G-values of the two proton sources are already different. For example, the G-value of ^•OH for N = 1 using the 10 MeV source is 1.92 ± 0.01, while it is 2.53 ± 0.04 applying the 100 MeV source. This is consistent with other studies in which an LET-dependence of the G-value was observed (Burns and Sims 1981). Additionally, for the high LET source the G-value of ^•OH decreases in a greater extend with increasing N than the G-value of ^•OH using the lower LET proton source. Since the density of secondary electrons is higher in the higher LET radiation, and thus also the density of chemical species of different primary particles, here, more inter-track interactions occur in the chemical phase with increasing N, which leads to the large variation of the G-values.

The variation of the G-value of e_aq , ^•OH , H₂ and H₂O₂ resulting from inter-track interactions (see figure 3) was also observed by Kreipl et al (2009) in a similar way simulating inter-track interactions of He²⁺ ions, C⁶⁺ ions and protons. Although they investigated inter-track interactions of only two primary particles in comparison to simulations without inter-track interactions, comparable to N = 1 and N = 2 in our simulations, inter-track interactions of numerous secondary particles contribute to the considerable change of the G-value. However, systematically increasing the number of inter-track interactions as done by our work group, illustrates effects of different particle and dose rates as described in section 2.2. The higher a dose rate, the more inter-track interactions can occur since an overlap of chemical tracks of different primary particles becomes more probable. Comparing this work with simulation studies that directly induced dose rates by varying the initial time of primary particles in a fixed pulse width, the same trend in the change of the G-value can be observed: Ramos-Méndez et al (2020) simulated different dose rates by simulating a constant number of primary particles delivered in pulses of different width. A maximum change in the G-value was found for a zero width pulse, which is comparable to our method in which all 60 primary particles were initialized at the same time. Even larger pulse widths of 1 ns, 1 μs and 10 μs showed significant changes in the G-value compared to conventional dose rates. Nevertheless, larger pulse widths, which are equivalent to lower dose rates, had less effects on the G-value which is in agreement with our results allowing less inter-track interactions. Comparable results were obtained by Lai et al (2021). In their simulation study, higher dose rates contributed to a smaller amount of ^•OH molecules. While the results of the simulation studies show a similar trend of the G-values with an increasing number of inter-track interaction or the dose rate, experimental results are not straightforward. On the one hand, Kusumoto et al (2020) observed a decrease of ^•OH with increasing dose rate, as investigated in the simulation studies, but on the other hand, a decrease of H₂O₂ with increasing dose rate was observed in some experimental trials (Montay-Gruel et al 2019, Blain et al 2022). This is in direct conflict with the results of our study and the mentioned simulation studies. Here, an increase of H₂O₂ has been observed for FLASH dose rates in comparison to conventional dose rates. Differences could, for example, be caused as the simulations mostly run only up to 1μs after irradiation, whereas the experimental measurements are recorded at a later point in time. In addition, dissolved molecular oxygen is not included in some of the simulations. Since ^•OH often forms a compound with oxygen, less H₂O₂ can thus be generated by ^•OH . However, since the amount of oxygen decreases for ultra-high dose rates, less ^•OH radicals might react with oxygen in comparison to low dose rates resulting in a higher production of H₂O₂ . For a more detailed discussion of the contrary observations of H₂O₂ with increasing dose rate of experimental and simulation studies we here refer to commentary by Wardman (2020).

As already observed when regarding the G-value of ^•OH for all r_phsp (see figure 2), the variation of the G-value of all chemical species with N shown in figure 3 is more pronounced for r_phsp = 1 nm than for r_phsp = 100 nm. The influence of the spatial separation between interacting tracks was also investigated by Kreipl et al (2009) using primary He²⁺ ions, C⁶⁺ ions and protons. In contrast to this work, whereby the spacing of the interacting tracks is not systematically controlled but influenced by the source radius r_phsp, the spatial separation in their research was systematically parameterized in the range of 0 nm and 1000 nm. Nevertheless, for all particle sources, the G-value was lowest for the smallest spatial separation of the interacting tracks as seen in our results (see figure 2). In addition, Kreipl et al observed that the number of interacting tracks had no influence on the yield of radicals for a spatial separation of 1000 nm. In comparison, in our study, the spatial separation for the largest radius (r_phsp = 100 nm) cannot be larger than 200 nm, which would only be the case for N = 2 if the particles are emitted directly opposite on the sphere surface ⁹ . To investigate whether the saturation effect regarding the spatial separation can be reproduced with our simulations, we simulated the G-value using a radius of r_phsp= 750 nm, corresponding to an approximately average spatial separation of 1000 nm of the initial electrons placed on the surface of the sphere. As shown in figure B1 in the appendix, varying N in the range of 1–100 has no influence on the G-value. Generally, the effect can be explained due to the fact that diffusion ranges of the chemical species up to 1000 nm are extremely improbable as shown by Ramos-Méndez et al (2018). Thus, chemical species of different primary particles are unlikely to react with each other at such large spatial separations. However, if N would be increased by a multiple more than 100, inter-track interactions could be triggered by the high density and thus smaller distances of the primary tracks. However, the systematic arrangement of the various primary tracks as in this study and in that of Kreipl et al (2009), do not reflect the situation of realistic extended beam fields. For particles sources with a different source geometry as well as with a different LET and hence, a corresponding different spatial arrangement of inelastic interactions, the situation of inter-track interactions might be different. However, in case of simulation time we investigated the spatial separation using electrons of low energy and not using the proton sources showing exemplary that the spatial arrangement of the source geometry has an effect on inter-track interactions. The results illustrate that the denser the primary tracks are, the more likely inter-track interactions occur reducing the amount of chemicals which is illustrated by a reduced G-value. In summary, the occurrence of inter-track interactions can be controlled by the spatial arrangement of the source geometry and by the number N of primary particles simulated simultaneously in one event, i.e. broadly the dose rate, all influencing the spatial distribution of the physical particle track.

In this study, inter-track interactions of chemicals produced by different primary particles were investigated whereas the primary particles were initialized simultaneously, that means at the same initial time. In fact, using ultra-high dose rates as in FLASH-RT, the particles are delivered at different initial times in a given pulse width as simulated for example by Ramos-Méndez et al (2020) and Lai et al (2021). However, simulating all interacting tracks simultaneously, as done in this study, actually has certain advantages: first, the impact of inter-track interactions on their own can be analyzed without other parameters affecting the results. Parameters could for example be the time structure of the pulses, the dose-per-pulse, the oxygen level and the type of the surrounding tissue cells, since these parameters obviously have an influence on the dynamics of chemical reactions. However, we first wanted to show that we are able to simulate inter-track interactions in TOPAS-nBio and investigate their influence on a first principal basis. Second, Abolfath et al (2022) pointed out that in TOPAS-nBio inter-track interactions are underestimated. By simulating all tracks simultaneously, the maximized effect of inter-track interactions can be studied.

Due to limitations in our study, a quantitative relation between N and a dose rate cannot be established. One the one hand, the temporal profile of the G-value (see figure 5) clearly illustrates that the time factor is one of the key factors for the G-value. Since in our study the particles are simulated simultaneously, that means without a separation in time, a conversion to a realistic dose rate is not feasible. One the other hand, and most importantly, the simulations are performed on microscopic scale which disables a conversion into dose, a macroscopic quantity, as discussed in detail by Abolfath et al (2020). For the purpose of calibrating the number N of particles simulated simultaneously against dose rates, the simulations will be adapted in a subsequent study. Simulations will be performed under the same conditions as experiments so that by comparing both G-values a reference to clinical dose rates can be obtained.

In this initial study regarding inter-track interactions we did not considered the role of oxygen in the simulations, which will, in turn, be taken into account in subsequent studies. Additionally, in further steps, scavengers should be included in the simulations since Ramos-Méndez et al (2021) showed that including scavengers in the chemistry simulations improves the agreement to experimental data. To improve simulations of the chemical stage and linking them to experimental conditions, the GEANT4 developers are constantly updating the code with new molecules and scavengers (Chappuis et al 2023). In addition, Wardman (2022) clarified the importance of scavenger systems since in biological tissue intracellular molecules can compete with the chemical reactions in 'pure water'. This could change the effect of inter-track interactions. But as we could show that inter-track interactions occur also in the picosecond regime, we expect that the effect of inter-track interactions is still significant in the presence of scavengers as in cellular mediums. Nevertheless, this should be investigated in further studies including the mechanisms of living cells. Until now, however, there are no established simulation codes that can cover the needs of those simulations as pointed out by Wardman (2022) and will be a quite challenging task.

Furthermore to the presented results, since we could show that inter-track interactions influence the yield of chemicals, we now apply our approach to investigate the resulting DNA damage yield in more realistic simulation set-ups e.g. including a nucleus model and optimize the spatial and temporal arrangement of initial particles.

As an outlook, TOPAS-nBio has been updated to version 2.0, which provides the simulation of inter-track interactions. Therefore, we will evaluate in a following study if comparable results can be achieved by simulations applying our approach and simulations using TOPAS-nBio 2.0.

This research showed that increasing the number of histories from which the chemical species can interact with each others, representing higher dose rates, reduces the total G-value (see figures 3 and 6) and hence more chemical reactions occur. Increasing the number N of histories simulated simultaneously in one event has a comparable effect on the G-value like increasing the dose rate as observed in an experimental study by Kusumoto et al (2020). With higher dose rates, and correspondingly more inter-track interactions, a reduction of ^•OH is observed. Since ^•OH is considered to be mainly responsible for indirect damages on the DNA, we hypothesize that this has an impact on the amount and distribution of DNA damages. This in turn, could play a crucial roll in the FLASH effect. To provide an explanation of the FLASH effect, these results are surely not sufficient enough since a difference between tumor and normal tissue response cannot be specified. But as the FLASH effect is 'the result of a multi-parameter situation' (Rothwell et al 2021), each individual variation is relevant. The exact explanation of the FLASH effect is beyond the scope of this paper but inter-track interactions may be a parameter to be considered in the mechanism.