**1. Introduction**

Understanding biological responses to ionizing radiation is of crucial importance for cancer treatment using radiotherapy. Mechanistic Monte Carlo (MC) simulation of radiation's effect on DNA in a water medium is a promising tool for relevant studies after decades of development [1]. The central idea of such an approach is to obtain the initial DNA damage spectrum via mechanistic modeling of the radio-biological interactions at the atomic or molecular levels. This includes the development of track-structure codes [2–15] and the subsequent computation of DNA damages by incorporating DNA models [5,7,16–18]. The track-structure simulation can be divided into the simulations of the physical stage and the chemical stage. The physical-stage simulation deals with the ionization, excitation, and elastic scattering processes between the ionizing radiation particles and the water media and records the 3D coordinates of energy deposition events. The chemical-stage simulation computes how the chemical radicals, produced after the physical-stage simulation, diffuse and react mutually with the recording of the residual radicals' positions. The positions of these energy deposition events and radicals are then utilized to compute the initial DNA damage sites, followed by an analysis to characterize DNA strand break patterns.

**Citation:** Lai, Y.; Jia, X.; Chi, Y. Recent Developments on gMicroMC: Transport Simulations of Proton and Heavy Ions and Concurrent Transport of Radicals and DNA. *Int. J. Mol. Sci.* **2021**, *22*, 6615. https://doi.org/ 10.3390/ijms22126615

Academic Editor: Małgorzata Borówko

Received: 18 May 2021 Accepted: 16 June 2021 Published: 21 June 2021

**Publisher's Note:** MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Copyright:** © 2021 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/).

Although many developments have been performed to generate state-of-the-art mechanistic MC simulation tools, it is still necessary to further improve the simulation methods to accommodate different scenarios [8]. For instance, to make the code versatile for studies on the Oxygen Enhancement Ratio (OER) [19] and the Fenton reaction effect [20], it is desired to include more types of molecules other than free radicals generated by the initial radiation into the chemical-stage simulation. However, due to the computational complexity of the "many-body" problem and the long temporal duration of the chemical stage, a step-by-step simulation of these relevant processes on conventional CPU computational platforms can be extremely time consuming [21]. Under the constraint of computational resources, studies typically suffer from a restricted simulation region or a shortened temporal duration [22–24], limiting their broad applications.

To overcome these obstacles, Graphical Processing Unit (GPU)-based parallel computing can be a cost-effective option [25,26]. We developed an open-source, GPU-based microscopic MC simulation toolkit, gMicroMC [9], with the first version available on GitHub (https://github.com/utaresearch/gMicroMC (accessed on 17 June 2021)). We initially focused on boosting the chemical-stage simulation for radicals produced from water radiolysis, achieving a speedup of several hundred folds compared to CPU-based packages [27]. Later, we supported the physical track simulation for energetic electrons and implemented a DNA model of a lymphocyte cell nucleus at the base-pair resolution for the computation of electron-induced DNA damages [9]. Recently, we also included oxygen molecules in the chemical-stage simulation in a step-by-step manner, which enabled the study of the radiolytic depletion effect of dissolved oxygen molecules [28]. With these efforts, we were able to quantitatively study multiple critical problems that are computationally demanding. For example, we performed comprehensive simulations with gMicroMC to answer how uncertainties from the simulation parameters affect the accuracy of the final DNA damage computations [29]. We also studied the radiolytic depletion of oxygen under ultra-high dose rate radiation (FLASH) to investigate the fundamental mechanism behind FLASH radiotherapy with the developed oxygen module in gMicroMC [28].

In this work, we reported our recent progress on two new and important features that we recently introduced to gMicroMC, namely: (1) enabling the physical-stage simulations of protons and heavy ions; and (2) considering the presence of the DNA structure and its chemical reactions with radicals in the chemical-stage simulation. It was expected that the first feature would contribute greatly to the mechanistic study of particle irradiation, such as particle radiotherapy [30,31]. The presence of the DNA structure in the chemicalstage simulation will allow us to realistically describe the indirect DNA damage process. With the GPU acceleration, we were able to afford computationally challenging simulations that included detailed physics modeling and chemical reactions that spanned over a large temporal duration, enabling more realistic simulations of the relevant processes.

#### **2. Materials and Methods**

#### *2.1. Cross-Sections for the Transport Simulation of Protons and Heavy Ions*

When a proton or heavy ion moves through a medium, it interacts with the atomic electrons inside the medium [32]. Considering that there have been various models developed and implemented to describe this process, in this work, we focused on a novel implementation of existing models on GPU parallel computing platforms. Specifically, we only considered the interactions between particles and water molecules because this is representative of modeling the cell environment. We employed the Rudd model [33] to compute the ionization of a water molecule by a proton in the energy range from 10 eV to 1 TeV. We implemented the Plante model [34] and Dingfelder's model [35] to simulate the excitation of a water molecule for protons with an energy above and below 500 keV, respectively. We also applied Booth's empirical formula [36] to include the charge effect on the cross-section computation. Lastly, we used the charge scaling rule [37] to obtain the cross-sections for heavy ions based on those for a proton. To make the manuscript easy to follow, we briefly introduce these models in the following subsections.

#### 2.1.1. Ionization for Protons

An energetic proton could eject a secondary electron from different atomic subshells when it ionizes a water molecule. In the Rudd model [33,38], the partial Singly Differential Cross-Section (SDCS) can be described as:

$$\frac{d\sigma\_i^{ion}}{dw} = \frac{S\_i}{B\_i} \frac{\left(F\_1(\nu) + wF\_2(\nu)\right)}{\left(1 + w\right)^3} \frac{1}{1 + \exp\left[\frac{a\left(w - w\_i\right)}{\nu}\right]}\tag{1}$$

Here, *i* refers to the subshells of the water molecule, namely <sup>1</sup>*b*1, <sup>3</sup>*a*1, <sup>1</sup>*b*2, <sup>2</sup>*a*1, and <sup>1</sup>*a*1. *Bi* is the binding energy for electrons on shell *i*. *w* = *Ee*/*Bi*, and *Ee* is the energy of the secondary electron. *Si* = <sup>4</sup>*πa*20 ∗ *Ni* ∗ (*ER*/*Bi*)2, where *a*0 = 5.3 × 10−<sup>11</sup> *m* is the Bohr radius, *ER* = 13.6 eV is the Rydberg energy, and *Ni* is the number of electrons on shell *i*. *ν* = √*T*/*Bi* denotes the scaled velocity of the projectile, with *T* = *mM* ∗ *Ek*. m and M are the masses of the electron and proton, while *Ek* is the kinetic energy of the proton. *wi* = 4*ν*<sup>2</sup> − 2*ν* − *ER*/4*Bi* is the scaled cutoff energy, and *α* is a numerical parameter related to the relative size of the target molecule. The specific values for *Bi*, *Ni* and *α* are listed in Table 1. *<sup>F</sup>*1(*ν*) and *<sup>F</sup>*2(*ν*) are fitting functions, defined as:

