Comparison of two methods simulating inter-track interactions using the radiobiological Monte Carlo toolkit TOPAS-nBio

Objective. To compare two independently developed methods that enable modelling inter-track interactions in TOPAS-nBio by examining the yield of radiolytic species in radiobiological Monte Carlo track structure simulations. One method uses a phase space file to assign more than one primary to one event, allowing for inter-track interaction between these primary particles. This method has previously been developed by this working group and published earlier. Using the other method, chemical reactions are simulated based on a new version of the independent reaction time approach to allow inter-track interactions. Approach. G-values were calculated and compared using both methods for different numbers of tracks able to undergo inter-track interactions. Main results. Differences in the G-values simulated with the two methods strongly depend on the molecule type, and deviations can range up to 3.9% (H2O2), although, on average, the deviations are smaller than 1.5%. Significance. Both methods seem to be suitable for simulating inter-track interactions, as they provide comparable G-values even though both techniques were developed independently of each other.


Introduction
The radiobiological Monte Carlo toolkit TOPAS-nBio (Schuemann et al 2019), an extension of TOPAS (Perl et al 2012, Faddegon et al 2020), has recently been updated to version 2.0, introducing new features.One feature is the simulation of inter-track interactions between chemical radicals from different primary tracks and their effects on the yield of chemicals during the chemical stage (Ramos-Méndez et al 2020).
Inter-track interactions occur when the chemical stage (t 1 ps) of two or more primary particles overlaps in time, which may happen at ultra-high dose rates emitting a large number of particles in very short time intervals.This way, molecules generated through the radiolysis of water cannot only interact within their track (intra-track interactions), but can also perform chemical reactions with molecules from another track (intertrack interactions).In TOPAS-nBio version 1.0, inter-track interactions are not taken into account by default since under conventional conditions and dose rates of under 0.03 Gy s −1 in radiotherapy inter-track interactions are highly improbable (Lai et al 2021).Thus, all physical and chemical interactions of one track are supposed to be independent of other tracks and are simulated sequentially, meaning that chemicals from different primaries cannot interact with each other.However, applying new techniques in radiotherapy like the use of ultra-high dose rate, it cannot be assured that the events occur independently of each other.Even though the time scale of physical interactions is quite small (1 fs) (Ramos-Méndez et al 2018), the chemical stage of different primary particles could overlap, making inter-track interactions possible (Lai et al 2021, Baikalov et al 2022).However, in TOPAS-nBio 1.0, there is no feature included allowing for the modelling of inter-track interactions.Thus, in a preliminary work (Derksen et al 2023), we developed a method enabling inter-track interactions and investigated their effect on the dynamics of the chemical stage following radiation using TOPAS-nBio 1.0.
Recently, an updated version of TOPAS-nBio 2.0 was released, in which inter-track interactions can be modeled by default using a specific scorer.
The aim of this study is to compare these two independently developed methods for simulating inter-track interactions in TOPAS-nBio.Therefore, we examined the yield of chemicals following the irradiation of a waterfilled sphere using an isotropic electron source of 4.5 keV, modelling different amounts of inter-track interactions.

