**3. Results**

#### *3.1. Validation of Development for Protons and Heavy Ions*

Figure 1 presents the total and partial cross-sections for ionization and excitation as a function of incident proton energy. From Figure 1a, the total cross-section for ionization from our simulation agreed well with that from Plante et al.'s work [34]. From Figure 1b, for a proton energy > 500 keV, our simulated total cross-section for excitation matched that from Plante et al.'s work [34]. As for the slow proton, it followed that from Dingfelder et al.'s work [35]. The results revealed that the ionization model and the two-stage excitation model were successfully implemented as expected.

In Figure 1b, we noticed a dramatic drop-off of the total excitation cross-section at around 10 keV for the Plante model. This is due to the cross-section formula shown in Equation (15) depending on the scaled energy *T* = *m M Ek*. When *Ek* drops below 10 keV, *T* is too small to excite even the lowest excitation channel (*j* = *A*<sup>1</sup>*B*1). After replacing it with Dingfelder's model (Equation (19)) at the low-energy region, the excitation cross-section drops much more slowly. Considering that the low-energy proton largely appears after the Bragg peak, a proper excitation model could be important for the distal dose computation in proton therapy.

**Figure 1.** Total and partial cross-sections of (**a**) ionization and (**b**) excitation channels for protons with different energies.

The calculated unrestricted LETs for different ions are plotted in Figure 2a. They agreed well with Plante et al.'s work for ions with an energy greater than 0.5 MeV per nucleon. At the low energy range, LETs from our simulation were lower than those from Plante et al.'s work, which can be explained by the different excitation models used in the two simulations. As shown in Figure 1b, the excitation cross-section from our work was higher than that from Plante et al.'s work at the low energy range, resulting in a higher sampling rate of excitation interactions in our simulation. Considering that the energy loss from an excitation event was typically smaller than that from an ionization event (Table 1), a higher sampling of excitation could result in a lower LET. In Figure 2b, we show the proton range from our simulation and its comparison with the NIST data. As is shown, our simulation result matched well with the NIST PSTAR data (https://physics.nist.gov/PhysRefData/Star/Text/PSTAR.html (accessed on 17 June 2021)), with the relative difference smaller than 1%.

**Figure 2.** (**a**) The unrestricted LETs for different ions with different energies. The unit amu<sup>−</sup><sup>1</sup> means per nucleon. Solid lines represent data extracted from Plante et al.'s work, while data with diamond symbols are from our simulation with gMicroMC. (**b**) The simulated proton range for different energies.

Figure 3 shows the energy spectrum of secondary electrons generated from a 5 MeV proton and a 4 MeV alpha particle. From the figure, the yielding rates of secondary electrons dropped quickly along with the increase of the electron energy. For the entire plot, our simulated results agreed well with that from GEANT4-DNA. We did not compare the spectrum for electron energy greater than 200 eV due to a too low yielding rate and the consequent large uncertainty.

**Figure 3.** Secondary electron spectrum for (**a**) a 5 MeV proton and (**b**) a 4 MeV alpha particle.

In Figure 4, we see the radial dose distributions for 10 and 50 MeV protons from ours and Wang et al.'s work (Equations (1)–(7) in [44]) under the same setup. As is shown, the two curves matched quite well, although the curve from our simulation suffered a relatively large statistical fluctuation for the regions > 1000 nm from the primary track axis. Figures 3 and 4 together validated the energy spectrum and angular distribution of secondary electrons from our simulation, furthering proving the successful implementation of the transport models for protons and heavy ions.

**Figure 4.** Radial dose distributions for (**a**) 10 MeV and (**b**) 50 MeV protons.

We then present a track structure for a 5 MeV primary proton and its produced secondary electrons in liquid water in Figure 5. For simplicity, we only present the entrance (Figure 5a) and Bragg peak (Figure 5b) regions. At the entrance region, the secondary electron tracks were quite sparse. In contrast, they were much denser in the Bragg peak region. In addition, the electron track lengths were shorter in the Bragg peak region. This was mainly due to the kinetic energy of the proton being much higher at the entrance than the Bragg peak region. This led to a smaller total cross-section and a longer interval between the production of secondary electrons. Plus, high-energy electrons (Equation (23)) would be favored when the proton energy is high. In general, most of the electrons travel a tiny distance before being locally deposited, forming the dense blue area around the central proton line, and hence a high radial dose distribution in the regions with small radii (Figure 4).

Finally, in Table 4, we report the DNA damages in the form of DSBs induced by 0.5 and 0.9 MeV protons. The results from our simulation were compared with those from Nikjoo et al.'s work with the KURBUC model [45]. The difference was found within 10%.

**Figure 5.** A representative track structure for a 5 MeV proton at the entrance part (**a**) and in the Bragg peak region (**b**). The proton was emitted along the positive *Z* direction. Red and blue dots represent the energy depositions by the proton and secondary electrons, respectively. Note: in the two subplots, we kept the same aspect ratio between the z and x/y axes, but plotted them with different ranges.

**Table 4.** The DSB yields (number per Gy per Gbp) obtained under the overlay method for two proton energies.


#### *3.2. Validation of Concurrent Transport*