$$F\_1(\nu) = L\_1 + H\_1 \tag{2}$$

$$F\_2(\nu) = \frac{L\_2 H\_2}{L\_2 + H\_2} \tag{3}$$

where:

$$L\_1 = \frac{C\_1 \nu^{D\_1}}{1 + E\_1 \nu^{D\_1 + 4}} \,\prime \tag{4}$$

$$H\_1 = \frac{A\_1 \ln\left(1 + \nu^2\right)}{\nu^2 + B\_1/\nu^2},\tag{5}$$

$$L\_2 = \mathbb{C}\_2 \nu^{D\_2},\tag{6}$$

$$H\_2 = \frac{A\_2}{\nu^2} + \frac{B\_2}{\nu^4}.\tag{7}$$

The values for the nine basic parameters *A*1, ... , *E*1 and *A*2, ... , *D*2 used in Equations (4)–(7) can be seen in Table 1. These values differed for inner shell orbitals and external orbitals, and an orbital was regarded as an inner one when its binding energy exceeded twice the binding energy of the least-tightly bound orbital [33].

**Table 1.** Parameters used in this work for Equations (1)–(7). Data were extracted from [35,38].


From Equation (1), we can calculate the total cross-section for subshell *i* as:

$$
\sigma\_i^{\rm ion} = \int\_0^{w\_m} \frac{d\sigma\_i^{\rm ion}}{dw} dw\_\prime \tag{8}
$$

where *wm* = *EmBi* and *Em* is the scaled maximum transferable energy from the proton to the ejected electron. The relativistic expression of *Em* was given by Plante et al. [34] as:

$$E\_m = \frac{2mc^2\left(\gamma^2 - 1\right)}{1 + 2\gamma\left(\frac{m}{M}\right) + \left(\frac{m}{M}\right)^2},\tag{9}$$

with:

$$\gamma = 1 + \frac{E\_k}{Mc^2} \tag{10}$$

and *c* is the speed of light. The relativistic format for the scaled velocity *ν* is then written as:

$$\nu^2 = \frac{mc^2}{2B\_i} [1 - \frac{1}{\gamma^2}].\tag{11}$$

With the ionization model and parameters determined under both the relativistic and nonrelativistic formalism, we could integrate Equation (8) numerically to obtain the ionization cross-section table for different subshells of a water molecule in a broad proton energy range. In our implementation, we computed the table for proton energies ranging from 10 eV to 1 TeV with a 0.01 increment along the logarithmic scale. It is worth mentioning that the computation of the cross-section table only needed to be computed once in an offline manner and was stored in a data file that could be loaded to GPU memory for the query of the cross-sections of any incident energy. More details of this usage are given in Section 2.3.1.

#### 2.1.2. Excitation for Proton