Materials and methods
2.1.The Monte Carlo toolkit TOPAS-nBio Performing track structure simulations, we used TOPAS-nBio version 2.0 with TOPAS version 3.9 to simulate physical and chemical interactions following the irradiation.Based on the open-source Monte Carlo Code GEANT4/ GEANT4-DNA (version 10.07.p03) (Agostinelli et al 2003, Allison et al 2006, Incerti et al 2010), the toolkit TOPAS and its extensions are freely accessible. 5With the TOPAS extension library TOPAS-nBio, in addition to simulations on the macroscopic scale, simulations can also be performed on nanoscopic scales, like cells and DNA.Here, not just every physical interaction of a particle track is simulated step-by-step until the primary particle has a residual kinetic energy of a few eV,6 but also the production and reactions of chemical species following the radiolysis can be simulated.Following the update of TOPAS-nBio, the code offers two methods of simulating the chemical stage: the step-by-step method (SBS) (Karamitros et al 2011) and the independent reaction time method (IRT) (Green et al 1990, Schuemann et al 2019, Ramos-Méndez et al 2020).While using the SBS method, the reaction kinetics are simulated step-by-step based on their simulated diffusion steps.Using the IRT method, the diffusion of chemicals is not calculated step-by-step, but chemical reactions are simulated by calculating reaction times based on the initial positions of pairs of radiolytic species and probability functions.Even though the spatial information of the chemicals is not provided in contrast to the step-by-step method, this approach features much higher computing efficiency.A detailed description of this approach and a comparison to the SBS chemistry is given by Ramos-Méndez et al (2020).In previous studies, both TOPAS and TOPAS-nBio have been well-validated against experimental data and other simulation studies (Perl et

Simulation of inter-track interactions
Investigating the effect of inter-track interactions was performed by calculating the G-value in dependence on the number of tracks with inter-track interactions.Here, the G-value is defined as the total number M(t) of radiolytic species that are generated or consumed in the case of chemical reactions at a given time t normalized to 100 eV of deposited energy E, formalized by Karamitros et al (2011) as: In the following, both methods to enable the modelling of inter-track interactions will be described.

Method 1: phase space approach
This method was developed in a preliminary work (Derksen et al 2023) to investigate inter-track interactions with TOPAS-nBio 1.0 using the SBS approach.However, it is also possible to apply this method using the IRT model available in TOPAS-nBio 2.0.To maintain consistent simulation conditions for both methods, we utilized TOPAS-nBio version 2.0 in both cases.This approach is based on a manipulation of the particle source using a phase space file.Here, secondary particles produced by different primary particles are appointed to one primary particle.This is achieved by modifying the label 'Flag to tell if this is the First Scored Particle from this History' in the phase space file.This label can be 1 or 0. A label of 1 identifies the first particle scored of a given history, which is not necessarily a primary particle depending on the settings of the scorer.A label of 0 signifies particles belonging to the same history as the previously scored particle.In our case, only primary particles are scored, meaning that all particles are initially labeled with 1.By changing the label of a specific number of sequential particles to 0, those particles are considered as secondary particles from the previous particle with label 1.Thus, all these particles are simulated simultaneously including the subsequent created chemical species enabling inter-track interactions of different originally primary particles.This approach allows to directly set the number N of primary particles whose secondary particles can interact with each other.For details of this method, we refer to Derksen et al (2023).In the following, this method is referred to as the phsp method.

Method 2: new inter-track scorer in TOPAS-nBio 2.0
To simulate the yield of radicals following radiolysis with inter-track interactions, TOPAS-nBio version 2.0 includes a scorer called TsIRTInterTrack.Applying this scorer, the user can set the number of tracks with intertrack interactions, but also the relative position and time delay of the primary tracks.However, this method only works with the IRT method.An application of this scorer can be found in Ramos-Méndez et al (2020).This method will further be referred to as the TsIRTInterTrack method.

Simulation setup
An isotropic electron source of 4.5 keV positioned in the center of the coordinate system was used as the source in the simulations.This energy was chosen since electrons with this kinetic energy are important considering the formation of chemical spurs and DNA damage segments (Ward 1988).
Since for the phsp method a phase space file needs to be scored, which will be post-processed and used in the following simulations with inter-track interactions as explained in section 2.2.1, the same phase space file was used as a basis for the source in the simulations applying both methods to ensure consistent simulation setups.While using the TsIRTInterTrack method, the phase space file was not modified in its fundamental way, but using the phsp method, the phase space file had to be modified before simulating G-values with inter-track interactions.This phase space file is scored on the surface of a water-filled7 sphere with a radius of 0.5 nm placed around the source.
Using both methods, the G-value was scored in a sphere with 5 μm radius filled with water.For the phsp method, the scorer TsIRTGvalue was applied to score the G-value of each molecule type making use of the IRT chemistry method.Here, the IRT chemistry model is employed for simulating the chemical stage of the simulation since the subsequent TsIRTInterTrack method solely works with the IRT chemistry approach, and consistent simulation setups were required for comparison of the two methods enabling inter-track interactions.For the TsIRTInterTrack method, the scorer TsIRTInterTrack was applied, which enables the simulation of inter-track interactions by specifying the number N of tracks with inter-track interactions as well using the IRT method.Using this scorer, no additional offset and no time delay was applied and all N primary particles were simulated at the same time.When applying the TsIRTInterTrack method, the deposited energy was also scored, since the TsIRTInterTrack scorer does not normalize the number of chemicals to 100 eV of deposited energy, which is necessary when calculating G-values.In general, the objective of both scoring methods is to provide a G-value.However, one scorer, the TsIRTInterTrack scorer, incorporates the generation of inter-track interactions, and the other, the TsIRTGvalue scorer, does not, as inter-track interactions are achieved through a different mechanism.
To vary the number of inter-track interactions, the number N of primary particles, from which the chemical species can react in inter-track interactions with each other, was set to N = 2, 3, 4, 5, 6, 10, 12, 15, 20, 30, 60.When using the phsp method, 60 primary particles were simulated in total in one simulation run for all N.For the phsp method, this is valid since all N are divisors of 60, and the phase space file is modified so that N original primary particles are assigned to one primary particle in the phase space file.Contrary, applying the TsIRTInterTrack method, exactly N primary particles need to be simulated in one simulation run.For statistical reasons, the simulations were repeated 300 times with different random seeds.
Physical interactions were simulated by using the GEANT4 physics constructor G4EmDNAPhysics_option2, which is described in detail by Incerti et al (2018) and has previously been investigated using TOPAS-nBio by Derksen et al (2021).For simulating chemical kinetics and radiolytic species, we used the TsEmDNAChemistry module with the corresponding list of chemical specifications given in TOPASDefaultReaction.txt. 8The end of the chemical stage was set to 1 μs.

Results
In figures 1 and 2, the G-values for various chemical species and the sum of all species (G all ) are shown in dependence of the number N of tracks with inter-track interactions for both the phsp and TsIRTInterTrack methods.In addition, the yield of recombined water molecules, formed through a reaction of • OH and H • or H 3 O and OH − , is shown as it contributes to G all .The general development of the G-value with increasing N agrees for both methods.With increasing N, the overall G-value decreases when using the phsp method as well as using the TsIRTInterTrack method.While the G-value decreases for most molecule types ( • OH, H 3 O, OH − , H 2 O 2 and e aq9 ) with N, for some molecules (recombined H 2 O, H 2 ) it increases.In comparison to our previous work (Derksen et al 2023) the trends of individual molecules with N in this study differ slightly because on the one hand, a different source with a different LET was used and on the other hand, a different list of chemical reactions was included and the SBS approach was used instead of the IRT approach in this study.For a detailed analysis of G-values in dependence of N, we refer to Derksen et al (2023).Regarding the summarized G-value of all molecule types, G all , the G-values agree very well with a maximum difference of 0.2%.Nonetheless, examining solely G all is insufficient, as the compensation of consumed and generated radicals may not expose differences between the methods.Consequently, a comparison of the G-values of individual molecules is essential.Regarding the individual radical species, the maximum difference of the G-value between both methods is observed for H 2 O 2 (3.9%) and the minimum difference is seen for OH − , for which the mean difference is less than 0.5%.Except for • OH and recombined water, the majority of G-values using the phsp method are smaller than those obtained using the TsIRTInterTrack method.This, and in combination with the fact that the amount of recombined water constitutes approximately 20%-35% of G all , proves our initial hypothesis that the positive difference of the G-value of • OH and recombined water using the TsIRTInterTrack method compared to the phsp method, can compensate the negative difference of the other chemical species.This way, the differences between phsp and TsIRTInterTrack for G all are very low and present an overstated positive agreement.Since the deviations between the methods are always either negative or positive for each molecule type, we suspect that these are systematic deviations.
In fact, it may initially seem confusing why these variations arise at all, considering that the simulations employ the same toolkit version, identical physics and chemistry lists, as well as the same dissociation tables, reactions, and diffusion constantsessentially under identical conditions.The sole distinction lies in the generation of inter-track interactions, implying that there are divergent operations in the code attributed to the generation of inter-track interaction affecting the chemical stage of the simulation.We investigated two potential factors that could be contributing to the deviations in the number of molecules: First, we examined the number of inelastic processes from which the chemical species are generated.Since inter-track interactions should not impact the physical stage, there should be no difference between the two methods for generating inter-track interactions.Nevertheless, we conducted this examination excluding any issues in the physical stage that might affect the yield of chemicals.Additionally, we investigated the number of chemical species processed directly at the beginning of the chemical stage.This analysis was carried out because in the pre-chemical stage, which connects the physical and chemical stages, not only the production of chemicals occurs, but also chemical reactions, albeit relatively infrequently.This could be handled differently between the two methods generating inter-track interactions resulting in a variation of the G-value between the methods.
In figure 3 (left), the number of inelastic processes normalized to 100 eV of deposited energy is shown for both methods generating inter-track interactions.Here, only ionizations, electronic excitation and electron attachment are considered for the inelastic processes, since these processes generate the chemical species that are processed during the chemical stage.For all N, differences between the number of inelastic processes of both methods are smaller than 0.1% and with an average difference of 0.01(4)%, there is no significant difference between both methods noticeable.
In figure 3 (right), the number of chemical species that are initially processed within the chemical IRT method normalized to 100 eV of deposited energy is shown for both methods, as well as the relative difference of the TsIRTInterTrack method to the phsp method.Using the IRT chemistry method, this number is given by default in the output of each simulation for every history and considers only the chemical species processing at the beginning of the chemical stage, not the pre-chemical stage.Interestingly, the number of chemical species processing is on average 1.13(4)% higher with the phsp method compared to the TsIRTInterTrack method.As this investigated number of chemical species is the same as the G-value scored at 1 ps, we further studied timedependent G-values comparing both methods for generating inter-track interactions.In figure 4, timedependent G-values are shown for N = 2 for G all (top), for • OH (bottom left) and for H 2 O 2 (bottom right) exemplary.As it can be seen in the bottom panel of each graph, relative differences of the TsIRTInterTrack method to the phsp method are always the highest at the beginning of the chemical stage and become smaller towards the end of the chemical simulation at 1 μs.Different G-values at the beginning of the chemical stage could be explained by the pre-chemical stage being treated differently in terms of chemical reactions between the two methods.
In combination with the results of the number of inelastic processes, we conclude that the differences in the G-value of some molecule types are caused by processing the pre-chemical stage differently between the phsp and TsIRTInterTrack method, as the G-values at 1 ps differ the most.Since detailed processes in the prechemical stage cannot be studied with TOPAS-nBio, and there is no reference for these data, it is not possible to decide which method provides the most accurate results.Nevertheless, since the differences are acceptable with regard to deviations of experimental and simulated data of radiolytic species (Ramos-Méndez et al 2018), both methods seem to be suitable for the simulation of inter-track interactions.
However, we noticed some differences in the handling of the two methods.Using the TsIRTInterTrack method, the number of primary particles set in one simulation run is defined by N, which means that exactly the number of tracks with inter-track interactions needs to be simulated as primary particles to obtain accurate Gvalues.Otherwise, invalid results will be obtained, which is further described in the appendix.Thus, a higher number of simulation runs is required to receive precise statistical results.In contrast, using the phsp method, multiple of N primary particles can be simulated in one simulation run, reducing the total number of runs required to obtain comparable statistical results.However, applying the phsp method, two separate simulations are necessary to create inter-track interaction because a phase space file, which is further modified, needs to be scored in a first simulation.The TsIRTInterTrack method has the advantage that, in principle, only one simulation is sufficient to model inter-track interactions since a predefined scorer can be used.For the phsp method, some programming effort is needed when modifying the phase space file.Nevertheless, the user knows exactly how inter-track interactions are generated since the simulations are manipulated accordingly.Another benefit of this method is that the G-value is scored via the TsIRTGvalue scorer, so no post-processing is necessary.In contrast, with the TsIRTInterTrack scorer, only the number of molecules, not normalized to the deposited energy, is scored, and therefore the deposited energy must also be scored in order to normalize the number of molecules to 100 eV of deposited energy to obtain the G-value.Furthermore, it is crucial to acknowledge that the phsp method is compatible with both the IRT and the SBS chemistry approach, while the TsIRTInterTrack method is exclusively compatible with the IRT method.However, it is worth noting that transitioning from the SBS to the IRT approach can result in differences of approximately 5% (Ramos-Méndez et al 2020).

Conclusion
G-values for different numbers of tracks with inter-track interactions were calculated in TOPAS-nBio and compared using two different and independently developed methods for generating inter-track interactions of N primary tracks.It was shown that the G-values of both methods agree well, and differences regarding individual molecule types are, on average, 1.5%.A reason for this could be a different consideration of the prechemical stage influencing the dynamics of the following chemical stage.Nevertheless, within the limits of uncertainties of experimental data of radiolysis products and deviations to simulated data at conventional dose rates, both methods seem suitable for simulating inter-track interactions, as they provide comparable G-values, even though both techniques were developed independently of each other.Furthermore, this comparison proved particularly informative by identifying a limitation in the TsIRTInterTrack method, as this method can only accommodate a maximum of N primary particles in one simulation.

Appendix. Additional information simulating inter-track interactions applying the TsIRTInterTrack scorer
Figure A1 shows that simulating more primary particles than the defined number N of tracks with inter-track interactions leads to invalid G-values using the TsIRTInterTrack scorer.In figure A1(left), G all is shown for enabling inter-track interactions of N = 2 primary tracks in dependence of the total number of simulated primary particles in one simulation run.In comparison to simulations with two primary particles, which equals N, in simulations using a higher number of primary particles than N, G all decreases significantly.In figure A1(right), G all in dependence of time is illustrated simulating 2, 10 and 50 primary particles in one simulation run with inter-track interactions of N = 2 tracks.As it can be seen, differences in the G-value between the number of simulating primary particles become mainly visible at the end of the chemical stage.Figure A2 shows an extract of the TOPAS simulation output using 10 primary particles for inter-track interactions of N = 2 tracks.Here, the number of simulating chemicals at the beginning of the chemical stage are shown after each pair of primary tracks (indicated by a new history) with inter-track interactions.Interestingly, the number of processing species does not remain roughly the same for each pair of simulated histories as expected, but increases significantly with further simulated histories.It seems like a pair of simulated histories considers also every chemical species created from all previous simulated histories, which would explain the increase of chemical species with increasing number of successive histories.This in turn explains the decreasing G-value for an increasing number of primary particles simulated for the same number of tracks with inter-track interactions shown in figure A1.This is due to the fact that the averaged G-value is actually calculated from G-values considering inter-track interactions of an increasing number of primary particle tracks (multiple of N, see figure A2) even though the number of tracks with inter-track interactions is defined in the scorer.Therefore, when using the TsIRTINterTrack scorer, care should be taken to ensure that the number of primary particles is equal to the number of tracks with inter-track interactions.
Figure A1.Left: G-value G all using the TsIRTInterTrack method simulating inter-track interactions of two tracks (N = 2) in dependence of the total number of primary particles simulated in one run.Right: Time-dependent G all for N = 2 simulating 2, 10 and 50 primary particles in one simulation run.For both graphs, statistical uncertainties are smaller than the symbol size and hence not depicted.

Figure A2
. Extract of the output simulating 10 primary particles with inter-track interactions of N = 2 primary tracks.Here, a history refers to a new primary particle.

Figure 1 .
Figure 1.Comparison of the G-value G all considering all molecules (upper left), • OH (upper right), recombined H 2 O (bottom left) and H 3 O (bottom right) using the phsp and TsIRTInterTrack method in dependence of the number N of tracks with inter-track interactions.In each upper panel, the G-value using both methods is depicted (phsp: blue triangles, TsIRTInterTrack: red diamonds) and in each bottom panel, the relative deviation of the TsIRTInterTrack results to phsp results are shown.Statistical uncertainties are represented by vertical error bars and correspond to one standard deviation.

Figure 2 .
Figure 2. Comparison of the G-value considering OH − (upper left), H 2 (upper right), H 2 O 2 (bottom left) and e aq (bottom right) using the phsp and TsIRTInterTrack method in dependence of the number N of tracks with inter-track interactions.In each upper panel, the G-value using both methods is depicted (phsp: blue triangles, TsIRTInterTrack: red diamonds) and in each bottom panel the relative deviation of the TsIRTInterTrack results to phsp results are shown.Statistical uncertainties are represented by vertical error bars and correspond to one standard deviation.

Figure 3 .
Figure3.Left: number of inelastic processes creating chemical species normalized to 100 eV of deposited energy in dependence of N shown for both methods (phsp: blue triangles, TsIRTInterTrack: red diamonds).Right: number of chemical species processing normalized to 100 eV of deposited energy in dependence of the number N of tracks with inter-track interactions using the phsp (blue triangles) and TsIRTIntertrack (red diamonds) method.Each bottom panel shows the relative difference of TsIRTIntertrack to phsp results.Statistical uncertainties are represented by vertical error bars and correspond to one standard deviation.

Figure 4 .
Figure 4. Comparison of time-dependent G-values for N = 2 using the phsp and TsIRTInterTrack method (top: G all , bottom left: Gvalues of • OH, bottom right: G-values of H 2 O 2 ).In each upper panel, the G-values are depicted for both methods and in each bottom panel, the difference of TsIRTInterTrack results relative to phsp results are shown.Statistical uncertainties represented by vertical error bars correspond to one standard deviation.