A study of indirect action’s impact on simulated neutron-induced DNA damage

Objective. The risk of radiobiological stochastic effects associated with neutrons is strongly energy dependent. Recent Monte Carlo studies simulating neutron-irradiated nuclear DNA have demonstrated that this energy dependence is correlated with the relative biological effectiveness (RBE) of neutrons to inflict DNA damage clusters that contain difficult-to-repair double-strand breaks. However, these previous investigations were either limited to modeling direct radiation action or considered the effects of both direct and indirect action together without distinguishing between the two. In this study, we aimed to quantify the influence of indirect action in neutron irradiation scenarios and acquire novel estimations of the energy-dependent neutron RBE for inducing DNA damage clusters due to both direct and indirect action. Approach. We explored the role of indirect action in neutron-induced DNA damage by integrating a validated indirect action model into our existing simulation pipeline. Using this pipeline, we performed track-structure simulations of monoenergetic neutron irradiations (1 eV to 10 MeV) in a nuclear DNA model and analyzed the resulting simple and clustered DNA lesions. We repeated the irradiation simulations for 250 keV x-rays that acted as our reference radiation. Main results. Including indirect action significantly increased the occurrence of DNA lesions. We found that indirect action tends to amplify the damage due to direct action by inducing DNA lesions in the vicinity of directly-induced lesions, resulting in additional and larger damage clusters. Our neutron RBE results are qualitatively similar to but lower in magnitude than the established radiation protection factors and the results of previous similar investigations, due to the greater relative impact of indirect action in photon-induced damage than in neutron-induced damage. Significance. Although our model for neutron-induced DNA damage has some important limitations, our findings suggest that the energy-dependent risk of neutron-induced stochastic effects may not be completely modeled alone by the relative potential of neutrons to inflict clustered lesions via direct and indirect action in DNA damage.