Due to the lack of experimental data, different models could have differential configurations of the excitation pathways. In our implementation, we adopted the threepathway model [39,40] containing *<sup>A</sup>*<sup>1</sup>*B*1 and *<sup>B</sup>*<sup>1</sup>*A*1 and plasma mode for protons with energy >500 keV and the model with the excitation channels of *A* <sup>1</sup>*B*1 and *<sup>B</sup>*<sup>1</sup>*A*1, Ryd A + B and Ryd C + D, and the diffusion band [35] for a proton energy < 500 keV. Specifically, in the three-pathway model, the differential cross-section for the excitation channel *j* is expressed as:

$$\frac{d\sigma\_{\hat{j}}^{\text{ext}}}{d\mathcal{W}} = \rho(\mathcal{W})\mathcal{W}f\_{\hat{j}}(\mathcal{W})\ln\left(\frac{\mathcal{W}}{Q\_{\text{min}}}\right),\tag{12}$$

where:

$$
\rho(W) = \frac{8\pi Z^2 a\_0^2}{m u^2} \frac{R y^2}{W^2} \, , \tag{13}
$$

$$Q\_{\rm min} = 2T \left(\frac{M}{m}\right)^2 \left(1 - \frac{1}{2}\frac{m}{M}\frac{W}{T} - \sqrt{1 - \frac{m}{M}\frac{W}{T}}\right). \tag{14}$$

Here, *j* denotes different excitation channels, namely *A* <sup>1</sup>*B*1 and *<sup>B</sup>*<sup>1</sup>*A*1, and plasma mode. *u*, *Z*, and *W* are the velocity, charge, and energy loss of the incident proton. Other parameters were the same as those used in Section 2.1.1. When *mM WT* = *WEk* 1, Equation (12) can be simplified as:

$$\frac{d\sigma\_{\bar{\jmath}}^{\rm ext}}{d\mathcal{W}} = \rho(\mathcal{W})\mathcal{W}f\_{\bar{\jmath}}(\mathcal{W})\ln\left(\frac{4T}{\mathcal{W}}\right). \tag{15}$$

In the relativistic situation, *mu*<sup>2</sup> in Equation (13) can be expressed as:

$$
tau^2 = mc^2 \left[1 - \gamma^{-2}\right].\tag{16}$$

*fj*(*W*), as a function of the excitation pathway *j*, has the form of:

$$f\_{\hat{f}}(\mathcal{W}) = \begin{cases} f\_{\hat{\boldsymbol{\beta}}}^0 \sqrt{a\_{\hat{\boldsymbol{\beta}}}/\pi} e^{\left[-a\_{\hat{\boldsymbol{\beta}}} \left(\mathcal{W} - w\_{\hat{\boldsymbol{\beta}}}\right)^2\right]}, & \text{if } \boldsymbol{j} = \hat{A}^1 B\_1, \ \hat{B}^1 A\_1 \\ f\_{\hat{\boldsymbol{\beta}}}^0 a\_{\hat{\boldsymbol{\beta}}} e^{\mathbf{x}} / \left(1 + e^{\mathbf{x}}\right)^2, & \text{otherwise} \end{cases} \tag{17}$$

where *x* = *<sup>α</sup>j*(*<sup>W</sup>* − *wj*) and *f* 0*j* , *<sup>α</sup>j*, and *wj* are parameters with their values summarized in Table 2. Under the assumption that the proton only loses a small portion of its kinetic energy to excite a water molecule, i.e., *WEk* 1, Equation (15) can be used to calculate the total cross-section for excitation channel *j* as:

$$
\sigma\_j^{\text{exc}} = \int\_{\mathcal{W}\_{\text{min}}}^{\mathcal{W}\_{\text{max}}} \frac{d\sigma\_j^{\text{exc}}}{d\mathcal{W}} d\mathcal{W}. \tag{18}
$$

In principle, the upper and lower boundaries of the integration can be *Ek* and zero. However, in practical usage, it is common to set *Wmax* =50 eV and *Wmin* = 2 eV. The reason is that *dσexc j dW* (*W*) drops to a negligible value when *W* ∈/ [*Wmin*, *Wmax*], and the boundary cutoffs also ensure a positive and convergen<sup>t</sup> integrated total cross-section.

**Table 2.** Parameters used in Equations (17).


When a proton's energy is smaller than 500 keV, the Born approximation is no longer a good approximation [35], and Equation (15) may have problem in evaluating the crosssections. We then applied the semi-empirical model [35] to obtain the excitation crosssection for a low-energy proton as:

$$
\sigma\_j^{\text{exc}}(E\_k) = \frac{\sigma\_0 (Za)^{\Omega} \left(E\_k - E\_j\right)^v}{J^{\Omega + v} + E\_k^{\Omega + v}}.\tag{19}
$$