As for the validation of the concurrent DNA transport module, we first studied the influence of different checking time intervals *ti* and chemical stage durations *tc* on the yields of different radicals. As shown in Figure 6, at a fixed *ti*, the yields of the *eh* and *OH*· radicals reduced when *tc* increased. This was because longer *tc* enabled more reactions among radicals and DNA. For the *eh* and *OH*· radicals, these reactions were mainly consumption channels, resulting in a reduced yield rate when *tc* increased. In contrast, although the presence of DNA also consumed *H*· and *H*2*O*2 radicals, reactions among radicals could contribute positively to the yields of these two radicals. Hence, the production of the *H*· and *H*2*O*2 radicals could be dependent on *tc* in a more complex way. In addition, at a fixed *tc*, varying *ti* from 1 ps to *tc* transformed the simulation from the concurrent method to the overlay method. All lines were connected smoothly, and the G value with *ti* = *tc* matched with that in our previous publication under the overlay method, indicating the self-consistency of the concurrent DNA transport in gMicroMC.

**Figure 6.** The yields of (**a**) *eh*, (**b**) *OH*·, (**c**) *H*·, and (**d**) *H*2*O*2 chemical species at different checking time intervals *ti*and chemical stage durations *tc*.

After examining the self-consistency of the developed concurrent DNA transport method, we comprehensively studied the DSB yields as a function of proton energy *Ek* and chemical stage duration *tc*. The results were compared to Zhu et al.'s work [42] and are shown in Figure 7. From the figure, all data points had relative DSB yields >1, and the R(10) values were larger than R(2.5) for the same *Ek*. This indicated the DSB yields increased when the chemical stage expanded from 1 ns to 10 ns under the concurrent transport frame. The reason was that the longer the chemical stage lasted, the more checks between radicals and DNA were performed, and hence, more DSBs could be formed. Along with the increase of the proton energy, the relative DSB yields exhibited a maximum in the middle energy range. Comparing the data from our simulation to that from Zhu et al.'s work, the trends generally agreed, especially for the R(2.5) data. Some larger discrepancies existed for the R(10) values, which could be explained partially by the different radical diffusion rates and different DNA geometries applied in the two works. For example, the diffusion rate of the *OH*· radical was larger in our package. This could make *OH*· diffuse longer and experience a higher scavenging rate from the histone protein within one diffusion step. In addition, a larger diffusion rate could result in a smaller reaction radius between *OH*· and the DNA sugar-phosphate moiety. Both led to a reduced DSB yield. The longer the *tc* was, the more reduction effect it could create, such that we would obtain a smaller relative DSB yield than that from Zhu et al.'s work for the R(10) data than for the R(2.5) data.

**Figure 7.** The relative DSB yields at different chemical stage durations *tc* and different proton energies with *R*(*t*) = *DSB*(*tc*=*<sup>t</sup>*) *DSB*(*tc*=<sup>1</sup> *ns*). The data from gMicroMC simulation were compared to that from Zhu et al.'s work [42].

## *3.3. Computational Efficiency*

After evaluating the two newly developed features of gMicroMC by comparing to the NIST data or simulation results from other packages, it was important to evaluate the time performance of the new modules for practical applications. In Table 5, we show the simulation time for the physical transport of 1, 10, and 100 protons at 1 and 10 MeV with gMicroMC on a single GPU card (Nvidia V100). As can be seen, it only took 2 and 4 s to transport a single proton with an initial kinetic energy of 1 and 10 MeV. In contrast, it could take multiple hundreds to thousands of seconds to perform similar simulations with single-CPU-based packages, showing the high efficiency of gMicroMC. It is also important to point out that when raising the proton numbers from one to one-hundred, the simulation time only increased by <5 and 10 folds for 1 and 10 MeV protons, respectively. This feature was against the linearly increasing behavior for typical CPU-based simulations, making gMicroMC especially suitable for large-scale simulations. Actually, when the proton number was small, the parallel computing capacity of the GPU was far from being exhausted when launching the kernel for the primary particle transport, such that increasing the proton number, the running time would not significantly increase.


**Table 5.** The simulation time (s) of physically transporting 1, 10, and 100 protons for proton energies at 1 and 10 MeV by gMicroMC on a single GPU card.

As for the simulation time of the concurrent transport (*ti* = 1 ps) in the chemical stage, it could be affected by many parameters, such as the number of radicals, the chemical stage duration, and the DNA complexity. Here, we reported the simulation time for a representative case. Considering that the number of radicals produced from one a 100 keV electron or a 1 MeV proton was roughly 10<sup>5</sup> and the longest chemical stage duration used in the overlay method was ∼1 μs, we set the initial number of radicals *Ni* = 10<sup>5</sup> and chemical stage duration *tc* = 1 μs in our simulation. We also included our DNA model containing of ∼6.2 × 10<sup>9</sup> bps for a human lymphocyte cell nucleus in the simulation, the complexity of which is high. The simulation time was found to be 470 s with gMicroMC on a single GPU. Compared to the simulation time of 31 s under the overlay scheme for gMicroMC, the concurrent transport frame was still quite efficient since the simulation became much more sophisticated with the presence of DNA. This indicated that gMicroMC can be applied

in simulations with realistic settings. In comparison, restrictions on the reaction region and time duration, etc., are typically required in other packages for memory- or time-saving purpose [22,23].
