# Assessment of DNA damage with an adapted independent reaction time approach implemented in Geant4-DNA for the simulation of diffusion-controlled reactions between radio-induced reactive species and a chromatin fiber

Hoang Ngoc Tran, José Ramos-méndez, Wook-geun Shin, Yann Perrot, Bruce Faddegon, Shogo Okada, Mathieu Karamitros, Marie Davídková, Václav Štěpán, Sébastien Incerti, et al.

► To cite this version:

Hoang Ngoc Tran, José Ramos-méndez, Wook-geun Shin, Yann Perrot, Bruce Faddegon, et al. Assessment of DNA damage with an adapted independent reaction time approach implemented in Geant4-DNA for the simulation of diffusion-controlled reactions between radio-induced reactive species and a chromatin fiber. Medical Physics: The international journal of medical physics research and practice, 2021, 48 (2), pp.890-901. 10.1002/mp.14612. hal-03149166

HAL Id: hal-03149166

https://hal.science/hal-03149166v1

Submitted on 9 Mar 2021

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# Assessment of DNA damage with an adapted independent reaction time approach implemented in Geant4-DNA for the simulation of diffusion-controlled reactions between radio-induced reactive species and a chromatin fiber

**Hoang Ngoc Tran<sup>a</sup>**

IRSN, Institut de Radioprotection et de Sûreté Nucléaire, BP17, 92262 Fontenay aux Roses France

**José Ramos-Méndez**

Department of Radiation Oncology, University of California San Francisco, San Francisco, CA 94115, USA

**Wook-Geun Shin**

Université de Bordeaux, CNRS/IN2P3, UMR5797, Centre d'Études Nucléaires de Bordeaux Gradignan, 33175 Gradignan, France
Department of Radiation Convergence Engineering, Yonsei University, 26493 Wonju, Korea

**Yann Perrot**

IRSN, Institut de Radioprotection et de Sûreté Nucléaire, BP17, 92262 Fontenay aux Roses France

**Bruce Faddegon**

Department of Radiation Oncology, University of California San Francisco, San Francisco, CA 94115, USA

**Shogo Okada**

KEK, 1-1, Oho, Tsukuba, Ibaraki 305-0801, Japan

**Mathieu Karamitros**

Radiation Laboratory, University of Notre Dame, Notre Dame, In 46556, USA

**Marie Davídková and Václav Štěpán**

Department of Radiation Dosimetry, Nuclear Physics Institute of the CAS, Prague, Czech Republic

**Sébastien Incerti**

Université de Bordeaux, CNRS/IN2P3, UMR5797, Centre d'Études Nucléaires de Bordeaux Gradignan, 33175 Gradignan, France

**Carmen Villagrasa**

IRSN, Institut de Radioprotection et de Sûreté Nucléaire, BP17, 92262 Fontenay aux Roses France

(Received 8 July 2020; revised 6 November 2020; accepted for publication 10 November 2020; published 20 December 2020)

**Purpose:** Simulation of indirect damage originating from the attack of free radical species produced by ionizing radiation on biological molecules based on the independent pair approximation is investigated in this work. In addition, a new approach, relying on the independent pair approximation that is at the origin of the independent reaction time (IRT) method, is proposed in the chemical stage of Geant4-DNA.

**Methods:** This new approach has been designed to respect the current Geant4-DNA chemistry framework while proposing a variant IRT method. Based on the synchronous algorithm, this implementation allows us to access the information concerning the position of radicals and may make it more convenient for biological damage simulations. Estimates of the evolution of free species as well as biological hits in a segment of DNA chromatin fiber in Geant4-DNA were compared for the dynamic time step approach of the step-by-step (SBS) method, currently used in Geant4-DNA, and this newly implemented IRT.

**Results:** Results show a gain in computation time of a factor of 30 for high LET particle tracks with a better than 10% agreement on the number of DNA hits between the value obtained with the IRT method as implemented in this work and the SBS method currently available in Geant4-DNA.