Here, *j* denotes the excitation channels *A* <sup>1</sup>*B*1 and *<sup>B</sup>*<sup>1</sup>*A*1, Ryd A+B and Ryd C+D, and the diffusion band. The corresponding energy loss *Ej* of the proton is discrete instead of continuous. Further details of the model can be found in [35].

With the excitation cross-section given in Equation (18) and the relevant parameters determined, we could integrate it numerically to obtain the excitation cross-section table for different subshells of a water molecule at proton energies above 500 keV. Meanwhile, we could rely on Equation (19) to handle protons with energies below 500 keV. To make the cross-section data computed from the two models smoothly connected at the proton energy of 500 keV, we adjusted the obtained cross-section data as follows. We applied coefficients of 1.23 and 3.5 to the cross-section data for the *A* <sup>1</sup>*B*1 and *<sup>B</sup>*<sup>1</sup>*A*1 channels obtained from Equation (18) to make them smoothly connected to that obtained from Equation (19) at 500 keV for these two modes. We then multiplied 0.339 in plasma mode to make the total cross-section also smoothly connected. Similar to the strategy applied to obtain the ionization cross-section table, we also only needed to compute the excitation cross-section table once and stored it in a data file. Its usage on a GPU is also given in Section 2.3.1.

## 2.1.3. Charge Effect

When charged particles travel through a water medium, except for ionizing or exciting the water molecule, they could also drag electrons from the medium to move with them, forming a reduced effective charge *Z*eff < *Z*. This effect is found reversely proportional to the kinetic energy of the incident particle. In our simulation, we adopted the empirical Booth model [36] to obtain the effective charge *Z*eff as:

$$Z\_{\rm eff} = Z \left( 1 - \exp\left[ -1.316y + 0.112y^2 - 0.065y^3 \right] \right),\tag{20}$$

where *y* = 100 *Z*−2/3 !1 − (1 + *Ek*/(*AMc*<sup>2</sup>)) −2 and *A* is the mass number of the particle. The correction is larger than 5% ( *Z*eff < 0.95 · *Z*) when *y* < 2.172, which gives *Ek* ∼ 18 MeV per nucleon for the Fe ion and 0.22 MeV for the proton.

#### 2.1.4. Cross-Section for Heavy Ions

Within the first-order plane-wave Born approximation, we could correlate the ionization and excitation cross-section for bare and sufficiently fast heavy ions to that of the proton by the scaling law. Specifically, for a heavy ion with velocity *u* and charge number *Z*, the doubly differential cross-section can be scaled from that of the proton with the same velocity *u* by a factor of *Z*<sup>2</sup> [37]:

$$\frac{d^2\sigma\_{ion}}{d\mathcal{W}d\underline{Q}}(\mathfrak{u}) = Z^2 \times \frac{d^2\sigma\_{proton}}{d\mathcal{W}d\underline{Q}}(\mathfrak{u}),\tag{21}$$

where *W* and *Q* refer to the energy transfer and the recoil energy, respectively. After integrating over *Q*, we could obtain the SDCS as a function of *W*. Considering that an ion of mass number *A* and kinetic energy *Ek* has the same velocity with a proton of kinetic energy *Ek*,*<sup>p</sup>* = *Ek M Mion*≈ *Ek A* , we could rewrite the scaling formula as a function of *Ek* as:

$$\frac{d\sigma\_{\rm ion}}{dW}(E\_k) = Z^2 \times \frac{d\sigma\_{\rm proton}}{d\mathcal{W}}(E\_k/A),\tag{22}$$

It holds for ions for both the nonrelativistic and relativistic formats. The electron attachment effect can be more significant for a heavy ion than for a proton of the same velocity since a heavy ion typically carries more charge than a proton. With the electron attachment effect considered, we replaced *Z* with *Z*eff when scaling the cross-section from a proton to a heavy ion using Equation (22), with *Z*eff calculated from Equation (20).

#### *2.2. Concurrent Transport Method*

Due to the computational challenge, existing MC tools compute the DNA damage formed by radical attachment typically via two successive steps. First, the radicals are diffused and mutually reacted in the chemical stage without DNA. Second, the coordinates of *OH*· radicals obtained at the end of the chemical stage are overlapped with the DNA geometry such that DNA damages caused by radicals can be computed [9,29]. We refer to this approach as the "overlay method". This approach is effective for simple applications. However, it can be problematic for those scenarios sensitive to radical evolution. In reality, DNA could be present and react with radicals during the radical diffusion. This could affect the radical yields and the damage site distribution on the DNA chain, consequently impacting the final characterization of the DNA strand breaks. To model this effect, in this work, we included the simulation of the reactions between radicals and DNA at the time of transporting the radicals in the chemical stage and refer to this approach as the "concurrent transport method".

In our previous development of gMicroMC without considering DNA in the chemical stage, we simulated the diffusion and reactions among radicals in a step-by-step fashion. The relevant parameters were the diffusion coefficient for each radical species and the