Introduction
The risk of neutron-induced stochastic radiobiological effects, such as radiation-induced cancer, in the human body is known to be strongly energy dependent (ICRP 2007, Sato et al 2013a, US NRC 2021. Monte Carlo studies modeling neutron-irradiated nuclear DNA have demonstrated that this energy dependence is correlated with the relative propensity of neutrons to inflict clusters of DNA lesions containing double-strand breaks (DSBs) (Baiocco et al 2016, Montgomery et al 2021. However, these investigations either only modeled direct radiation action on DNA substructures (Montgomery et al 2021) or modeled both direct and indirect action, but did not consider their individual effects (Baiocco et al 2016). Therefore, a study modeling and quantifying the potentially influential effects of indirect action on the energy dependence of neutron-induced stochastic effects has been outstanding, and is the subject of this report. Humans are exposed to neutron radiation in a variety of scenarios including space travel (Benton et al 2001, Koshiishi et al 2007, nuclear incidents (Shuryak et al 2020), and radiotherapy procedures such as fast neutron therapy (Jones 2020), neutron capture therapy (Sauerwein et al 2012), and high-energy (8 MeV) radiotherapy (Howell et al 2006, Maglieri et al 2015. Cancer patients treated with high-energy radiotherapy are exposed to a (relatively) low absorbed dose of various types of non-therapeutic 'secondary' radiation, including a polyenergetic spectrum of neutrons that are generated by interactions of the high-energy primary radiation with matter. Because of this unavoidable exposure to secondary neutrons, radiotherapy patients, especially pediatric patients (Friedman et al 2010), are at risk of developing iatrogenic second cancers later in life. What makes these secondary neutrons particularly pernicious is (1) their whole-body area of effect (AoE, a term we have borrowed from the world of video games that refers to the area over which an action, in this case radiation damage, can be expected to have an effect (Heger et al 2009)) and (2) their elevated and energy-dependent risk for inducing stochastic effects.
Although a patient's absorbed dose from secondary neutrons is only a fraction of their therapeutic dose in the case of high-energy photon radiotherapy (Kry et al 2017), the former's AoE is the whole body, whereas the latter's is the irradiated volume containing the malignancy. Given the widely-accepted radiation protection paradigm that posits that stochastic radiobiological effects can occur at any dose, AoE may play a more important role in iatrogenic neutron-induced carcinogenesis than absorbed dose. Furthermore, compared to other types of radiation, neutrons are associated with a risk for inducing stochastic effects that is relatively higher in magnitude and energy dependent. This elevated and energy-dependent risk is typically quantified by the neutron weighting factors (w R ) (ICRP 2007) and the neutron quality factors (Q) (US NRC 2021). These independently-determined factors are widely-used radiation protection quantities for characterizing radiationrelated biological risk and were derived from aggregated data on the relative biological effectiveness (RBE) of neutrons for inducing stochastic effects. Thus, higher w R and Q values denote higher risk of stochastic effects. It is important to recognize, however, that although essential in radiation protection, these quantities were not intended for estimating carcinogenic risk. The energy dependence of neutron-induced stochastic effects is especially relevant during high-energy radiotherapy because the most abundant secondary neutrons permeating the treatment room (fast neutrons that peak at energies around 1MeV (Howell et al 2006, Maglieri et al 2015) have the highest w R and Q values (ICRP 2007, US NRC 2021.
The neutron weighting and quality factors w R and Q follow similar trends across energy, however, there is a noticeable discrepancy between their magnitudes as depicted in figure 1. This discrepancy can be attributed to the wide variety of neutron RBE data that were used to derive the weighting and quality factors: findings from epidemiological studies and radiobiological experiments investigating different stochastic biological endpoints such as the induction of dicentric chromosomal aberrations in cells, and carcinogenesis and life shortening in animals (ICRP 2007, US NRC 2021. The difference between the curves and the need to better understand the biophysics behind their energy dependence, have motivated parallel research efforts such as that of the ANDANTE group (Ottolenghi et   endpoint that may be modeled. Both groups used Monte Carlo simulations to model neutrons impinging on nuclear DNA and have demonstrated that the energy dependence of the w R and Q factors can be correlated with the RBE of neutrons to inflict DSB-containing DNA damage clusters. The interest in radiation-induced nuclear DNA damage stems from the most widespread mechanistic model of carcinogenesis, the more-than-a-century-old somatic mutation theory (Boveri 1914, Hanahan and Weinberg 2000, 2011, Vaux 2011, which posits that genetic alterations or mutations can lead to carcinogenesis. DNA damage clusters are pertinent in the radiation context because they rarely occur via endogenous processes but can be effectively induced by ionizing radiation. Furthermore, their associated damage repair mechanisms (especially for clusters containing DSBs) are prone to misrepairs that can result in genomic mutations. Although such DNA damage may simply result in cell death or cell cycle arrest, persistent mutations in specific gene segments of inflicted cells that survive and divide may cause them to develop traits characteristic to cancer cells. As such, many theorize that the main mechanism of radiation-induced mutagenesis (and eventual carcinogenesis) is via the infliction of DNA damage clusters by direct or indirect physico-chemico action (Goodhead 1994, Ward 1995, Magnander and Elmroth 2012, Georgakilas et al 2013, Sage and Shikazono 2017. It is important to realize that physico-chemico damage to nuclear DNA is considered to be but one pathway for ionizing radiation to induce mutagenesis. There is evidence that mutations can also arise from non-DNAtargeted effects of ionizing radiation such as bystander effects and genomic instability (Iyer andLehnert 2000, Little 2000). Baiocco et al (2016) from the ANDANTE group used the simulation frameworks PHITS (Sato et al 2013b) and PARTRAC (Friedland et al 2011) to study the neutron-induced emergence of DNA damage clusters shorter than 30 bp, with a particular focus on those containing at least two DSBs within 25 bp (sometimes referred to as DSB++ lesions (Friedland et al 2011)) due to both direct and indirect action. On the other hand, Montgomery et al (2021) from our group used TOPAS (Perl et al 2012) and TOPAS-nBio (Schuemann et al 2019) to investigate an extended scope of DNA damage clusters: those with and without DSBs (respectively referred to as complex DSB clusters and non-DSB clusters) where complex DSB clusters serve as the superset to DSB++ lesions. Thus, although the neutron RBE estimates of Baiocco et al (2016) and Montgomery et al (2021) exhibit similar qualitative dependence in energy, their magnitudes differ considerably (see figure 1). Importantly, although Montgomery et al (2021) used a more inclusive definition of clustered lesions that are believed to be heavily implicated in mutagenesis (Magnander and Elmroth 2012, Georgakilas et al 2013, Sage and Shikazono 2017, they did not model the indirect action of ionizing radiation that occurs promptly (∼1 ps (Le Caër 2011)) after direct action. Indeed, various in vitro (Roots and Okada 1972, 1975, Ito et al 2006 and in silico (Nikjoo et al 2002, Lampe et al 2018, Sakata et al 2020 studies have demonstrated that indirect action is the dominant damageinducing mechanism of ionizing radiation to the DNA molecule in the low-LET domain, while direct action dominates at higher LET (Hirayama et al 2009).
Results of condensed-history neutron simulations inside a human tissue phantom performed separately by Baiocco et al (2016) and our group (Lund et al 2020) have shown that neutrons interacting with the human body produce a polyenergetic fluence of multiple secondary particle types, many of which fall into the low-LET domain. Moreover, simulations modeling both direct and indirect action of various radiation types have reported that hybrid DSBs (having simultaneously a direct and an indirect strand break or SB) constitute a significant portion of all induced DSBs (de la Fuente Rosales et al 2018, Zhao et al 2020, Zhu et al 2020 and that the yield of complex DSB clusters increased when indirect action was included (Henthorn et al 2019). Thus, to achieve a more complete in silico model of neutron-induced nuclear DNA damage and, by extension, a more accurate assessment of the energy-dependent neutron RBE for inducing complex DSB clusters, indirect radiation action with DNA must be incorporated into the model.
In this work, we extended our existing simulation pipeline (Lund et al 2020, Montgomery et al2021) to include a model for indirect action. We validated our implemented indirect action model by benchmarking against published results and performed neutron simulations using the updated pipeline. Our goal was to elucidate and quantify the effects of indirect action on the quantity of neutron-induced DNA damage yields and on the properties of resulting damage clusters. We also aimed to recalculate our group's previous energydependent neutron RBE estimates related to the induction of DNA damage clusters (Montgomery et al 2021), but this time considering both direct and indirect action.

Overview
Our methodology for estimating the energy-dependent biological risk associated with neutrons consisted of two parts. In part 1 (Modeling and validation), we adopted and adapted the model for the indirect action of ionizing radiation by Zhu et al (2020) to extend our group's existing simulation pipeline (Montgomery et al 2021). We then tested the validity of this indirect action model by simulating proton irradiations on our DNA model (Montgomery et al 2021) and benchmarking the resulting DNA damage yields against published results from experimental and computational studies. In part 2 (Neutron simulations), we repeated the simulated neutron and photon irradiations of our previous work (Montgomery et al 2021), but this time including our validated model for indirect action (from part 1) in the simulation pipeline. The resulting DNA damage yields from these irradiations were used to estimate neutron RBEs for DNA damage due to direct and indirect action. An overview of the simulations and analyses performed in parts 1 and 2 of our work is presented in figure 2. Both parts are described in detail below.
2.2. Part 1: modeling and validation 2.2.1. DNA model As mentioned earlier, our previous work that estimated neutron RBE related to stochastic effects was limited to modeling direct action. For this current study, we used the same human fibroblast cell model that we developed in the previous one, the details of which are summarized in figure 3 and described in detail in Montgomery et al (2021). Consistent with the indirect action model of Zhu et al (2020), all products of water radiolysis available in the chemistry constructor (see below) were simulated in our work. Reactive species originating from inside DNA and histone volumes were immediately terminated and histones were considered as scavengers of hydroxyl (·OH), solvated electron (eaq ), and hydrogen (H·) radicals. Similarly, our simulated ·OH radical tracks were terminated after an indirect action event with a sugar-phosphate or nucleobase volume, because the ·OH radicals would have been rendered stable by the interaction. In reality, ·OH radicals interact with the sugar- phosphate backbone via the abstraction (i.e., removal) of a deoxyribose H-atom that then combines with the attacking radical (·OH) to produce a stable compound (H 2 O) (Balasubramanian et al 1998). As for nucleobases, ·OH radicals enter into an addition reaction with nucleobase C-H bonds where they effectively become absorbed (Talpaert-Borlè 1987, Balasubramanian et al 1998).

Indirect action model
In our simulations (part 1 and part 2), we used the G4EmDNAPhysics_hybrid2and4 physics constructor (a combination of the opt2 and opt4 constructors of Geant4) used by Montgomery et al (2021)

Types of DNA damage
The types of DNA damage investigated in this study were the same as those by Montgomery et al (2021), but further stratified according to their cause. That is to say, simple lesions, namely single-strand breaks (SSBs) and base lesions, could be direct or indirect, whereas DSBs, complex DSB clusters, and non-DSB clusters could be direct, indirect, or hybrid.  of 40%. We opted to limit the DNA damage infliction to ·OH radicals due to their significantly higher reactivity compared to other species (Roots and Okada 1975) and to facilitate the validation of our pipeline-integrated indirect action model through comparison with reports in the literature that followed the same approach.
The 40% damage probability was based on findings showing that two out of the five C-H bonds in deoxyribose molecules found in the sugar-phosphate backbone of the DNA are more topologically accessible to ·OH radical attacks, and are thus more prone to damage (Balasubramanian et al 1998). Given that backbone molecules constitute four out of six of the spherical volumes in the nucleotide base pairs in our DNA model, we estimated that about 27% of all ·OH interactions with DNA volumes in our simulations would result in a SB. This value lies between experimental estimates ranging from 14%-22% (Van Rijn et al 1985) and 32%-44% (Milligan et al 1993). In reality, about 80%-90% of all ·OH interactions with the DNA occur on nucleobases, while about 10%-20% occur on the sugar-phosphate backbone (O'Neill 2001, Greenberg 2016). However, due to the lack of published data for simulated base lesions due solely to indirect action that can be used as a baseline, we decided to use the same damage probability of 40% for simulated ·OH-nucleobase interactions, analogous to how our direct action implementation uses the same energy-deposition threshold of 17.5 eV for both the deoxyribose and nucleobase volumes (Montgomery et al 2021).

Model validation
To test the validity of our indirect action model, we simulated the irradiation of DNA by monoenergetic protons and compared our findings with published computational and experimental data. This exercise is relevant for two reasons. Firstly, protons are the dominant dose contributors for primary neutrons around 1 MeV (Lund et al 2020) where the neutron RBE peaks as observed in figure 1. Secondly, we wanted to determine the compatibility of our indirect action model with open-source physics and chemistry models in producing realistic DNA damage counts under comparable conditions that have been previously published by other groups.
In our irradiation scenario identical to that of Zhu et al (2020), the protons were generated randomly throughout the surface of the cell nucleus and directed inwards. The irradiation was terminated after 1 Gy of dose was delivered to the nucleus. This procedure was performed for five different proton energies ranging from 500 keV to 20 MeV.
Using our existing direct action model, our new indirect action model, and our updated clustering algorithm, we were able to extract yields for SSBs (Y SSB ; direct and indirect) and DSBs (Y DSB ; direct, indirect, and hybrid). The extracted damage yields were normalized to the dose and to the total number of base pairs in our DNA model (6.3 Gbp). The mean values obtained across 100 independent runs were then benchmarked against published data from simulated proton irradiations (Nikjoo et al 2001, Friedland et al 2003, 2011 (Roots and Okada 1972, Chapman et al 1973, Roots and Okada 1975, Frankenberg et al 1999, Belli et al 2000, Campa et al 2005. The standard uncertainty of the mean for each set of runs was also calculated. Further clustering beyond simple DSBs was not performed for this validation analysis, as there were no published results available for comparison. 2.3. Part 2: neutron simulations 2.3.1. Neutron simulation pipeline After validating our model of indirect radiation action, we simulated monoenergetic neutron irradiation of DNA using our updated simulation pipeline (blue block in figure 2). Previously, our group (Lund et al 2020) used Geant4 to determine the secondary particle spectra and relative dose contributions of secondary particles produced by neutrons in a human tissue-equivalent phantom, the ICRU-4 sphere (White et al 1989). This analysis was performed for three regions of increasing depth, as shown in figure 4. In the present work, similar to that of Montgomery et al (2021), we sampled these secondary particle spectra and used them as input to simulate the irradiation of our nuclear DNA model using TOPAS-nBio. This two-step approach, which was a feature of our previous work (Lund et al2020), served to allow us to reuse our previously-generated neutron secondary particle spectra and also to simplify the comparison of our RBE results with our previous direct-action only estimates that used the same process (Montgomery et al 2021).
Details regarding the process of obtaining these secondary charged particle spectra are discussed extensively elsewhere (Lund et al 2020, Montgomery et al 2021. In brief, 1 × 10 10 monoenergetic primary particles of uniform fluence were simulated to generate secondary particle species whose tracks were immediately terminated upon entering a scoring volume. However, due to limitations of Geant4-DNA transport models (Incerti et al 2018), electrons generated above 1 MeV were instead allowed to propagate and tracked down to 1 MeV (Lund et al 2020). The subsequent tracks produced during this slowing down process were included in the resulting charged particle spectra.
The sampled secondary particles were generated randomly throughout the entire volume of our full cell model (see figure 3) with isotropic direction and their direct and indirect action with the DNA volumes were modeled using track-structure simulations. For each simulated irradiation of monoenergetic neutrons, the total target dose to the nucleus, accounting for all species of secondary particles, was set to 1 Gy. To achieve this, the number of simulated particles for each secondary species was limited by the known relative dose contributions of each species as determined by Lund et al (2020). Note that the secondary particle species we considered in our simulations were limited to electrons, protons, and alpha particles due to the lack of complete physics models in Geant4-DNA for the transport of heavier ions at our energy range of interest (Incerti et al 2016).
The same procedure of particle sampling and irradiation was also performed for the 250 keV x-ray photons that served as our reference radiation for RBE determination. The yields for each type of damage discussed in section 2.2.3 and the properties of the clusters formed were extracted for data analysis and for comparison with our previous direct-action-only results (Montgomery et al 2021).  In these equations, the subscript n 0 refers to neutrons, while X refers to x-ray photons. The DNA damage yields induced by the reference 250 keV x-rays did not need to be summed across i since the only associated secondary particle species was electrons. Note that lesions inflicted by tertiary particles (i.e., particles generated by the initial secondary charged particles during a simulation run) were considered to be inflicted by the secondary species i in our calculations. Similar to our damage yields in part 1, our reported values for part 2 were normalized to the total number of base pairs and dose delivered, and were calculated as the mean of 100 independent runs. The standard uncertainty for each mean was also calculated and reported. The corresponding neutron RBE values were then calculated using the following formula: (1) cluster length, (2) cluster complexity, and (3) cluster density. Cluster length refers to the number of base pairs between the two most extreme lesions in each cluster, while cluster complexity is the total number of SSBs, base lesions, and DSBs (if any) per cluster. Lastly, cluster density is simply the cluster complexity divided by cluster length.

Results and discussion
3.1. Part 1: modeling and validation 3.1.1. Model validation As mentioned in section 2.2.5, we have benchmarked our results against various other simulation data. However, given the differences in the many simulation parameters and irradiation scenarios with competing effects on the resulting DNA damage yields, we opted to primarily compare our model to the most similar one, which is that of Zhu et al (2020). The complete comparison we have performed for part 1 can be found in appendix A.1.
In figure 5(a), which presents direct, indirect, and total + Y p SB ( damage due to protons), we see that our direct + Y p SB are similar to that of Zhu et al (2020), but our indirect and total + Y p SB are relatively lower (these observations can be extended to + Y p SSB ). In principle, the only differences between our proton irradiation simulations are the DNA model (size(s) and spatial distribution(s) of sensitive volumes), the physics constructor, and the scavenging of ·OH radicals by DNA and histone volumes regardless of damage induction.
The effects of DNA model and physics constructor in the context of direct action have been discussed before by Montgomery et al (2021). Due to the nature of sensitive volume distribution, damage events (direct and indirect) are more likely to occur in DNA models that are denser and contain more sensitive volume. The nuclear DNA model used by Zhu et al (2020) has an average density of 14.4 Mbp/μm 3 and sensitive DNA volume proportion of 1.5% (calculated from their nucleotide base pair volume of about 1.05 nm 3 ). Ours is about 8% less dense at 13.3 Mbp/μm 3 and has about 45% less proportion of sensitive DNA volume at only 0.8% (each base pair is about 0.6 nm 3 ). The fact that our total and indirect + Y p SB in figure 5(a) is only about 30% lower than theirs (instead of 45%) is likely due to the competing effect of larger histone sizes of about 195 nm 3 in the model of Zhu et al (2020) (ours: 105 nm 3 ) that can terminate more simulated ·OH tracks. Thus, overall, ·OH molecules are more likely to induce damage (leading to larger + Y p SB ) in the model of Zhu et al (2020) compared to ours. Furthermore in figure 5(a), our calculated proportion of indirect (green) relative to total (black) + Y p SB is approximately 68% in the low-LET domain. This value is in agreement with other published computational estimates (see figure A1(a)) that range from 62%-85%. On the other hand, experimental estimates for the percentage of SBs induced by radiolytic species with the DNA via indirect action vary from 65%-71% Okada 1972, 1975), with an estimated 62% arising from ·OH radicals specifically (Chapman et al 1973) and 3%-9% from other species. Thus, our implemented indirect action model produces yields of indirect damage that are also consistent with published experimental results. Figure 5(b) presents a stratification of our + Y p DSB according to damage cause as a function of LET. Whereas + Y p SB decreases with increasing particle LET (as seen in figure 5(a)), + Y p DSB increases with LET. Considered together, these trends indicate the expected clustering of SSBs at higher LET values. Comparing our results with that of Zhu et al (2020), we see that although our total + Y p DSB values are again on the lower end, our hybrid + Y p DSB values are consistent with theirs. The dominance of hybrid DSBs (red) compared to their strictly direct (blue) and strictly indirect (green) counterparts demonstrates the importance of modeling the synergetic effects of direct and indirect action when studying radiation-induced DNA damage.
Although our damage yields are consistently lower than that of Zhu et al (2020), we found that our + Y p DSB values are in strong agreement with most other simulations and are comparable with experimental data (see figure A2(a)). This result is important given that we are interested in the formation of complex DSB clusters, which are believed to play a key role in the emergence of radiation-induced mutations. Moreover, we also found that our + Y p SSB -to-+ Y p DSB ratios are comparable with many published computational data (see figure A2(b)). Given that we have reasonable amounts of + Y p DSB , our + Y p SSB -to-+ Y p DSB ratios suggests that our + Y p SSB are also reasonable. In brief, our simulated proton irradiation results in part 1 show that our indirect action model that has ·OH radicals as its sole damage-inducing chemical species with the damage probability of 40% in our DNA model, together with our direct action model, are compatible with open-source and validated physics and chemistry models in producing realistic + Y p DSB with realistic indirect-to-direct damage ratios. Therefore, we consider our simulation pipeline to be appropriate.

Open-source code release
In support of open science, the code of our updated simulation pipeline for radiation-induced DNA damage has been released open-source as an update and extension to our previous version of the TOPAS-nBio application TOPAS_Clustered_DNA_Damage (Montgomery et al 2022). Aside from minor changes in the nuclear DNA model (implemented as a custom TOPAS geometry component) to uniquely identify histone volumes, the main modifications to our previous simulation code involved our DNA damage scorer (implemented as a custom TOPAS scorer) and are summarized as follows: • Added: simulation of indirect action events and indirect damage scoring, and user-modifiable simulation parameters related to indirect action.
• Updated: DNA damage clustering algorithm to account for indirect and hybrid lesions.
The documentation accompanying our previous release has been updated and example parameter files related to our indirect action simulations have been made available.

Part 2: neutron simulations
In this section, only a subset of our neutron and photon simulation results is presented and discussed for maximum clarity and concision. Specifically, we used our calculated values in the intermediate depth (see figure 4) to representatively illustrate our most significant findings. This choice of radial depth is consistent with our previous paper (Montgomery et al 2021). Some results from the inner and the outer scoring regions are also presented for comparison.

DNA damage yields
The increase in neutron damage yields obtained by incorporating indirect action into our model is shown in figure 6(a) (for DSBs and simple lesions) and figure 6(b) (for clustered lesions). The largest absolute increase can be observed in Y n SSB 0 (red), followed by Y n BL 0 (cyan) and -Y n N DSB 0 (blue), for initiating neutron energies less than ∼100 keV, where the dominant dose contributor is electrons (Lund et al 2020). For initiating neutron energies higher than ∼100 keV, protons (having higher LET than electrons) become the dominant dose contributor (Lund et al 2020), which leads to an increase in clustered lesions (black and blue in figure 6(b)) and a decrease in simple lesions (red and cyan in figure 6(a)). On the other hand, across the full range of initiating neutron energies investigated, we remark a relative constancy of Y n DSB 0 (green in figure 6(a)). Considering these observations together with the secondary particle dose contributions reported in Lund et al (2020), we infer that, starting around 100 keV, the secondary protons inflict simple lesions that are close enough to form clusters, and thus many of the DSBs induced in this manner become complex DSB clusters.
Although similar observations can be made in the case of direct action alone (Montgomery et al 2021) (broken lines), the inclusion of indirect action serves to accentuate and amplify these trends. Indeed, separating the neutron-induced clustered DNA lesions according to their damage cause reveals that the proportion of hybrid clusters for both neutrons and reference 250 keV x-rays is approximately 80% for complex DSB clusters (see figure A3(a)) and about 50% for non-DSB clusters (see figure A3(b)). The abundance of hybrid DNA damage clusters suggests that indirect action increases the damage potential of direct action by inflicting additional lesions near direct damage sites.
The ratio of Y n SSB 0 to Y n BL 0 in figure 6(a) is approximately 2:1 when considering direct action alone (broken red and cyan lines). This is simply due to the ratio of sugar-phosphate backbone molecules to nucleobases in the DNA model (Montgomery et al 2021). On the same figure, with the inclusion of indirect action in the simulations, we observe a larger increase in Y n SSB 0 relative to Y n BL 0 (solid red and cyan lines) despite the identical ·OH damage probability (40%) towards backbone and base volumes. This effect is likely because the primary particles that induce direct damage can traverse and deposit energy across multiple DNA volumes, whereas each simulated ·OH radical can only interact once before being terminated, whether or not it is able to induce indirect damage. Thus, the observed preference for SSBs over base lesions is likely due to the geometry of the DNA model where backbone volumes are more exposed to radical attacks. This result may be unrealistic given that the kinetically-preferred ·OH interaction with DNA occurs at the nucleobases, as discussed in section 2.2.4. Thus, our Y n BL 0 in figure 6(a) (cyan) and consequently, -Y n C DSB 0 and -Y n N DSB 0 in figure 6(b) (black and blue), as well as their x-ray-induced counterparts, are likely underestimates.

Neutron RBE
In figure 7(a), our estimates of neutron RBE for various types of DNA damage that incorporate indirect action (solid lines) are lower than our previous estimates that considered direct action alone (broken lines). The most significant changes can be observed for the aggregate damage types: RBE C−DSB (black), RBE N−DSB (blue), and RBE DSB (green). This reduction in RBE was also observed across various scoring depths (see figure 4), as shown representatively for complex DSB clusters in figure 7(b). For this representative case, our current estimates for RBE C−DSB across penetration depth converge to a lower RBE value of about 3.5 at initiating neutron energies above 1 MeV. Nevertheless, from these two figures, we see that the qualitative trends in neutron RBE from our previous direct-action-only results are conserved.
A closer look at our yields of clustered lesions in figure 6(b) reveal that the relative increase in photoninduced yields (star to circle markers) is considerably higher than that of neutron-induced yields (broken to solid lines) when indirect action is considered. Remarkably, for complex DSB clusters, the relative increase in -Y X C DSB is up to two times larger than that of -Y n C DSB 0 . Thus, it is not surprising that our current RBE C−DSB estimates in figure 7(b) (solid lines) fall to about half of our previous estimates (broken lines) for neutron energies between 100 keV and 10 MeV.
The fact that the relative increase in Y C−DSB and Y N−DSB upon the inclusion of indirect action in figure 6(b) (broken to solid lines and star to circle markers) is considerably higher for lower-LET particles (i.e., reference photons and neutrons with energies below ∼100 keV as discussed in section 3.2.1) than for higher-LET particles (i.e., neutrons with energies above ∼100 keV) supports the hypothesis that indirect action serves to amplify direct damage by inducing DNA lesions close to direct damage sites. This damage amplification by indirect action increases the likelihood of damage grouping in the lower-LET case where direct lesions are more likely to be isolated, thereby increasing Y DSB , Y N−DSB , and Y C−DSB for photons and lower-energy neutrons (see figure 6). Although this effect is also present in the higher-LET case, it has less of an impact on the number of damage clusters since direct lesions induced in this case are already likely to form clusters (broken lines in figure 6(b)). Moreover, this damage amplification acts to increase the length and complexity of direct DNA damage clusters (see section 3.2.3).
Comparing our neutron RBE values at the intermediate scoring volume with published factors of neutron stochastic risk in figure 7(c), we see that our neutron RBE C−DSB (solid black line) peaks at a similar initiating neutron energy of around 1 MeV as the others, but is considerably lower in magnitude. Overall, our results obtained with a more comprehensive model of neutron-induced DNA damage qualitatively suggest that the capacity of neutrons to induce complex DSB clusters may indeed be correlated with published factors that describe stochastic effects in humans after exposure to neutron radiation. However, the quantitative differences in these data suggest that clusters having at least two DSBs (red) correlate more with radiobiological risk compared with clusters with at least one DSB (black and blue). Furthermore, such differences are indicative that there may be additional as-yet-unmodeled factors at play in neutron-induced carcinogenesis such as DNA damage repair mechanisms and the non-DNA-targeted effects of ionizing radiation mentioned in section 1.

Cluster properties
By including indirect action, we see a substantial increase in the length of clustered DNA damages in figure 8(a) and a minor increase in their complexity (lesion count per cluster) in figure 8(b) relative to the results obtained for direct action alone (broken to solid lines and star to circle markers). This finding further confirms the damage-enhancing property of indirect action discussed in sections 3.2.1 and 3.2.2.
Moreover, in figures 8(a) and (b), the impact of indirect action is indicated by the formation of a shallow peak (solid lines) in the lower energy region where lower-LET electrons contribute most of the dose. We remark the familiar increase of these cluster properties at around 100 keV where protons start to dominate the neutron dose, and another one for complex DSB clusters above 2 MeV with the onset of alpha particles as substantial dose contributors (Lund et al 2020). Given that the increase in cluster length is greater than the increase in cluster complexity when indirect action is incorporated, we observe an expected decrease in cluster density in figure 8(c) (broken to solid lines and star to circle markers), with the concavity of the trends inverted. Comparing the two types of clustered DNA lesions, we find that on average when indirect action is involved, complex DSB clusters are 40%-50% longer (solid lines in figure 8(a)) and contain 55%-80% more lesions than non-DSB clusters (solid lines in figure 8(b)). These values were previously estimated by our group to be 30%-50% and 40%-70% respectively with direct action alone (Montgomery et al 2021). Thus, indirect action has a comparable effect on both complex DSB clusters and non-DSB clusters. At least three simple lesions are required to form a complex DSB cluster, while a non-DSB cluster needs at least two. From our results in figures 8(a) and (b), we find that on average, neutron-induced complex DSB clusters are about 15 bp long with most having 3-5 lesions (black solid lines), while non-DSB clusters are about 10 bp long and contain 2-3 lesions (blue solid lines). With a total base pair count of about 6.3 Gbp in our model (Montgomery et al 2021), the chance of occurrence of overlap between two independent and randomly distributed clusters in our nuclear DNA model is negligible. We can also observe from the plots in figure 8 that the cluster properties of complex DSB (black lines) and non-DSB (blue lines) clusters follow similar trends across the neutron energy range we investigated. However, in relative terms, cluster length varies more with energy than does cluster complexity.
With increasing depth in the human tissue-equivalent phantom, we found that the local maximum at the lower LET range we have remarked earlier decreases and shifts to higher initiating neutron energy (see figure  A4). This trend is comparable to how the neutron RBE curves in figure 7(b) behave. Lastly, at initiating neutron energies above 500 keV, both cluster length and complexity tend to converge and become independent of depth of measurement (see figure A4), which again can be observed in the neutron RBE curves of figure 7(b).
It is likely that a true radiobiological RBE is dependent on these cluster parameters but a theoretical formalism for incorporating these parameters was beyond the scope of this work. Doing so will require DNA repair modeling that is influenced by these parameters. As it stands, we present information on these parameters to indicate that damage yields alone do not paint the full picture.

Limitations and future work
The physics-related and DNA-model-related limitations of our simulation pipeline have been discussed in our previous publication (Montgomery et al 2021). However, we would like to add that our DNA model represents the nucleus as a whole, and thus, we cannot relate the damage yields we obtained in our simulations to specific genes in specific chromosomes that might lead to specific mutations with known outcomes such as cell death or certain carcinogenesis. As for additional chemistry-related limitations, our simulations did not model the DNA damage induced by non-·OH species which is estimated to contribute to about 3%-9% of all SBs induced by low-LET radiation (see section 3.1.1). For future work with new available data, our updated pipeline allows for a convenient modification of the damage probabilities of other simulated chemical species interacting with DNA volumes.
Our simulations also did not take into account the process by which base lesions can turn into SBs: a base lesion in the form of a radical site induced in a nucleobase (via direct or indirect action) may lead to a SB if the radical site transfers to the sugar-phosphate backbone via H-atom abstraction from nearby deoxyribose C-H bonds (Dizdaroglu andJaruga 2012, Greenberg 2016). This effect will be important to consider in future work if the probability of base lesion induction by ·OH radicals (and other radiolytic species) is to be increased relative to that of SB induction.
Also, in TOPAS-nBio, an extended chemistry constructor is already available (TsEmDNAChemistryExtended) that includes a model for four additional oxygen species (O 2 , -O 2 , HO 2 , and HO 2 − ) compared to TsEmDNAChemistry that was used in our work (Schuemann et al 2019). Oxygen species are of interest because their presence has been demonstrated to influence the damage potency of indirect radiation action (Chapman et al 1973), often measured via the parameter known as the oxygen enhancement ratio (OER). For future studies, using TsEmDNAChemistryExtended would allow us to analyze the impact of the presence of these oxygen species in the yield of indirect DNA damage and resulting neutron RBE. The chemistry models available in TOPAS-nBio assume that the medium is made of pure liquid water at neutral pH and 25°C (Ramos-Méndez et al 2018). In reality, however, cells are not made of pure water, intracellular pH ranges from 4.5 to 8.0 (Asokan andCho 2002, Proksch 2018), and normal human body temperature is around 36°C-37°C (Geneva et al 2019). At the time of our study, there were no chemistry models in Geant4-DNA or TOPAS-nBio accounting for these factors, thus their influence on our DNA damage and RBE results are difficult to ascertain. However, once such models become available, the design of our simulation pipeline is robust enough to facilitate their incorporation.
As mentioned earlier, our simulations were limited to prompt DNA damage via direct and indirect action. However, mutagenesis and subsequent carcinogenesis occur much later. Thus, our next steps will be to simulate DNA damage repair mechanisms and investigate the emergence of neutron-induced genomic mutations and how they affect neutron RBE estimates.

Conclusion
In this work we used Monte Carlo simulations of radiation-induced DNA damage in order to carry out a novel investigation into the role of indirect action on the energy dependence of neutron RBE for stochastic effects via the induction of complex DSB clusters. To do so, we integrated a model of indirect action into our previouslypublished simulation pipeline that considered only direct effects, and validated it for protons against other published models and experimental data. Our simulated neutron irradiations of a nuclear DNA model revealed that including indirect action significantly increased DNA damage yields (both simple and clustered lesions), as well as the average length and complexity of the damage clusters thus formed. Furthermore, we found that the majority of clusters contain lesions from both direct and indirect action (i.e., hybrid clusters), demonstrating how indirect action serves to amplify the damage potential of direct action.
We found that our updated estimates of neutron RBE for inducing clustered DNA damage, via both direct and indirect action, exhibit a qualitatively similar energy-dependence to both our previous direct-action-only estimates and published neutron radiation protection factors, but that the magnitudes are much lower. This reduction in magnitude of RBE may be attributed to the fact that indirect action has a greater relative impact for our reference photons than for neutrons. However, to better understand these quantitative discrepancies, a number of outstanding limitations regarding the modeling of DNA damage must be addressed. Beyond this, DNA repair and non-targeted effects should be considered to ascertain a more complete picture of neutron-induced carcinogenic effects.