**Conclusion:** Offering in Geant4-DNA more efficient methods for the chemical step based on the IRT method is a task in progress. For the calculation of biological damage, information on the position of chemical species is a crucial point. This can be achieved using the method presented in this paper. © 2020 The Authors. Medical Physics published by Wiley Periodicals LLC on behalf of American Association of Physicists in Medicine [https://doi.org/10.1002/mp.14612]

**Key words:** diffusion-controlled reaction, DNA damage, Geant4-DNA, IRT
# 1. INTRODUCTION

Among the possible causes of the effects induced by ionizing radiation on living organisms, DNA damage is of particular interest¹ and many techniques are being developed to both measure and predict it. The creation of radiolytic species through ionization and excitation in the vicinity of the biological target significantly contributes via the so-called indirect damage. Investigations of indirect damage consider the mechanisms of free radical attack on the DNA target, including taking account of its chemical and geometrical structure. The accuracy of the assessment of indirect damage depends in part on the knowledge of the positions of radiolytic species in relation to the DNA molecule. This depends on the ability of each of these species to diffuse and react with other molecules in the environment, including DNA. For modeling purposes, the basic DNA elements (2-deoxyribose, phosphate, base) are in Geant4 often represented by static reactive spheres (or sinks) linked together to form more complex structures, from the nucleotide pair to the genome.²⁻⁶ The reaction mechanism is diffusion-controlled, as it is supposed to be,⁷ and triggered when radiolytic species diffuse and encounter a reactive site (generally reactive spheres) that are either representing the sugar-phosphate backbone or the bases of a nucleotide. This simple model still ignores the overlap of the reaction sites as well as the overlap of multiple reactive centers in the sugar-base system.⁸ However, it is well suited when focusing on the assessment of the DNA damage on nucleotide or base-pair level.

The theory of diffusion-controlled reactions⁹ describes solutions of the diffusion equation (Smoluchowski equation) with boundary conditions. Based on these solutions, stochastic simulation techniques, through step-by-step (SBS) or independent time reaction (IRT) methods, describe that the diffusion of molecules and the reactions between reactants¹⁰ are used to simulate the evolution of heterogeneous reactive species distributions from the initial radiolysis. The SBS method¹¹,¹² is able to provide the spatial positions of the diffusing species at a given time. This knowledge of the temporal evolution of trajectories calculated in discretized time steps may be advantageous for the assessment of indirect DNA damage, as discussed in Section 3. The accuracy of this method depends on the determination of the time step, requiring a compromise between accuracy, which makes necessary a sufficiently small time step, and calculation time, which increases dramatically with a decreasing time step. Indeed, the SBS method shows very long computation times for radiobiology applications. In that frame, IRT method has been developed to save a considerable amount of computing time. However, this method does not provide exact positions of the radicals in time.¹⁰ Indeed, the IRT method is based on the “Independent Pair Approximation”; thus, reactive pairs are assumed independent, that is, the reaction time between two reactants does not depend on the other reactants present in the medium, and the diffusion time can be calculated directly using Green’s functions.¹³,¹⁴ The predictions of the IRT method agree with SBS for systems of two particles in

diffusive motion for totally and partially diffusion-controlled reactions.¹⁰,¹⁵

An important question is how the “Independent Pair Approximation” method could be applied to reactions of radicals with biological molecules which are considered static in the simulated environment. Based on a simple model of biological system as DNA elements, which are represented by two reactive spheres linked together at a fixed distance, Bluett and Green⁸ have reported the competitive effect as a function of the intersphere distance of reaction probability of diffusion-influenced reactions and found an overestimation of this probability with the IRT method compared to an exact analytical solution. Due to the complexity of many reactive sites in a complex DNA geometry model, solving the problem as an exact solution is intractable. In this work, an IRT method is used in Monte Carlo simulation of free radical attacks on DNA.

Among the Monte Carlo codes, Geant4¹⁶⁻¹⁸ is a general purpose toolkit first dedicated to high energy and nuclear physics. Developments have enabled other applications in medical or space fields. Geant4-DNA¹⁹⁻²² is an extension of Geant4 for the modeling of biological damage induced by ionizing radiation at the nanometer scale. To complement comparisons of our results, we have used published data obtained with TOPAS-nBio.²³ TOPAS is a Monte Carlo²⁴ system that wraps and extends Geant4 to facilitate the use of Monte Carlo simulation of radiotherapy by medical physicists. TOPAS-nBio is the extension of TOPAS to model radiobiological at the cellular and subcellular scale. Since the release of Geant4 version 10.1, Geant4-DNA provides a computational framework for the physicochemical and chemical stage in which the SBS method, that includes Brownian motion and chemical reactions between molecules resulting from water radiolysis, is available for users.¹¹ Applications are mainly related to the validation of radiolytic yields and the evaluation of radiation-induced DNA damage. To enhance the potential of Geant4-DNA, developments are ongoing in different groups with the goal to provide IRT methods to Geant4-DNA users.²⁵,²⁶ In particular, an implementation of the original IRT method has been developed and will be released soon.²⁵ In this work, an implementation of an IRT variant is proposed. This method is implemented in the current chemistry module of Geant4-DNA with the idea of keeping the spatiotemporal information of the radiolytic species to calculate DNA damage. The proposed method brings the possibility to choose a reasonable compromise between the gain in calculation time over the SBS method and the accuracy of the assessment of the DNA damage using a molecular geometrical model of a segment of chromatin fiber. To do so, this new IRT method is complemented with the synchronous algorithm described in Karamitros et al.¹¹ Prior to the public release of this method, this implementation is compared to the current SBS technique, and the differences are assessed in terms of calculation time and quantities of interest. The radiochemical yields calculated with our IRT method are compared to those obtained with the SBS method of Geant4-DNA, to those obtained with TOPAS-nBio, using the IRT method described in Schuemann et al.,²³ and to
published experimental data. Finally, a DNA reaction model of free radicals based on reactive sinks of a heterochromatin fiber geometry is used to compare DNA damage calculated in Geant4-DNA using the newly implemented IRT method and the existing SBS method.

## 2. MATERIALS AND METHODS

### 2.A. Diffusion-controlled reactions

The reactions of radicals with each other or with DNA molecules can be described as diffusion-controlled reactions.²⁷ Diffusion-controlled reactions are reactions occurring immediately when molecules encounter each other. The reaction rate is equal to the rate of diffusion of the reactants through the medium when the reaction is totally diffusion-controlled and takes into account the steady-state rate constant from the encounter to their reaction when the reaction is partially diffusion-controlled. In this section, we briefly present the probability distribution of these reactions and sampling techniques for an independent two-body system.

#### 2.A.1. Totally diffusion-controlled reaction

Suppose that the two reactants considered are labeled 1 and 2. The solution of the Smoluchowski equation with an absorption boundary condition (so-called totally diffusion-controlled reaction) transformed to radial Green's function in a spherical coordinate system is used to deduce the probability of reaction $p(t|r_0)$ as Ref. [13]:

$$p(t|r_0) = \frac{R}{r_0} \text{erfc} \left[ \frac{r_0 - R}{\sqrt{4Dt}} \right] \quad (1)$$

where R is the reaction radius which is sum of reaction radius of both reactants 1 and 2, $r_0$ is the initial distance between reactants 1 and 2, D is the relative diffusion coefficient and equals $D_1 + D_2$ where $D_1$ and $D_2$ are the diffusion coefficients of reactant 1 and 2, respectively. The experimental reaction rate $k_{obs}$ is used to define R by the relation $k_{obs} = 4\pi RD$.

This probability can be used for charged species by replacing effective distances $d_{eff}$ to $r_0$:

$$d_{eff} = - \frac{r_c}{1 - e^{(r_c/R)}}$$

where $r_c$ is the Onsager radius.¹³ In this work, $r_c = 0.71$ nm (particles with charge equal to 1, in liquid water at 25°C).

#### 2.A.2. Partially diffusion-controlled reaction

In the radiation boundary condition (or partially diffusion-controlled reaction), Smoluchowski equation gives a well-known distribution for calculating the distribution probability:

$$p(t|r_0) = P_0 \left( \text{erfc}(y) - e^{(x^2+2xy)} \text{erfc}(x+y) \right) \quad (2)$$

where $P_0 = \frac{R_{eff}}{r_0}$, $R_{eff} = R \frac{k_{act}}{k_{dif}+k_{act}}$, $x = \frac{(k_{dif}+k_{act})\sqrt{Dt}}{k_{dif}R}$ and $y = \frac{r-R}{\sqrt{4Dt}}$, $k_{dif} = 4\pi RD$ is the reaction rate constant of transport “to encounter” and $k_{act} = 4\pi R^2 v$ is the steady-state rate constant “from encounter” to the reaction which can be defined by the velocity of reaction $v$ but usually calculated using the experimental reaction rate $k_{obs}$ in relation:

$$\frac{1}{k_{obs}} = \frac{1}{k_{dif}} + \frac{1}{k_{act}} \quad (3)$$

For charged species, x and y can be calculated as follows¹³:

$$x = \frac{4R^2\alpha'}{r_c^2} \sqrt{A} \frac{t}{D} \sinh^2\left(\frac{r_c}{2R}\right), \quad (4)$$

and

$$y = \frac{r_c \left[ \coth\left(\frac{r_c}{2r_0}\right) - \coth\left(\frac{r_c}{2R}\right) \right]}{2\sqrt{4Dt}}, \quad (5)$$

where

$$\alpha' = v + \frac{r_c D}{R^2 \left[ 1 - e^{(-r_c/R)} \right]}, \quad (6)$$

In this case,

$$R_{eff} = \frac{-r_c}{1 - \left[ 1 + D r_c / v R^2 \right] e^{(r_c / R)}}, \quad (7)$$

but practically is calculated from the experimental reaction rate $k_{obs}$ with the relation $k_{obs} = 4\pi R_{eff} D$.

**SBS method:** One of the methods that can be implemented for the simulation of diffusion-controlled reactions is SBS method. A detailed description of the SBS model using dynamic time step implemented in Geant4-DNA can be found in Karamitros et al.¹¹ Briefly, this method proposes a time step model that allows the choice of virtual time steps during which the reaction cannot occur with at least 95% (by default) confidence (named dynamic time step). One can visualize this as creating a protection domain that surrounds the particle, ensuring that this particle will not react with any other particle up to its border with 95% confidence. Therefore, in the protection domain, the particle is considered approximately independent and may be able to take longer diffusion steps. The function in Eq. (1) is used to evaluate this time step. This process is repeated many times until a chemical interaction takes place. Thus, we may have one or many time steps before the reaction occurs. To avoid the scenario of many small time steps, the Minimum Time Steps and the Brownian bridge technique have been added to limit the number of time steps to an encounter. While Minimum Time Steps constrain the minimum time step allowed for each reactant pair, the Brownian bridge technique computes the probability of encounter during their Minimum Time Steps and thus compensates for the “missed” reactions. Detail of these techniques can be found elsewhere.¹⁰,¹¹
To save computation time, it can be decided to increase the Minimum Time Step value as function of the virtual time used in the simulation using the predefined static time steps proposed by Kreipl et al.,²⁸ since radical spurs become more and more sparse. This is the SBS-dynamic time steps using variable Minimum Time Step values. It is also important to note that the current SBS-dynamic time step method assumes that all reactions are totally diffusion-controlled.

**IRT method:** IRT is another method for simulating diffusion-controlled reactions. By simplifying the multiple particle problem to the two-particle problem in an approximation, the IRT method is based on the relative comparison of random times to reaction calculated for all possible reactant pairs independently of the particle system. The minimum reaction time obtained is selected such that the corresponding reaction will occur in the next step. These random times are sampled from the reaction probability distributions of the reactant pairs which are considered independent of each other and in infinite space. Therefore, only the initial spatial distribution of reactant pairs plays an important role on these probability distributions.

The IRT method determines the minimum time to the next reaction. Reactive products created by reactions that have occurred can undergo reactions with other reactants. These new reactions need to be considered and included in the possible reaction times, depending on their initial positions. Based on an approximation proposed by Clifford et al.,¹⁰ models for reactive products allow to calculate their reaction times. We briefly summarize these models as follows:

* **Diffusion and time approaches:** These approaches propose the deduction of the reaction time of reactive products based on their parent pair position and the difference of diffusion coefficient between them and their parents. While the time approach uses a rescaling of the reaction time of parent reactant pair to calculate the reaction time of products, the diffusion approach recalculates the random time corresponding to the distance between products and their reactants. Thus, the position of reactants is not determined.

* **Position approach:** In this approach, the minimum required time for one possible reaction in the system is considered as a “time step.” In the meantime, all remaining reactive radicals are diffused in a single step to this time step following a Gaussian distribution of their positions as Brownian objects. In this approach, the position of reactive products following the reactions are determined explicitly and the reaction times between the products and all the remaining radicals involving the products are simulated.

## 2.B. IRT method proposed in this work

### 2.B.1. Synchronous event-driven algorithm

As explained in the introduction, in the traditional IRT method,¹⁰,¹³ due to the fact that space is assumed infinite, only relative distances of these reactive pairs are considered

to calculate their reaction probability. Under this assumption, IRT is an exact method for two diffusive particle systems in an infinite and homogeneous space.⁸ Based on this assumption, the reaction pairs are simulated independently and asynchronously in time and space. While this is a considerable advantage in terms of computing time, the spatial-temporal information of the system is not simulated explicitly. As a complementary approach, our particular implementation of the IRT method is designed to be adapted to the chemical module of Geant4-DNA. This implementation uses the sampled random time given by the reaction probabilities in Section 2.A [by using Eqs. (1) and (2)] as a function of time step (or time slice) to the next reaction that should occur. In other words, instead of optimizing the time step to the next reaction as Geant4-DNA SBS-dynamic time step method does, we calculate directly the time to the reaction of an independent reactive pair. Therefore, each step of the simulation consists of two stages: (a). the reaction times of all possible reactant pairs are calculated; (b). the minimum reaction time and corresponding reaction are selected to trigger. The reactive product positions created in the reactions and the remaining molecules are considered explicitly together to diffuse for the time step using the position approach (see Section 2.A.2). Then, based on their new positions, the new random reaction times are re-evaluated sequentially for all the radicals in the system and the new minimum reaction time and corresponding reaction is selected for next time step. This procedure is repeated until the end time of simulation. Consequently, for each time step, we need to update the time and position of all molecules in diffusion. This is the main drawback of this approach. However, this provides explicitly spatiotemporal information of the reactive species after each time step which can then be coupled with information on the geometrical boundaries or the biological target.

When reactions with DNA or other biomolecules are considered, the reaction times are compared for each reactive pair, keeping the spirit of independent pair approximation. The resulting minimum reaction time is then compared with the minimum reaction time between radicals themselves. The smaller of the times will define the type of reaction that will occur: radical/radical (Table I) or radical/DNA (Table III).

It is worth noting that due to the free diffusion of all remaining radicals for each time step, the inter-reactant distance of some reactant pairs may be smaller than the reaction radius and they may cause extra reactions with the probability of reaction²⁹:

$$ \phi = \frac{e^{-\left(\frac{r_c}{R}\right)} - e^{-\left(\frac{r_c}{R+r_s}\right)}}{e^{-\left(\frac{r_c}{R}\right)} - e^{-\left(\frac{r_c}{R+r_s}\right)} - \left(\frac{k_{dif}}{k_{act}}\right)\left[1 - e^{-\left(\frac{r_c}{R}\right)}\right]}, \quad (8) $$

For radicals without charge, $\phi$ is:

$$ \phi = \frac{R_s}{R_s + \left(\frac{k_{dif}}{k_{act}}\right)\left(R + R_s\right)}, \quad (9) $$

The interpretation may be addressed by an unsuccessful encounter where $R_s$ is the separation distance of the encounter, equals 0.3 nm.
<table><caption>TABLE I. Reactions and reaction rate coefficients used in this work. For IRT, the list of reactions are assigned by type: partially diffusion-controlled, totally diffusion-controlled, and spin statistical factor.³⁰ For SBS, all the reactions in the table are considered totally diffusion-controlled reactions.</caption><thead><tr><th>Reaction</th><th>k<sub>obs</sub> (×10<sup>10</sup> M<sup>-1</sup>s<sup>-1</sup>)</th><th>Partially diffusion-controlled</th><th>Totally diffusion-controlled</th><th>Spin statistical factor</th></tr></thead><tbody><tr><td>H<sup>*</sup> + e<sup>-</sup><sub>aq</sub> + H<sub>2</sub>O → OH<sup>-</sup> + H<sub>2</sub></td><td>2.5</td><td></td><td>X</td><td>X</td></tr><tr><td>H<sup>*</sup> + <sup>*</sup>OH → H<sub>2</sub>O</td><td>1.55</td><td>X</td><td></td><td></td></tr><tr><td>H<sup>*</sup> + H<sup>*</sup> → H<sub>2</sub></td><td>0.503</td><td></td><td>X</td><td>X</td></tr><tr><td>H<sub>2</sub>O<sub>2</sub> + e<sup>-</sup><sub>aq</sub> → OH<sup>-</sup> + <sup>*</sup>OH</td><td>1.1</td><td>X</td><td></td><td></td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + e<sup>-</sup><sub>aq</sub> → H<sup>*</sup> + H<sub>2</sub>O</td><td>2.11</td><td>X</td><td></td><td></td></tr><tr><td>H<sub>3</sub>O<sup>+</sup> + OH<sup>-</sup> → 2H<sub>2</sub>O</td><td>11.3</td><td></td><td>X</td><td></td></tr><tr><td><sup>*</sup>OH + e<sup>-</sup><sub>aq</sub> → OH<sup>-</sup></td><td>2.95</td><td>X</td><td></td><td></td></tr><tr><td><sup>*</sup>OH + <sup>*</sup>OH → H<sub>2</sub>O<sub>2</sub></td><td>0.55</td><td>X</td><td></td><td></td></tr><tr><td>e<sup>-</sup><sub>aq</sub> + e<sup>-</sup><sub>aq</sub> + 2H<sub>2</sub>O → 2OH<sup>-</sup> + H<sub>2</sub></td><td>0.636</td><td></td><td>X</td><td>X</td></tr></tbody></table>

For totally diffusion-controlled reactions for which ($k_{act} \to \infty$), $\phi=1$, the reaction immediately occurs. Moreover, a spin statistical factor for any encounter of $H^*$ and $e_{aq}^-$ or their combinations (see Table I) is added to consider the fact that only the singlet configuration of their spins allows the occurrence of the reaction.¹⁴

### 2.B.2. Sampling method for partially diffusion-controlled reactions

Following the description of Green et al.,¹³ the IRT method requires the generation of random times from these reaction probability distributions (see Section 2.A) for each of the reactant pairs at any initial distance between them. For total diffusion-controlled reactions, the inversion method¹³ provides a simple way of calculating these reaction times by using the inverse of the error function in Eq. (1). The same method is difficult to apply when considering partially diffusion-controlled reactions [Eq. (2)], because of the complexity of this equation. In this case, to sample time $t$, a technique previously proposed by Bluett and Green⁸ can be applied, in which an imaginary trajectory of one reactant (particle 1 which is considered moving) from the initial position to the reaction point with other reactant (particle 2 which is considered not moving) is split into two parts of the time step. The first part of time step corresponds to the time spent from the initial position of reactant 1 to encounter at the reaction radius $R$. This time step can be sampled using the distribution of Eq. (1) as a totally diffusion-controlled reaction. After the encounter, the reactant is assumed to be standing at the boundary of the reaction radius where $r_0 = R$ and is substituted in Eq. (2):

$$p(t|R) = \frac{k_{act}}{k_{dif} + k_{act}} (1 - e^{x^2} \text{erfc}(x)) \quad (10)$$

The second part of the time step corresponds to the time between the first encounter to the final reaction between both reactants and is sampled from this distribution [Eq. (10)]. We summarized the second sampling as follows⁸:

1. Generate two uniform [0,1] random numbers $U_1$ and $U_2$

2. If $U_1$ is greater than $\frac{k_{act}}{k_{dif}+k_{act}}$, then the reaction will not take place

3. Otherwise, generate an absolute value of normally distributed random Y with mean 0 and standard deviation $\sqrt{2}$

4. Calculate $t = \left\{ \frac{1}{D} \left( \frac{Rk_{dif}\ln(U_2)}{(k_{dif}+k_{act})Y} \right)^2 \right\}$ as a random time for the second sampling.

The reaction time is the sum of times of two parts of the time step. Figure 1 presents the sampled times using the sampling algorithm for four reactions among those shown in Table I combining both totally and partially diffusion-controlled reactions. These four reactions are sampled using $10^6$ histories with an initial distance in this example being 1.0 or 1.5 nm for each reactant pair. In all cases, the sampling distributions are as predicted by the probability density functions that are mentioned by Ref. [30]. The error bars in the figure show the square of weight deviation for each point.

### 2.C. IRT and SBS reaction schemes

Two reaction schemes for IRT were studied in this work. The first scheme combined reactions and reaction rate coefficients for partially diffusion-controlled, totally diffusion-controlled and spin statistical factor correction (Table I) and indicated as “this work” in the results section. The second scheme only considered totally diffusion-controlled reactions, indicated as “this work-TDC” in the results section. This second scheme had been considered for a “fair” comparison with the SBS method, currently available in Geant4-DNA, where all the reactions shown in Table I are handled as totally diffusion-controlled.

### 2.D. Description of the geometrical models and simulation setups

The time evolution of the number of molecular species created or lost per 100 eV of deposited energy (G-value)
<figure><img src="image_2.png" /></figure>

<figure><img src="image_3.png" /></figure>

<figure><img src="image_4.png" /><figcaption>FIG. 1. Sampled times to reaction from the probability densities of reactions ("theory")$^{30}$ using the IRT method presented in this work. H$^{•}$ + •OH → H$_{2}$O (initial distance of reactants 1.0 nm, partially diffusion-controlled), H$_{3}$O$^{+}$ + e$^{-}_{aq}$ → H$^{•}$ + H$_{2}$O (initial distance of reactants 1.5 nm, partially diffusion-controlled), H$_{3}$O$^{+}$ + OH$^{-}$ → 2 H$_{2}$O (initial distance of reactants 1.5 nm, totally diffusion-controlled), H$^{•}$ + e$^{-}_{aq}$ + H$_{2}$O → OH$^{-}$ + H$_{2}$ (initial distance of reactants 1.5 nm, totally diffusion-controlled). [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

was obtained with the different reaction schemes of Sec-
tion 2.B. The simulation setup consisted of a box of homo-
geneous liquid water (1 g/cm³) with dimensions
representing a “pseudo infinite” region (cube of 1 km side).
For computational time and memory reasons, the particle
source consisted of 80 keV electrons shot from the center of
the box. The initial chemical species are obtained from
physical track interactions for each primary electron consid-
ered independently of each other and depositing a total
energy of 1 keV using the physics list

<figure><img src="image_5.png" /></figure>

<figure><img src="image_6.png" /></figure>

<figure><img src="image_7.png" /></figure>

G4EmDNAPhysics_option8 in order to facilitate compar-
isons with previously published results from TOPAS-nBio.²³
The corresponding process classes, model classes, and
energy ranges of this physics list for electrons and protons
are shown in Table II. It is equivalent to the physics list used
by default except for the elastic scattering of electrons where
the CPA100 code model is used for energies below
256 keV. More details can be found elsewhere. ³¹ The tem-
poral evolution of G-values is recorded over the time range
from 1 ps to 1 µs.
<table><caption>TABLE II. Content of the physics list G4EmDNAPhysics_option8 used in this study: processes, models, and energy ranges for electrons and protons.</caption><thead><tr><th>Particle</th><th>Process</th><th>Model</th><th>Energy range</th></tr></thead><tbody><tr><td rowspan="6">Electron</td><td>G4DNAElastic</td><td>G4DNACPA100ElasticModel</td><td>11 eV–256 keV</td></tr><tr><td>G4DNAExcitation</td><td>G4DNAChampionElasticModel</td><td>256 keV–1 MeV</td></tr><tr><td>G4DNAIonisation</td><td>G4DNABornExcitationModel</td><td>9 eV–1 MeV</td></tr><tr><td>G4DNAVibExcitation</td><td>G4DNABornIonisationModel</td><td>11 eV–1 MeV</td></tr><tr><td>G4DNAAttachment</td><td>G4DNASancheExcitationModel</td><td>2 eV–100 eV</td></tr><tr><td>G4DNAIonisation</td><td>G4DNAMeltonAttachmentModel</td><td>4 eV–13 eV</td></tr><tr><td rowspan="5">Proton</td><td>G4DNAElastic</td><td>G4DNAIonElasticModel</td><td>100 eV–1 MeV</td></tr><tr><td>G4DNAExcitation</td><td>G4DNAMillerGreenExcitationModel</td><td>10 eV–500 keV</td></tr><tr><td>G4DNAIonisation</td><td>G4DNABornExcitationModel</td><td>500 keV–100 MeV</td></tr><tr><td>G4DNAChargeDecrease</td><td>G4DNAUddIonisationModel</td><td>0 eV–500 keV</td></tr><tr><td>G4DNAChargeDecrease</td><td>G4DNABornIonisationModel</td><td>500 keV–100 MeV</td></tr><tr><td></td><td></td><td>G4DNAingfelderChargeDecreaseModel</td><td>100 eV–100 MeV</td></tr></tbody></table>

In a practical way, in next releases we will modify the recent Geant4-DNA example dnadamage1 in order to offer the possibility to use both chemistry approaches: the current SBS method and the new IRT method presented in this work and this only by changing the chemistry list. More details about the C++ classes and their structure related to this method will be given to the users through the README file of the example.

In that example, as well as in this work, the DNA model was generated with the DNAFabric tool.⁶ DNAFabric is a stand-alone software which enables complex 3D geometries to be generated and visualized, in particular DNA structures from the nucleotide constituents (2-deoxyribose, phosphate, and DNA bases) to the whole genome representation. The geometrical model generated consisted of a piece of a 40 nm heterochromatin fiber including 3640 nucleotide pairs in a cubic voxel of 40 nm centered in a box of 2 µm side made of liquid water (1 g/cm³). These DNA elements representing the 2-deoxyribose, phosphate, and DNA bases (adenine, guanine, cytosine, and thymine) were considered as static spherical traps. Using Protein Data Bank information,³² atom's volumes and positions of each of these molecules were assembled in a unique sphere. From these spheres, nucleotide pairs were constituted and served as the basic units forming a double helix of the *B-DNA* configuration. More information about the generation of this geometrical model can be found elsewhere.²,⁶ Some segments of this double helix were then twisted around histone proteins represented by a sphere

acting like a histone with a 2.4-nm radius, building a nucleosome.³³ After that, nucleosomes were linked together to create the chromatin fiber (see Fig. 2). Additionally, a volume corresponding to water layers wrapped around each of the DNA elements was modeled as an outer hydration shell. This geometry can also integrate histones whose scavenging capacity would be simulated by introducing a specific reaction. Nevertheless, in this study, histones were not considered in order to maximize the probability of interaction and reaction with the DNA geometry.

In order to maximize this probability as well, the simulation was performed using 500 keV protons of normal incidence with respect to a side of the voxel and homogeneously distributed in a circle of radius 20 nm. To compute DNA damage, we used the criteria presented elsewhere.²,³³ In short, the total number of damages was computed as the sum of direct and indirect damage. A direct damage is scored if the cumulative deposited energy from ionizations and excitations in the individual volumes of a nucleotide backbone (i.e., the volumes representing a group of the phosphate, the 2-deoxyribose, and the hydration shell) is >17.5 eV for a given incident particle.³³ These volumes are assumed to be filled with liquid water for physical track interactions using the physics list G4EmDNAPhysics_option8, used for consistency with the calculation of G-values presented in this work.

For the DNA damage calculations, radicals undergo immediately the reactions (in Table III) with static DNA elements to produce indirect damage. We suppose that these

<figure><img src="image_8.png" /><figcaption>FIG. 2. 40 nm heterochromatin fiber (including in 18 nucleosomes and 19 linkers for a total of 3640 nucleotide pairs) in a cubic voxel of 40 nm generated with the DNAFabric tool.² [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>
<table><caption>TABLE III. Reaction rates between the $e_{aq}^{-}$, H*, *OH radicals and the DNA bases or 2-deoxyribose.⁷</caption><thead><tr><th>Reaction</th><th>Reaction rate (10<sup>9</sup> M<sup>-1</sup> s<sup>-1</sup>)</th></tr></thead><tbody><tr><td>2-deoxyribose + *OH</td><td>1.8</td></tr><tr><td>Adenine + *OH</td><td>6.1</td></tr><tr><td>Guanine + *OH</td><td>9.2</td></tr><tr><td>Thymine + *OH</td><td>6.4</td></tr><tr><td>Cytosine + *OH</td><td>6.1</td></tr><tr><td>2-deoxyribose + e<sup>-</sup><sub>aq</sub></td><td>0.01</td></tr><tr><td>Adenine + e<sup>-</sup><sub>aq</sub></td><td>9.0</td></tr><tr><td>Guanine + e<sup>-</sup><sub>aq</sub></td><td>14.0</td></tr><tr><td>Thymine + e<sup>-</sup><sub>aq</sub></td><td>18.0</td></tr><tr><td>Cytosine + e<sup>-</sup><sub>aq</sub></td><td>13.0</td></tr><tr><td>2-deoxyribose + H*</td><td>0.029</td></tr><tr><td>Adenine + H*</td><td>0.10</td></tr><tr><td>Thymine + H*</td><td>0.57</td></tr><tr><td>Cytosine + H*</td><td>0.092</td></tr></tbody></table>

reactions are totally diffusion-controlled. Thus, the random reaction times between radicals and DNA elements are calculated by Eq. (1). A reaction between a radical and a DNA element is counted as a primary damage event. Thereafter, the radical is killed and the damaged DNA element is no longer available for further reaction. Note that our damage are primary damage events that are not necessarily transformed into single-strand break (SSB) or double-strand break (DSB). Going out of the voxel volume, the radical species will be discarded. Table III shows the reactions used between the $e_{aq}^{-}$, H*, *OH radicals and the DNA bases or 2-deoxyribose and their reaction rates. For DNA damage calculations, we also perform a comparison with *SBS-dynamic time step using 0.1 ps Minimum Time Steps* for overall virtual simulation time. This is *SBS-dynamic time step using Minimum Time Steps of 0.1 ps*.

# 3. RESULTS AND DISCUSSION

## 3.A. Results on radiolytic yields

First, the reliability and performance of our IRT method implementation in Geant4-DNA are examined by using this method in the simulation to calculate time-dependent G-values of the chemical species created by water radiolysis, and comparing the results with other simulated data and experimental data from the literature. These quantities were calculated with our IRT method implementation (with both reaction schemes described in Section 2.B), the SBS-dynamic time step of the current public release of Geant4-DNA (Geant4 10.6), and the IRT method of TOPAS-nBio code (Topas-IRT).²³ Note that the physical stage and physicochemical stage³¹ are the same for these three approaches.

Figure 3 shows G-values as a function of time for different species (*OH, H₃O⁺, H*, e⁻ₐq, H₂, and H₂O₂) produced by 80 keV electrons depositing a total energy of 1 keV. Our

results follow the general trend of the experimental and calculated data. Our results are closer to Topas-IRT than to the SBS-dynamic time step results. Note that for both IRT simulations assigned, the same reaction types (see Section 2) and the reaction rates from Plante and Devroye³⁰ are used, whereas SBS-dynamic time step considers that all reactions are totally diffusion-controlled. The differences about a few percent between our results and Topas-IRT may be attributed to the re-estimation of the random reaction times for each time step made by our method.

For each selected random reaction time, besides the calculation of the reactive products with the position approach (see Section 2.A.2), the position of all the other chemical species is also recalculated with the synchronous algorithm (Karamitros et al.¹¹). In this way, both the evaluation of contact reactions and interaction with geometry can be performed (see Section 2.B.1). Nevertheless, the recalculated positions of chemical species forced to re-evaluate the random reaction times which may in turn lead to a higher reactivity in comparison with TOPAS-nBio, which implements the traditional IRT method. Our results show a difference of about 15% with SBS-dynamic time step G values at 1 µs for *OH. Note that G-value of H₂ shows an underestimation due to the inclusion of the spin effect. The balance material Eq. (11)¹² was also verified, that is, the number of reduced species agreed with the number oxidative species within 0.23%, in the entire time domain.

$$G(-H_2O) = G(e_{aq}^{-}) + 2G(H_2) + G(H) = G(OH) + 2G(H_2O_2). \qquad (11)$$

## 3.B. Results on damage yields

Damage yields presented here represent the addition of direct and indirect damage per particle to the DNA backbone and bases (see Table III). The Geant4-DNA results of our IRT implementation were compared to Geant4-DNA SBS public method by setting the parameters related to the time step (see Section 2): the current SBS-dynamic time step with varying time steps as shown in Table IV, the SBS-dynamic time step with the same minimum time step value of 0.1 ps. For these time step models, the same parameters for diffusion coefficients and reaction rate coefficients, as available in Geant4-DNA (see Table I), are used and the Brownian bridge technique is used.

The starting point that motivated the developments around the IRT method is the reduction of the otherwise long computation times, as it is known to be much faster than SBS methods. Our implementation of the IRT, coupled with the synchronous algorithm, has resulted in a gain in calculation time for 500 keV protons and one voxel of complex DNA geometry by a factor of 30 relative to *SBS-dynamic time step using Minimum Time Steps of 0.1 ps* and a factor of 15 relative to *SBS-dynamic time step using variable Minimum Time Steps*, the latter currently the Geant4-DNA default.
<figure><img src="image_9.png" /><figcaption>FIG. 3. The G-values as a function of time calculated by IRT (this work), Topas-IRT, Geant4-DNA-SBS-dynamic-time-step. Experimental data: OH*, ☐ Lav-erne, 2000,³⁶ ■ Jay-Gerin and Ferradini, 2000,³⁷ ● Jonah and Miller, 1977,³⁸ El Omar et al., 2011,³⁹ $e_{aq}^-$: Shiraishi et al., 1988,⁴⁰ ■ Sumiyoshi and Katayama, 1982,⁴¹ ● Hunt et al., 1973⁴² and Wolff et al., 1973,⁴³ ● Buxton, 1972,⁴⁴ Muroya et al., 2005,⁴⁵ $H_3O^+$: Pikaev et al., 1977,⁴⁶ ■ Cercek and Kongshaug, 1969,⁴⁷ Anderson et al., 1985,⁴⁸ ● Schmidt and Ander, 1969,⁴⁹ $H_2O_2$: × Laverne, 2000,³⁶ Elliot et al., 1993,⁵⁰ ● Appleby and Schwarz, 1969,⁵¹ $H_2$: Δ Draganic and Draganic 1975,⁵² ■ LaVerne and Pimblott (1991),⁵³ $H^*$: Δ Draganic and Draganic 1972.⁵⁴ [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

Figure 4 shows the damage yield produced per incident proton, at each value of the virtual time (not cumulated). For this calculations, $10^4$ initial protons of 500 keV were used traversing the voxel containing the chromatin fiber. Note that SBS-dynamic time step using Minimum Time Steps of 0.1 ps is selected as the reference data in this study. While the first maximum shows the direct damage and indirect damage at 2 ps (therefore just right after the physicochemical stage), the second one shows the “active” reaction time of indirect damage of the simulation. We observe a reasonable agreement between the three models from 1 ps to 2 ns and from 20 to

100 ns. The discrepancies raise up in the range 2–20 ns between SBS-dynamic time step using Minimum Time Steps of 0.1 ps and SBS-dynamic time step using variable Minimum Time Steps. The latter shows a maximum difference of about 15% at 6 ns. This can be explained by inaccuracies in the counting of damage due to a too long minimum time step (10 ps in the range 1–10 ns).

For the IRT method (both schemes), we observe a difference of about 35% at 6 ns in comparison with SBS-dynamic time step of 0.1 ps. The difference may be attributed to the higher reactivity of reactants and faster decrease than the SBS method for radical species as $e_{aq}^-$, $H^*$, $\cdot OH$ over the virtual time that we observe, as well, by their G-values in Fig. 3. As they participate in causing DNA damage (see Table III), their low concentration from 2 ns leads to the lower yield of DNA damage in comparison with SBS-dynamic time step of 0.1 ps.

In order to compare the methods, we chose the total number of damages as the key quantity. Figure 5 shows the cumulative damage expressed in terms of damage per primary particle as a function of virtual time obtained by integrating the data of Fig. 4. Consequently, there is a good agreement

<table><caption>TABLE IV. Variable minimum time steps proposed by Kreipl et al.²⁸</caption><thead><tr><td>Slice of simulated time (ps)</td><td>Time steps (ps)</td></tr></thead><tbody><tr><td>1–10</td><td>0.1</td></tr><tr><td>10–100</td><td>1</td></tr><tr><td>100–10<sup>3</sup></td><td>3</td></tr><tr><td>10<sup>3</sup>–10<sup>4</sup></td><td>10</td></tr><tr><td>Above 10<sup>4</sup></td><td>100</td></tr></tbody></table>
<figure><img src="image_10.png" /><figcaption>FIG. 4. DNA damage yields (damages/particle) including direct (plotted at 2 ps) and indirect damage including DNA backbone and bases (second part of the curve starting from 4 ps) as a function of virtual time. Results are shown for SBS-dynamic time step using 0.1 ps minimum time step (thin black line), SBS-dynamic time step using minimum time step as a function of virtual time from Table IV (dash-dot black line), IRT (dashed red line) from this work, and IRT using only totally diffusion-controlled reaction are shown for comparison. [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

<figure><img src="image_11.png" /><figcaption>FIG. 5. Cumulative damage yield as a function of virtual time corresponding to the results. [Color figure can be viewed at wileyonlinelibrary.com]</figcaption></figure>

of the cumulative damage between all methods for simulation times up to 3 ns. If we consider the IRT method we developed, the difference is of the order of 10% for higher simulation times due to the underestimation of reactions in the time

range 2–20 ns as explained above. This difference is less than the uncertainty generally obtained in experimental DNA damage assessment some minutes after irradiation (as experimental data on indirect effects some ns after the irradiation are not possible to perform³⁴). These results should be interpreted considering the fact that in concrete applications (e.g., cellular irradiations), the simulation time of the chemistry is an adjustable parameter set to reproduce scavenging effects that limit the diffusion of chemical species. For example, the simulation time of the chemical step with the SBS method of Geant4-DNA was limited to 1 ns for the calculation of damage to a chromatin fiber³³ or 2.5 ns for damage to a complete cell nucleus.²,³³ In the latter case, it was suggested that a time of 10 ns would be more appropriate.³⁵ Under these conditions, it must therefore be considered that the IRT method presented in this study gives results that are reasonably close to those of the method we have considered as a reference for Geant4-DNA calculations about 10% at 10 ns (SBS-dynamic time step using 0.1 ps minimum user time step).

## 4. CONCLUSIONS

The IRT-based time step method developed in this work was successfully implemented in the chemistry module of Geant4-DNA. This implementation allows calculation of the positions of the particles at each time step. It may be used in the simulations of complex distributions of molecules and radiolytic species that evolve with time as DNA is damaged. The development succeeded in our objective: keeping the efficiency in computation time of the IRT method while taking into account the constraints related to the calculation of DNA damage by having access to the position of reactants in the DNA geometrical frame. The obtained G-values as function of virtual time for *OH, $e_{aq}^-$, H₂, and H₂O₂ produced by 80 keV electron tracks are in generally good agreement with experimental data, considering the experimental uncertainty, and with results from IRT and SBS simulations, validating the developed method. By using a piece of heterochromatin fiber, consisting of 3640 nucleotide pairs, we have demonstrated the use of our IRT implementation in a complex DNA geometry. We observe that IRT shows a difference of about 35% decrease in DNA damage yields at 6 ns in comparison with SBS-dynamic time step (Minimum User Time Step of 0.1 ps). Despite this difference, the total number of damages calculated over a typical simulation time of the chemical part not exceeding 10 ns agrees within 10% with the SBS-dynamic time step method currently available in Geant4-DNA. This difference is acceptable, considering the uncertainty in the SBS-dynamic time step results and the large gain in computational time of a factor of 30. This gain is particularly valuable for radiobiology applications. In the near future, the potential of the simulation of the chemical step of Geant4-DNA will be extended not only by the addition of this IRT variant but also by the addition of the traditional IRT method made available by the Geant4-DNA collaboration. In this frame, two other IRT versions have been developed to cover several types of applications and will be released after their publication.
Investigations on the conditions of validity of each method and comparisons should be made available to allow Geant4-DNA users to choose, improve, or propose the method that best suits their needs. A comparison of simulations with experimental data on plasmid geometries could be a first step.

## ACKNOWLEDGEMENT

JRM and BF were supported by the National Institutes of Health/National Cancer Institute (NIH/NCI grant no. R01 CA187003: "TOPAS-nBio: A Monte Carlo tool for radiation biology research").

## REFERENCES

1. Mavragani IV, Nikitaki Z, Kalospyros SA, Georgakilas AG. Ionizing radiation and complex DNA damage: from prediction to detection challenges and biological significance. *Cancers*. 2019;11:1789.

2. Meylan S, Incerti S, Karamitros M, et al. Simulation of early DNA damage after the irradiation of a fibroblast cell nucleus using Geant4-DNA. *Sci Rep*. 2017;7:11923.

3. Lampe N, Karamitros M, Breton V, et al. Mechanistic DNA damage simulations in Geant4-DNA part 1: a parameter study in a simplified geometry. *Phys Med*. 2018;48:135-145.

4. Lampe N, Karamitros M, Breton V, et al. Mechanistic DNA damage simulations in Geant4-DNA part 2: electron and proton damage in a bacterial cell. *Phys Med*. 2018;48:146-155.

5. Sakata D, Lampe N, Karamitos M, et al. Evaluation of radiation DNA damage in a fractal cell nucleus model using Geant4-DNA. *Phys Med*. 2019;62:152-157.

6. Meylan S, Vimont U, Incerti S, Clairand I, Villagrasa C. Geant4-DNA simulations using complex DNA geometries generated by the DnaFabric tool. *Comput Phys Commun*. 2016;204:159-169.

7. Buxton GV, Greeanstock CL, Helman WP, Ross AB. Critical review of rate constants for reactions of hydrated electrons, hydrogen atoms and hydroxyl radicals (·OH/·O⁻) – in aqueous solution. *J Phys Chem Ref Data*. 1988;17:513-886.

8. Bluett VM, Green NJB. Competitive diffusion-influenced reaction of a reactive particle with two static sinks. *J Phys Chem A*. 2006;110:4738-4752.

9. Agmon N. Diffusion with back reaction. *J Chem Phys*. 1984;8:2811.

10. Clifford P, Green NJB, Oldfield MJ, Pilling MJ, Pimblot SM. Stochastic models of multi-species kinetics in radiation-induced spurs. *J Chem Soc Faraday Trans*. 1986;82:2673-2689.

11. Karamitros M, Luan S, Bernal MA, et al. Diffusion-controlled reactions modeling in Geant4-DNA. *J Comput Phys*. 2014;274:841-882.

12. Uehara S, Nikjoo H. Monte Carlo simulation of water radiolysis for low-energy charged particles. *J Radiat Res*. 2006;47:69-81.

13. Green NJB, Pilling MJ, Pimblot SM, Clifford P. Stochastic modeling of fast kinetics in a radiation track. *J Phys Chem*. 1990;94:251-258.

14. Frongillo Y, Goulet T, Fraser M-J, Cobut V, Patau JP, Jay-Gerin J-P. Monte Carlo simulation of fast electron and proton tracks in liquid water —II. *Nonhomogen Chem Radiat Phys Chem*. 1998;51:245-254.

15. Clifford P, Green NJB, Pilling MJ. Monte Carlo Simulation of diffusion and reaction in radiation-induced spurs. Comparison with analytic models. *J Phys Chem*. 1982;86:1322-1327.

16. Agostinelli S, Allison J, Amako K, et al. Geant4—a simulation toolkit. *Nucl Instrum Methods Phys Res A*. 2003;506:250-303.

17. Allison J, Amako K, Apostolakis J, et al. Geant4 developments and applications. *IEEE Trans Nucl Sci*. 2006;53:270.

18. Allison J, Amako K, Apostolakis J, et al. Recent developments in GEANT4. *Nucl Instrum Methods Phys Res A*. 2016;835:186-225.

19. Incerti S, Baldacchino G, Bernal M, et al. The Geant4-DNA project. *Int J Model Simul Sci Comput*. 2010;1:157-178.

20. Incerti S, Ivanchenko A, Karamitros M, et al. Comparison of Geant4 very low energy cross section models with experimental data in water. *Med Phys*. 2010;37:4692-4708.

21. Bernal MA, Bordage MC, Brown J, et al. Track structure modeling in liquid water: a review of the Geant4-DNA very low energy extension of the Geant4 Monte Carlo simulation toolkit. *Phys Med*. 2015;31:861-874.

22. Incerti S, Kyriakou I, Bernal MA, et al. Geant4-DNA example applications for track structure simulations in liquid water: a report from the Geant4-DNA Project. *Med Phys*. 2018;45:722-739.

23. Schuemann J, McNamara AL, Ramos-Méndez J, et al. TOPAS-nBio: an extension to the TOPAS simulation toolkit for cellular and sub-cellular radiobiology. *Radiat Res*. 2018;191:125-138.

24. Perl J, Shin J, Schümann J, Faddegon B, Paganetti H. TOPAS: an innovative proton Monte Carlo platform for research and clinical applications. *Med Phys*. 2012;39:6818-6837.

25. Ramos-Méndez J, Shin WG, Karamitros M, et al. Independent reaction times method in Geant4-DNA: implementation and performance. *Med Phys*. 2020;47:5919-5930.

26. Karamitros M, Brown J, Lampe N, et al. Implementing the Independent Reaction Time method in Geant4 for radiation chemistry simulations. *arXiv:2006.14225*.

27. Udovičić LJ, Mark F, Bothe E, Udovicić LJ. Yields of single-strand breaks in double-stranded calf thymus DNA irradiated in aqueous solution in the presence of oxygen and scavengers. *Radiat Res*. 2006;140:166.

28. Kreipl MS, Fridland W, Paretzke HG. Time- and space-resolved Monte Carlo study of water radiolysis for photon, electron and ion irradiation. *Radiat Environ Biophys*. 2009;48:11-20.

29. Goulet T, Jay-Gerin J-P. On the reactions of hydrated electrons with OH- and H3O+. Analysis of photoionization experiments. *J Chem Phys*. 1992;96:5076.

30. Plante I, Devroye L. Considerations for the independent reaction times and step-by-step methods for radiation chemistry simulations. *Radiat Phys Chem*. 2017;139:157-172.

31. Shin W-G, Ramos-Mendez J, Faddegon B, et al. Evaluation of the influence of physical and chemical parameters on water radiolysis simulations under MeV electron irradiation using Geant4-DNA. *J Appl Phys*. 2019;126:104301.

32. Delage E, Pham QT, Karamitros M, et al. PDB4DNA: implementation of DNA geometry from Protein Data Bank (PDB) description for Geant4-DNA Monte-Carlo simulations. *Comput Phys Commun*. 2015;192:282-288.

33. Tang N, Bueno M, Meylan S, et al. Influence of chromatin compaction on simulated early radiation-induced DNA damage using Geant4-DNA. *Med Phys*. 2019;46:1501-1511.

34. Kinner A, Wu W, Staudt C, Iliakis G. γ-H2AX in recognition and signaling of DNA double-strand breaks in the context of chromatin. *Nucleic Acids Res*. 2008;36:5678-5694.

35. Tang N, Bueno M, Meylan S, et al. Assessment of radio-induced damage in endothelial cells irradiated with 40 kVp, 220 kVp, and 4 MV x-rays by means of micro and nanodosimetric calculations. *Int J Mol S*. 2019;20:6204.

36. LaVerne JA. OH radicals and oxidizing products in the gamma radiolysis of water. *Radiat Res*. 2000;153:53-196.

37. Jay-Gerin J-P, Ferradini C. A new estimate of the radical yield at early times in the radiolysis of liquid water. *Chem Phys Lett*. 2000;317:388-391.

38. Jonah CD, Miller JR. Yield and Decay Of The Hydroxyl Radical From 200ps to 3ns. *J Phys Chem*. 1977;81:1974-1976.

39. Omar AK, El S, Uli J, et al. Time-dependent radiolytic yield of OH• radical studied by picosecond pulse radiolysis. *J Phys Chem A*. 2011;115:12212-12216.

40. Shiraishi H, Katsumura Y, Hiroishi D, Ishigure K, Washio M. Pulse-radiolysis study on the yield of hydrated electron at elevated temperatures. *J Phys Chem*. 1988;92:3011-3017.

41. Sumiyoshi T, Katayama M. The yield of hydrated electrons at 30 Picoseconds. *Chem Lett*. 1982;11:1887-1890.

<footnote>a) Author to whom correspondence should be addressed. Electronic mail: ngochoang.tran.vn@gmail.com.</footnote>
42. Hunt JW, Wolff RK, Bronskill MJ, Jonah CD, Hart EJ, Matheson MS. Radiolytic yields of hydrated electrons at 30 to 1000 picoseconds after energy absorption. *J Phys Chem.* 1973;77:425-426.

43. Wolff RK, Bronskill MJ, Aldrich JE, Hunt JW. Picosecond pulse radiolysis. IV. Yield of the solvated electron at 30 picoseconds. *J Phys Chem.* 1973;77:1350-1355.

44. Buxton GV. Nanosecond Pulse Radiolysis of Aqueous Solutions Containing Proton and Hydroxyl Radical Scavengers. Proceedings of the Royal Society of London A: Mathematical, Physical and Engineering Sciences; 1972: 328. https://doi.org/10.1098/rspa.1972.0065

45. Muroya Y, Lin M, Wu G, et al. A re-evaluation of the initial yield of the hydrated electron in the picosecond time range. *Radiat Phys Chem.* 2005;72:169-172.

46. Pikaev AK, Kabakchi SA, Zansokhova AA. Yields and reactions of hydrogen ions on radiolysis of water and aqueous solutions. *Faraday Discuss Chem Soc.* 1977;63:112-123.

47. Cercek B, Kongshaug M. Hydrogen ion yields in the radiolysis of neutral aqueous solutions. *J Phys Chem.* 1969;73:2056-2058.

48. Anderson RF, Vojnovic B, Michael BD. The radiation-chemical yields of H<sub>3</sub>O<sup>+</sup> and OH<sup>-</sup> as determined by nanosecond conductimetric measurements. *Radiat Phys Chem.* 1977;1985:301-303.

49. Schmidt KH, Ander SM. Formation and recombination of the hydrogen ion (H<sub>3</sub>O<sup>+</sup>) and hydroxide in irradiated water. *J Phys Chem.* 1969;73:2846-2852.

50. Elliot AJ, Chenier MP, Ouellette DC, et al. Temperature dependence of G values for H<sub>2</sub>O and D<sub>2</sub>O irradiated with low linear energy transfer radiation. *J Chem Soc Faraday Trans.* 1993;89:1193-1197.

51. Appleby A, Schwarz HA. Radical and molecular yields in water irradiated by gamma-rays and heavy ions. *J Phys Chem.* 1969;73:1937-1941.

52. Draganic ZD, Draganic IG. Formation of primary reducing yields (Geaq- and GH2) in the radiolysis of aqueous solutions of some positive ions. *Int J Radiat Phys Chem.* 1975;7:381-386.

53. LaVerne JA, Pimblott SM. Scavenger and time dependences of radicals and molecular products in the electron radiolysis of water: examination of experiments and models. *J Phys Chem.* 1991;95:3196-3206.

54. Draganic ZD, Draganic IG. Formation of primary hydrogen atom yield (GH) in the gamma radiolysis of water. *J Phys Chem.* 1972;76:2733-2737.