reaction rates for possible radical–radical reaction types. To include DNA in this transport frame, we need to know the diffusion coefficient of DNA and the reaction rates between radicals and DNA. Considering the relatively large mass of DNA, we assumed that the whole DNA molecule was static during the chemical transport and took its diffusion coefficient as zero to simplify the simulation. As for the reaction rates between radicals and DNA, we considered two types of reactions based on the DNA model for a whole lymphocyte cell nucleus [9]. The DNA was described in a voxel-based format with each voxel of side length 55 nm. The voxel was either empty or filled with a DNA chain that connected two faces of the voxel. The DNA chain consisted of a group of spheres representing the basic structures of the DNA: the base pair, the sugar-phosphate group, and the histone protein. With it, we considered the first reaction type as the damaging effect of *OH*· and *eh* radicals on the DNA bases and sugar-phosphate groups. The reaction rates are listed in Table 3. Here, although there were four types of DNA bases, that is Adenine (A), Guanine (G), Cytosine (C), and Thymine (T), associated with four different reaction rates with the *OH*· or *eh* radical, we used the average reaction rate in our simulation since our DNA model had no differentiation among the base types. The second type was the scavenging reaction for all radical species by the histone protein. In this reaction, the radical was assumed to be fully absorbed once it was within the histone protein volume.

**Table 3.** Reaction rates (×10<sup>9</sup> *L* · *mol*−<sup>1</sup> · *s*<sup>−</sup>1) used in gMicroMC for concurrent DNA transport [41].


\* A negative value means no reaction between the radical and the DNA substructure.

After introducing DNA into the chemical-stage simulation, two consequences required special attention. First, radicals were not supposed to be produced inside the DNA region; hence, at the beginning of the chemical stage, we needed to exclude those radicals produced inside the chromatin zone from the subsequent diffusion [42] without recording any damages on DNA. Second, as there were a huge number of DNA basic structures in our DNA model, for instance 6.2 × 10<sup>9</sup> base pairs, checking for reactions between radicals and DNA after each diffusion step would generate numerous computations. To circumvent the problem, we defined a time interval *ti* to control the frequency of checking for reactions between radicals and DNA. During the simulation of the chemical stage, as the time evolved, we only checked for radical–DNA reactions every *ti*. In the limit of a small *ti*, the frequent inspection for reactions ensured simulation accuracy. In the other extreme of a large *ti*, the DNA-related reactions would be less frequently checked, which eventually converged to the overlay picture. We study the impact of *ti* in later sections.

#### *2.3. GPU Implementation*

#### 2.3.1. Physical Transport for Protons and Heavy Ions

Before transporting protons and heavy ions on the GPU, we prepared lookup tables on the CPU host to store the tabulated cross-sections for a proton, as stated in Section 2.1. The lookup tables were then transferred to the GPU texture memory such that we could employ the GPU built-in fast interpolation technique to obtain cross-section data for particle transport. We supported two types of source particle generation: reading from a Phase Space File (PSF) or sampling from a given distribution. We sorted the source particles in descending order based on their charge number *Z*. If the particles were protons or heavy ions, we transported the sorted particles into groups using the GPU kernel we developed in this work dedicated to proton and heavy ion transport. If they were electrons, we transported them with our previously developed kernel for electron transport [9]. The purpose of particle sorting and grouping was to minimize the thread divergence on the GPU, and hence to improve simulation efficiency [25,26].

For the GPU kernel in charge of the transport of protons and heavy ions, each thread took care of one primary particle. For a particle with charge number *Z*, atomic mass number *A*, and incident energy *Ek*, the thread sampled its free travel distance *s* in water and its interaction with the water molecules in the iteration. Specifically, we first calculated the effective charge number *Z*eff according to Equation (20) and the kinetic energy *Ep* = *Ek A* for a proton with the same velocity as that of the primary particle. Based on the logarithmic value of *Ep*, we interpolated the lookup tables to obtain the cross-section *<sup>σ</sup>i*(*Ep*) for the proton. Here, *i* represents all ionization and excitation channels listed in the tables. We then scaled and summed *σi* to obtain the total cross-section for the primary particle as *σt* = *Z*2eff ∑ *σi* based on Equation (22). With *σt*, we sampled the free travel distance *s* in water as *s* = − *Mw ρ*·*σt*·*NA* ln *ζ*, where *ρ* and *Mw* are the density and atomic mass of water. *NA* is the Avogadro constant, and *ζ* is a random number uniformly distributed between zero and one. We forwarded the particle position by *s* along the momentum direction followed by a sampling of the interaction type based on the relative cross-section distribution *σi* ∑ *σj*.

 If the sampled interaction *i*0 was an ionization event, we then sampled the energy *Ee* and the emission angle of the ejected secondary electron, along with updating the kinetic energy of the primary particle. Noticing that the partial SDCS in Equation (1) had the form of *dσion i dw* ∝ *f*(*w*)*φ*(*w*), with *f*(*w*) = (*<sup>F</sup>*1(*ν*) + *wF*2(*ν*))/(<sup>1</sup> + *w*) 3, which was relatively easy to integrate, and *φ*(*w*) = 1/(1 + exp [*α*(*w* − *wi*)/*ν*]), we took *f*(*w*) as a sampling function and *φ*(*w*) as the rejection function to effectively sample *Ee*. Specifically, for a given proton energy *Ep*, *ν* can be solely determined based on Equation (11), and hence, *<sup>F</sup>*1(*ν*) and *<sup>F</sup>*2(*ν*) are just numbers. We wrote them as *F*1 and *F*2 for simplicity in the following equations. Applying the direct inversion method, we could sample *ws* from *f*(*w*) as:

$$w\_s = \left(-F\_1 + 2N\_c\zeta + \sqrt{F\_1^2 + 2F\_2N\_c\zeta - 2F\_1N\_c\zeta}\right)/(F\_1 + F\_2 - 2N\_c\zeta). \tag{23}$$

Here, *ζ* ∈ [0, 1) is a randomly sampled number and *Nc* has the form of:

$$N\_c = \int\_0^{w\_m} f(w) dw = \left(w\_m(F\_2 w\_m + 2F\_1 + F\_1 w\_m)\right) / \left(2(1 + w\_m)^2\right). \tag{24}$$

We repeatedly sampled *ws* with Equation (23) and updated *φ*(*ws*) until we obtained *φ*(*ws*) > *ζ* with *ζ* a random number ∈ [0, 1/(1 + *<sup>e</sup>*<sup>−</sup>*<sup>α</sup>wi*/*v*)]. Once reaching this stopping criterion, we accepted *ws* and computed *Ee* as:

$$E\_{\mathfrak{c}} = w\_{\mathfrak{s}} \* B\_{i\_{\mathfrak{l}}}.\tag{25}$$

The polar scattering angle *θe* of the electron relative to the moving direction of the primary particle satisfied cos *θe* = ! *Ee* 4∗*T* for *Ee* > *Bi*0 and was uniformly distributed between zero and *π* otherwise [40]. The azimuth scattering angle was uniformly sampled between zero and 2 *π*. The residual energy of the primary particle after ionization was *E k* = *Ek* − *Ee* − *Bi*0 , and its polar scattering angle was zero.

If the sampled interaction *i*0 belonged to the excitation category, there was no secondary electron emission, and we only needed to sample the energy loss *W* of the primary particle. In this case, the polar scattering angle for the primary particle was zero as well. When *Ep* > 500 keV, we sampled *W* based on Equation (15). Noticing that *dσexc i* 0 *dW* ∝ *fi*0 (*W*)*g*(*W*), where *g*(*W*) = 1 *W ln*( 4*T W* ), we then used *fi*0 (*W*) for the sampling of *W* and *g*(*W*) for rejection. For the *<sup>A</sup>*<sup>1</sup>*B*1 and *<sup>B</sup>*<sup>1</sup>*A*1 channels, *Ws* can be directly sampled as:

$$\mathcal{W}\_{\mathfrak{s}} = w\_{\mathfrak{i}\_0} + \frac{1}{\sqrt{2\alpha\_{\mathfrak{i}\_0}}} \mathcal{J}\_{\mathfrak{n}\_{\mathfrak{i}}} \tag{26}$$

where *ζn* is a random number following the standard normal distribution. As for plasma mode, we applied the direct inversion method to obtain *Ws* as:

$$\mathcal{W}\_{\mathfrak{s}} = w\_{\dot{i}\_0} + \frac{1}{\varkappa\_{\dot{i}\_0}} \ln \left( \frac{u\_1}{1 - \varkappa\_{\dot{i}\_0} u\_1 \mathcal{N}\_{\mathfrak{c}} \zeta} - 1 \right). \tag{27}$$

Here, *ζ* is a random number ∈ [0, 1) and *Nc* has the form of:

$$N\_{\mathbb{C}} = \int\_{W\_{\text{min}}}^{W\_{\text{max}}} f\_{\mathbb{i}\_{\mathbb{0}}}(\mathcal{W}) d\mathcal{W} = \frac{1}{\mathfrak{a}\_{\mathbb{i}\_0}} \left( \frac{1}{\mathfrak{u}\_1} - \frac{1}{\mathfrak{u}\_2} \right), \tag{28}$$

with *u*1 = 1 + *eαi*0 (*Wmin*−*wi*0 ) and *u*2 = 1 + *eαi*0 (*Wmax*−*wi*0 ). We repeated the sampling of *Ws* and updating *g*(*Ws*) until we obtained *g*(*Ws*) > *ζ* with *ζ* a random number ∈ [0, 1 *Wmin* ln 4*T Wmin* ]. The residual energy of the primary particle after excitation was then *Ek* = *Ek* − *Ws*. When *Ep* ≤ 500 keV, the energy loss *Ei*0 was directly obtainable from the discrete energy spectrum [37]. The residual energy of the primary particle was then *Ek* = *Ek*− *Ei*0 .

After transporting the primary particle with one step and simulating its interaction with one water molecule, we updated *Ek* with *Ek* and started the next round of transport sampling until the kinetic energy of the primary particle reached the cutoff energy or ran out in the simulation region. During the process, all secondary electrons were stored in a global stack to be further simulated using our previously developed kernel in charge of electron transport, the physics models that covered the electron spectrum as low as a few eV [9]. The ionized and excited water molecules were also tagged for further analysis.

#### 2.3.2. Concurrent Transport

In the concurrent transport picture, we simulated the reactions among radicals and DNA and the diffusion of the radicals in a step-by-step manner. Considering the complex structure of DNA and the possibly different checking frequencies for radical–DNA interactions and radical–radical reactions depending on the value of *ti*, we utilized two GPU kernels for the chemical stage transport in the presence of DNA. One GPU kernel was responsible for the interactions between radicals and DNA, and the other kernel was in charge of the radical–radical reactions and the diffusion of the radicals.

For the GPU kernel managing the reactions between radicals and DNA, each GPU thread was responsible for one radical. To obtain the possible reaction and reaction type between the radical and DNA, we needed to search the DNA geometry and compute the distances from the radical to the centers of the DNA basic structures (DNA base, sugarphosphate group, and histone protein). The smallest distance *dmin* was then compared to the reaction range of *R* + *Rc* with *R* the radius of that DNA structure. *Rc* = *k*/4*πNAD* for reactions between the radical and DNA base or sugar-phosphate group, where *k* is the reaction rate, *NA* is the Avogadro constant, and *D* is the diffusion rate for the radical. For all considered radical species, *Rc* was < 1 nm. Due to the lack of experimental data for the reaction between radicals and the histone protein, we assumed that the radical was only absorbed when it hit the histone, and hence, we set *Rc* = 0 for this case. If *dmin* < *R* + *Rc*, a reaction was recorded. Otherwise, the Brownian bridge method [9] was applied to compensate for possible reactions between the radical and DNA during the diffusion. As our DNA model was constructed with a huge amount of basic structures, it would be too time consuming to search the entire space to obtain the smallest distance from the radical to the DNA. To reduce the searching burden, we relied on the voxelized geometry of the DNA model and only performed the search at most on two voxels. Specifically, noticing that the outer boundary of the DNA chain was > 2 nm away from all edges of the voxel it occupied [9] and *Rc* < 1 nm for all reactions between the radical and DNA, this indicated that a radical could only react with those DNA structures in two special voxels: the current voxel in which it was located and the adjacent voxel having the surface closest to it. The latter voxel was not considered unless the radical was < 2 nm from its closest surface. In this way, we reduced the searching burden significantly. Once a reaction was recorded, the radical was removed from the reactant list. If the reaction was with the

DNA base or sugar-phosphate group, the reaction site was stored in a global stack for further analysis.

For the radical–radical reactions and radical diffusion, we continued to employ the GPU kernel developed in our previous work [9]. Each thread was in charge of one radical. To reduce the searching burden for mutual reactions, the entire space was divided into small grids with the grid size twice the largest reaction radius. This ensured that each radical only reacted with other radicals in the same or adjacent grids. The distances from the radical to other radicals were then computed and compared to the reaction radii to obtain whether a reaction would occur. If a reaction happened, the new products were placed, and radical–radical reactions were checked again. Otherwise, the radical was diffused by one step followed by the check for radical–radical reactions based on the Brownian bridge method.

At the beginning of the chemical stage, the GPU kernel for the DNA–radical reactions was executed to remove the radicals within the chromatin region from the subsequent chemical-stage simulation. This was followed by the launch of the GPU kernel in charge of the radical–radical reaction and radical diffusion. After that, we compared *<sup>t</sup>*elap, the time elapsed from the last execution of the former GPU kernel, to *ti*. If *<sup>t</sup>*elap ≥ *ti*, the former and latter kernels were called in sequence. Otherwise, only the latter kernel was executed. The process was repeated until reaching the end of the chemical stage.

## *2.4. Simulation Setup*

2.4.1. Simulation Setup for the Transport of Protons and Heavy Ions

We performed a series of simulations to validate the physical-stage transport for protons and heavy ions. These included: (1) the computations of the cross-section, linear energy transfer (LET), and traveling range; (2) the validation of the energy spectrum of secondary electrons, the radial dose distribution, and the track structure; and (3) the evaluation of the DNA damage spectrum.

We first calculated the total cross-section for both the ionization and excitation channels according to Equations (8) and (18), under the relativistic formats. The results were compared to Plante et al.'s [34] and Dingfelder et al.'s [35] works, as shown in Figure 1. Based on this, we calculated the track-length-averaged unrestricted LET for different ion species with their energy ranging from 0.1 ∼ 10<sup>4</sup> MeV amu<sup>−</sup>1. For an ion with energy *Ek*, we sampled the energy loss of primary particles *εi* and the free-fly distance *si*. We then repeated the simulation *N* = 10<sup>5</sup> times and computed the length-averaged unrestricted LET as:

$$LET = \sum\_{i=1}^{N} \frac{\varepsilon\_i}{s\_i} \cdot \frac{s\_i}{\sum\_{j=1}^{N} s\_j} = \frac{\sum\_{i=1}^{N} \varepsilon\_i}{\sum\_{j=1}^{N} s\_j}. \tag{29}$$

We compared the LETs to those reported by Plante et al. [34]. After that, we simulated the proton range by tracking its starting and ending positions for a proton energy of 0.1 ∼ 100 MeV and compared this to the NIST data. We show both results in Figure 2.

As for the validation of the energy spectrum of secondary electrons, we simulated the interactions of a 5 MeV proton and 4 MeV alpha particles with a liquid water target, recorded the energy of the secondary electrons, and compared it to those obtained with the GEANT4-DNA simulation [43] (GEANT4 Version 10.5.1). The result is plotted in Figure 3. As for the radial dose distribution, we transported 5 and 10 MeV protons within a water slab of 10 μm in thickness and infinitely long at the other two dimensions and analyzed the dose distribution within a thin slice 4 to 6 μm away from the proton starting point along the thickness direction. We accounted for the dose distributed in an annular ring with inner and outer radii of *r* and *r* + Δ*r* as the dose at radius *r*. We set Δ*r* = 1 nm, the same as that used in Wang et al.'s work [44]. We repeated the simulation 10<sup>5</sup> times, averaged the obtained radial dose, and compared it with that reported in Wang et al.'s work [44]. Finally, we show a representative physical track structure for a 5 MeV proton in Figure 5, including the track for both the primary proton and secondary electrons.

We used the lymphocyte nucleus model developed in our previous work for this evaluation study [9]. We initiated a proton emission plane of 11 × 11 μ m<sup>2</sup> and 5.5 μm away from the center of the cell nucleus for two proton energies, 0.5 and 0.9 MeV. For each energy, we randomly sampled the proton position at the emission plane and its momentum towards the positive z direction, transported the proton until reaching a cutoff energy of 1 keV, and recorded the dose inside the cell nucleus. We repeated the simulation until reaching the accumulated dose of 2 Gy. After that, we simulated the physio-chemical and chemical stage with a chemical stage duration of *tc* = 1 ns. We then applied the overlay method to obtain the DNA damage sites and grouped them into DNA Single-Strand Breaks (SSBs) and Double-Strand Breaks (DSBs) [29]. The result was compared to Nikjoo et al.'s work with the KURBUC model [45] and shown in Table 4.

#### 2.4.2. Simulation Setup for Concurrent Transport

We studied the impact of *ti* on the radical yields. We simulated the cases with a chemical stage duration of *tc* = 1, 10 ns, and 1 μs and *ti* from 1 ps to 1 μs with an increment of one at the logarithmic ten scale. Again, the lymphocyte cell nucleus with a radius of 5.5 μm was used as the Region Of Interest (ROI). As for the radical yield, we transported a 4.5 keV electron with its initial position randomly sampled inside the ROI and its direction towards the ROI center. We then took the generated radicals as inputs for the chemicalstage simulation. The final G values for the *eh*, *OH*·, *H*·, and *H*2*O*2 radicals were recorded. We repeated the simulation 100 times to reduce the statistical uncertainties and reported the averaged G values over all the simulations. The results are shown in Figure 6.

We also computed the DNA damage as a function of the incident proton energy and the chemical stage duration under the concurrent DNA transport frame. A proton energy *Ek* of 0.5, 0.6, 0.8, 1.0, 1.5, 2, 5 10, 20, and 50 MeV and a chemical stage duration *tc* of 1, 2.5, and 10 ns were considered, following the parameters used in Zhu et al.'s work [42]. We initiated the proton on a spherical shell with a radius of 5.5 μm and shot it randomly towards the inner space of the sphere. We repeated the simulation until having the accumulated dose in the sphere of 1 Gy. We then simulated the chemical stage with DNA concurrent transport (*ti* = 1 ps) and computed the total DSB yield. For each proton energy, using the DSB yield at *tc* = 1 *ns* as a reference, we defined *R*(*t*) = *NDSB*(*tc* = *t*)/*NDSB*(*tc* = 1 *ns*) to represent the relative DSB yields at *tc* = *t*. For each pair of *Ek* and *tc*, we ran the simulation 20 times and computed the mean and standard deviation for the relative DSB yield. We then compared the data with *tc* = 2.5 ns (R(2.5)) and 10 ns (R(10)) to Zhu et al.'s work [42] and show the results in Figure 7.
