Track structure Monte Carlo simulations are a useful tool to investigate the damage induced to DNA by ionizing radiation. These simulations usually rely on simplified geometrical representations of the DNA subcomponents. DNA damage is determined by the physical and physicochemical processes occurring within these volumes. In particular, damage to the DNA backbone is generally assumed to result in strand breaks. DNA damage can be categorized as direct (ionization of an atom part of the DNA molecule) or indirect (damage from reactive chemical species following water radiolysis). We also consider quasi-direct effects, i.e., damage originated by charge transfers after ionization of the hydration shell surrounding the DNA. DNA geometries are needed to account for the damage induced by ionizing radiation, and different geometry models can be used for speed or accuracy reasons. In this work, we use the Monte Carlo track structure tool TOPAS-nBio, built on top of Geant4-DNA, for simulation at the nanometer scale to evaluate differences among three DNA geometrical models in an entire cell nucleus, including a sphere/spheroid model specifically designed for this work. In addition to strand breaks, we explicitly consider the direct, quasi-direct, and indirect damage induced to DNA base moieties. We use results from the literature to determine the best values for the relevant parameters. For example, the proportion of hydroxyl radical reactions between base moieties was 80%, and between backbone, moieties was 20%, the proportion of radical attacks leading to a strand break was 11%, and the expected ratio of base damages and strand breaks was 2.5–3. Our results show that failure to update parameters for new geometric models can lead to significant differences in predicted damage yields.
INTRODUCTION
DNA damage is the dominant mechanism by which ionizing radiation drives biological responses, with cells' fates determined by whether or not they can effectively repair these alterations (1). For example, Cornforth and Bedford observed correspondence between fibroblast cell death and asymmetrical chromosome aberrations (2), and Ballarini and Carante implemented this idea in the BIANCA biophysical model to predict the biological effectiveness of different hadrontherapy beams (3, 4). Detailed descriptions of the processes involved in the interactions between ionizing radiation and genetic material are key to advancing radiation therapy. In this sense, mathematical modeling can help to understand the radiobiology of these processes. Besides analytical approaches to modeling physical and biological effects (5–7), the Monte Carlo method is a powerful tool to characterize radiation effects due to the large number of interactions each particle undergoes (8). Monte Carlo track-structure (MCTS) simulation aims to fully characterize the interaction of ionizing radiation at nanoscopic scales (i.e., at DNA scales) using a full description of the ionizing-particle track structure (9). Algorithms were developed over three decades ago to study the energy imparted by physical processes in microscopic sites (10, 11) and the chemical species generated after radiolysis of water molecules (12, 13) by different ionizing particles. However, uncertainty remains in the translation of microscopic energy deposition to nanoscopic damage to the structure of DNA.
To translate the simulated properties of ionizing radiation into damage inflicted to the DNA, the DNA structure is typically superimposed on charged-particle tracks precalculated in liquid water. This method requires geometrical models of DNA. Three different approaches can be considered with increasing complexity (14): (a) DNA strands as linear segments or cylinders randomly placed in the nucleus volume (linear segment model); (b) DNA nucleotides and subcomponents as volumes with a given shape and size (volume model); and (c) DNA nucleotides decomposed into their atoms (atomistic model). Codes in which damage to DNA components is probabilistically assigned and clustered may be considered as a sub-category of the linear segment approach, such as the Monte Carlo damage simulation (MCDS) (15) or density-based spatial clustering of applications with noise (DBSCAN) (16). While these codes allow efficient simulations (17), they lack detailed representations of the processes involving different parts of the DNA molecule, as well as potential differences caused by changing the compaction of the chromatin. Atomistic models are on the opposite end: theoretically, atom-wise processes could be resolved, although the computational burden would considerably increase. For instance, DNAFabric is an attempt to provide atomistic geometries, allowing the definition of the volume and position of each constituent of a DNA nucleotide based on the Glactone project data (18). Nonetheless, the lack of detailed cross-sections for radiation and atomistic components interactions may challenge the applicability of atomistic models.
Volume models represent a compromise between the very fast but too simplistic linear models and the very detailed but computationally expensive atomistic models. Volume models allow for the multi-scale definition of volumetric structures for each subcomponent of the DNA. A nucleotide can be represented by the combination of a base (guanine, adenine, cytosine, and thymine), deoxyribose (sugar) attached to it, and a phosphate molecule. This configuration has been used in some studies with the MCTS toolkit, Geant4-DNA (19, 20). The deoxyribose and the phosphate molecules are usually combined in one single volume, called a backbone or sugar-phosphate moiety, as done in PARTRAC (21) or other Geant4-DNA applications (22, 23). According to their van der Waals radii field, these volumes can be generated as the superposition of spherical atoms (21, 24, 25). Other shapes have also been used, such as spheres (19) or modified cylinders (14). The intention of these volumes is to accumulate energy depositions from physical processes and serve as reaction sites for chemical radicals.
TOPAS-nBio (26, 27) is a platform wrapping the capabilities of the well-validated code Geant4-DNA (28–30), offering a user-friendly interface to produce Monte Carlo simulations at the nanoscopic scale. TOPAS-nBio provides three DNA volume models, based on prominent literature descriptions, to represent the base and sugar-phosphate moieties, such as spheres (31), the ‘half-cylinder’ model proposed by Charlton and Humm (32) or the “quarter-cylinder” model proposed by Bernal and Liendo (33, 34). In the Charlton model, a pair of bases is modeled as a cylinder, surrounded by two half-cylindrical sectors representing the backbones. Each model imposes differences in DNA volume size that affect the energy scored and, in turn the assumptions for the parameters to estimate DNA damage.
Codes such as TOPAS-nBio, Geant4-DNA, or PARTRAC use MCTS simulations and volume models to perform radiobiological mechanistic modeling. These codes incorporate ways to account for different types of radiation effects. During the physical stage (in the order of 10–15 s), primary and secondary radiation tracks, physical processes, and local energy depositions are simulated according to probabilities proportional to the cross-section of each process, including excitations and ionizations of the medium, in particular, molecules of liquid water. During the pre-chemical stage (up to 10–12 s), excited and ionized water molecules give way to initial chemical species after water radiolysis. Finally, the chemical stage (from 10–12 s to 10–6 s) is typically simulated by letting these species diffuse and react with each other leading to new products. However, unlike atomistic models, volume models do not consider the molecular or atomic composition of DNA, so that ionizations and subsequent loss of chemical stability for DNA molecules cannot be explicitly simulated. Instead, 1. the probability of direct ionization is assumed to be proportional to the accumulated energy imparted into the considered volume in the physical stage, and 2. the probability of indirect ionization is assumed to be proportional to the number of reactive radicals coming into the considered volume in the chemical stage. This work explores the different parameters to be considered in volume model MCTS simulations and their relationship with the selected geometrical shapes to represent these volumes. Values of those parameters are adjusted to agree with reported experimental results in radiation chemistry.
MATERIALS AND METHODS
DNA Geometry Models
Three different volume models have been used to simulate damage of incident protons and alpha particles on DNA. Figure 1 illustrates the three models available in TOPAS-nBio (26). A pair of nucleotides (also called base pair or bp) is linked to another following the DNA double helix structure, with a thickness of 0.34 nm and a rotation of 36° per nucleotide in all cases. A hydration shell of thickness 0.16 nm was added as an external layer for the backbone on the outside of the double helix, based on the specific geometry of each model (see Fig. 1).
FIG. 1
DNA volume models used in this work: half-cylinder (left), cylinder sector (middle) and spheres (right) for backbones. Half cylinders are used to represent the bases in the cylindrical models while spheroids are used to represent the bases in the spherical model. Bases in the cylinder sector model are the same as in the half-cylinder model. Top panels show a section of 21 nucleosides or base pairs, where red and green volumes represent backbones and pink and sky-blue volumes represent bases. A 0.16-nm thick hydration shell was added to match the corresponding geometrical shape. Bottom panels show the coiling of double helix DNA around the histone volume (gray volume), represented here as a cylinder of 5.7 nm length and 6.6 nm diameter.
In all cases, the double helix had a diameter of 2.3 nm, excluding the hydration shell. Both cylinder sector and half-cylinder models shared the same bases represented as half-cylinders with a radius of 0.5 nm. Backbones in both cases were cylinder sectors with an outer radius of 1.15 nm and inner radius of 0.5 nm; this ensured that no void space was left between the backbone and base volumes. Backbone volumes in the cylinder sector and half-cylinder models were swept by 105° and 180°, respectively. The sphere-spheroid model, or simply sphere model, was specifically tailored for this work and differed from the models provided in TOPAS-nBio v1.0 (35), which consisted of spheres with diameters of 0.28 nm and 0.20 nm for backbone and base, respectively, without any representation of the hydration shell. Equal volumes were intentionally assigned for both base and backbone volumes, but a trade-off was required between maximizing the volume for each structure and avoiding overlap between one nucleotide and the next one. Taking this into account, backbones were modeled as spheres with a radius of 0.271 nm and bases as oblate spheroids with semi-major axes of 0.328 nm and semi-minor axis of 0.185 nm. The hydration shell as a spherical cap of 0.16 nm thickness covering the outer part of the backbone volume was also introduced for this work for all models.
Figure 1 bottom panel shows a nucleosome using the corresponding models: 144 bp coiled around a cylindrical histone. We built a cell nucleus for our simulations in the same way as done by Zhu et al. with TOPAS-nBio (36). Briefly, nucleosomes were arranged following a ‘beads-on-a-string’ pattern to generate a chromatin fiber with six nucleosomes per turn and 120 nm length (i.e., 15.15 kbp per fiber). Fibers were then arranged in voxels following a 3D Hilbert space-filling curve. Finally, voxels were repeated in a rectangular grid so that 6.018 Gbps were contained in a spherical nucleus with a diameter of 4.65 µm. Further details and figures illustrating the geometry can be found in Zhu et al. (36).
Representing DNA Damage Types in Monte Carlo Volume Models
Damage induced to the DNA by ionizing radiation is traditionally classified into two categories depending on the mechanisms involved: direct and indirect damage. Direct damage occurs whenever the DNA is ionized by physical processes, i.e., energy deposition in the DNA molecule itself. On the other hand, indirect damage is mediated by products of water radiolysis, such as •OH, H• and other radicals or solvated electrons (e•–aq). Becker and Sevilla (37) proposed a variation of this classification, also distinguishing the so-called quasi-direct effect. DNA molecules are typically embedded in an aqueous medium, and the water molecules directly adjacent to the DNA form the so-called hydration shell or solvation shell, which can have special properties. For instance, when molecules in the hydration shell are directly ionized, a charge from the generated electron-hole pair can be transmitted to the base and/or sugar-phosphate moieties. Once these radical cations are transferred, the resulting damage and subsequent chemical reactions are indistinguishable from those happening from direct ionization, so that from an experimental point of view, the quasi-direct effect can be considered a subtype of the direct effect. From a radiation chemistry point of view, however, these charge transfers occur with different efficiency to different sites of the DNA molecule. Thus the quasi-direct effect can be classified as a different category of ionizing radiation effects on DNA (38).
Direct effects
Direct damage to either base or sugar-phosphate moieties happens after an ionization occurs in these subcomponents, which, in turn, leads to a chain reaction of radical production, consequently resulting in the breakdown of the chemical stability of the double helix. However, as volume models do not simulate reactions to this level of detail, direct damage is usually determined by the accumulation of energy imparted into the considered volumes. A frequent approach is to consider a threshold or minimum value for the energy deposited to produce a DNA strand break (SB). In this sense, Charlton and Humm (32) obtained the value 17.5 eV as the best fit to reproduce the experiment from Martin and Haseltine in which DNA was irradiated by incorporated iodine-125 (39). A different approach, adopted by PARTRAC, for example, consists of considering the probability of strand break induction increasing with the energy deposited in the backbone volume (40). This approach is based on the experimental finding that shows that strand breaks can be induced by photons and electrons with low kinetic energies (41). In the simulations performed for this work, we have investigated both approaches: a single threshold for energy deposited into the backbone volume of 17.5 eV; and a linearly increasing probability of strand break induction from 0 to 1 for deposited energies from 5–37.5 eV. These values have been previously used in other works (40, 42). These criteria, within the base volumes, are also used to consider base damages (base damages). Note that energy is accumulated in these volumes from the interactions of a primary particle and all its secondary electrons, but not from different primary particles.
Quasi-direct effects
Ionizing radiation can also ionize strongly bound water molecules around the DNA molecule and transfer the resulting hole to the base and/or sugar-phosphate moieties. If an atomistic model of the DNA molecule is employed, then water molecules can be assumed to fill the void space between nucleosides, as performed by Aydogan et al. (43). However, for volume models, such a detailed representation is not required. In fact, hydration shells can be modeled as a water layer around the DNA double helix, that is, around the backbone volumes. It has been shown that the quasi-direct effect, i.e., damage indistinguishable from the direct effect, happens when DNA is hydrated up to approximately 11–13 water molecules per nucleotide (44, 45), which has been represented as an extra layer of 0.16 nm thickness around the backbone volume (21, 36). In atomistic models, the actual number of water molecules per nucleotide pair can be specified instead of a representative volume. For instance, Meylan et al. (18) and Tang et al. (46) include 24 molecules of water per nucleotide pair. Another approach, also used in PARTRAC (21), consists of doubling the resulting van der Waals radius for the backbone volume to include quasi-direct effect as a part of the direct effect. However, charge transfer from ionizations in the hydration shell is not uniform among DNA subcomponents: it has been reported using paramagnetic resonance spectroscopy that approximately one-third of the holes formed are transferred to the backbone and two-thirds to the base moiety (47). Although energy transfers have been considered in atomistic DNA damage simulations (48), it is not usual to include them in simulations with volume models. In this work, accordingly, each ionization of water (using the G4Ionisation process of Geant4) occurring inside the hydration shell volume has probabilities equal to 1/3 of becoming a strand break and to 2/3 of becoming a base damage. All radical productions in the hydration shells were suppressed to reflect the fact that water molecules were no longer being ionized.
Indirect effects
Determination of indirect damages requires the simulation of chemical reactions, such as radiolysis of water molecules, as well as transport of the resulting chemical products. Excited and ionized water molecules can decay into new molecules, such as H2, H•,•OH, or H3O+, whereas sub-excited electrons are thermalized, producing e–aq. To ensure the validity of an MCTS, G values, i.e., the number of species produced or lost by 100 eV of energy deposited, need to be determined (49–51). Of all these chemical species, only radicals produce damage to DNA (38), and in particular, hydroxyl radicals (•OH) are the principal inductors of strand breaks. •OH can be generated not only from radiation but also from other cytotoxic agents (52). Other species such as e–aq and peroxyl radicals can induce base damages but likely do not cause strand breaks (52). Hydroxyl radicals can attack the sugar moiety by oxidizing one of the deoxyribose hydrogens, and the probability for each hydrogen in the sugar and the base moieties to be oxidized varies (53). Since Monte Carlo volume models do not resolve hydrogen atoms, probabilities are assigned to induce either a strand break or base damage for each •OH reacting with DNA. In the absence of direct measurement data, a reasonable approach is to adjust these probabilities based on radiochemical results, either experimental or modeled. In this sense, we used: (a) that the proportion of hydroxyl radicals attacking bases and backbones distributed as 80% and 20%, respectively, as shown in multiple experimental and modeling studies (37, 54); and (b) that about 11–13% of the •OH reacting with DNA molecules lead to strand breaks, as modeled from measurements in plasmids dissolved in unscavenging solutions by Milligan et al. (55). Based on these two results, Nikjoo et al. (13) proposed that a strand break is generated in linear DNA with a 13% probability whenever a hydroxyl radical interacts with a nucleotide as a whole, i.e., the base or backbone; whereas Friedland et al. (40) used a probability of 65% in PARTRAC to generate a strand break for hydroxyl radicals interacting only with the backbone volume, which implicitly assumed a 20% of probability of interaction with the backbone to give 13% of interactions leading to strand break. Later studies using Geant4-DNA adjusted this latter probability to 40% according to their geometry. (19, 22). Zhu et al. (56) quantified the impact of this in TOPAS-nBio by about 40% difference in strand breaks and 50–70% difference in double-strand breaks (DSBs). We propose that to meet the results (a) and (b), these values are, in fact dependent on the geometry and scoring approach used and should be therefore adjusted. Assuming that the number of hydroxyl radicals diffused within each DNA subcomponent is proportional to their volume, then the probability of a reaction with a backbone should be: firstly, four times lower than reacting with a base to meet (a); and secondly, proportional to the ratio of backbone and base volumes. Table 1 shows the volumes of the three DNA geometry models and the corresponding probability for a hydroxyl radical interacting with a backbone Pint,backbone.
TABLE 1
Volumes of the Subcomponents of a Nucleotide with the Three Models Shown in FIG. 1 and the Calculated Probability of a Backbone and •OH Reacting with Each Other
The probability for a hydroxyl radical to react with a base was set to 1 in all cases, so that the corresponding Pint,backbone for each model results in 20% of the reactions happen in backbones and 80% in bases. To meet (b), once a hydroxyl radical reacts with a backbone it is scavenged, and a probability of 55% to become a strand break was assigned, resulting in 11% of the total number of reactions with the DNA structures causing strand breaks. All interactions with bases were scored as base damage and their tracks were stopped and terminated. Furthermore, e–aq were also allowed to react with bases to cause base damages, but not with backbones to induce strand breaks. Such explicit considerations for strand breaks and base damages induction are not common in Monte Carlo volume model simulations, however, the validity of our assumptions can be compared to another reported result: the yield of base damages should be on the order of 2.5–3 times the yield of strand breaks in the absence of consideration of peroxyl radicals, which can attack bases but are unlikely to produce strand breaks (57, 58). Figure 2 illustrates our approach based on two independent steps to determine indirect damage. The scoring approach adopted in previous studies with TOPAS-nBio (36) and many others, for example Henthorn et al. (34) is also illustrated.
FIG. 2
Left side: Method to score damage from chemical tracks used by Henthorn et al. and previous studies using TOPAS-nBio. Right side: Method used in this work. Yellow stars represent hits within the backbone volume, i.e., processes susceptible to attack the backbone; while the blue star represents a hit within the base volume. On the right picture, not all hits in the backbone react with the backbone so the track continues to be transported.
Simulations Performed in this Work
We simulated monoenergetic particles traversing a cell nucleus. A uniform field of particles was generated at 4.65 µm distance from the nucleus center. For each new event, i.e., primary particle, the nucleus was randomly rotated to avoid dependencies of the result on the orientation of the voxelized geometry. To evaluate the dependence of the damage induced on the linear energy transfer (LET), we simulated protons and alpha particles of varying energy. Protons with kinetic energy 0.5 MeV, 1 MeV, 2 MeV, 5 MeV, 10 MeV and 20 MeV (i.e., LET ranging from 2.6 keV/µm to 41.3 keV/µm) and alpha particles with energies 1 MeV, 2 MeV, 4 MeV, 6 MeV and 8 MeV (or LET from 63 keV/µm to 219.0 keV/µm at the nucleus entrance) were simulated. In every case, tracks were sampled within parallel beams of width equal to the nucleus cross-section. Physics processes were handled by the default constructor of TOPAS-nBio v1.0 (51), based on Geant4-DNA (59), whose detailed models and cross-sections for each particle and energy range can be found elsewhere (60).
The chemical stage was simulated using a step-by-step approach with the default diffusion coefficients and reaction rates of TOPAS-nBio v1.0 (51), using TOPAS v3.6.1, which in turn is built on top of Geant4 v10.06.p01. TOPAS-nBio takes the chemical species resulting from water radiolysis after the pre-chemical stage and simulates their diffusion by sampling the direction and distance traveled (defined by a given time step) according to Brownian motion, based on the Smoluchowski equation in 3D. Diffusion coefficients are temperature-dependent and are obtained from fits to experimental data (51). Reactions between species in TOPAS-nBio and Geant4-DNA are diffusion-controlled, meaning that reaction times are negligible compared to the time for both species to diffuse nearby (61). Then, chemical species are allowed to interact with each other if they are located within a given interaction radius, and different reactions are assigned with probabilities to happen to match reaction rate constants after radiolysis of water at different temperatures. Detailed explanations of these mechanisms, values used and validation of these can be found in the chemistry-specific TOPAS-nBio publications (51, 62, 63). In this work, tracks of chemical species were simulated up to 1 ns with a time step of 0.5 ps for the species to diffuse, following other similar studies (19, 36). The entire nucleus was composed of liquid water. The histones were composed of liquid water of density 1.407 g/ cm3. Histones were considered as scavengers, i.e., all chemical species entering the histone volume had their track stopped and terminated.
Simulations were repeated for the three DNA volume models shown in Fig. 1. We studied the effects of: (a) varying the method to score direct damage using an energy threshold for the probability and a linear probability increasing with the energy deposited; and (b) varying the probability Pint,backbone or a hydroxyl radical to react with the backbone. Statistical uncertainties are indicated by error bars in the figures, showing one standard error in the mean value. Besides the total number of strand breaks we also considered their association in DSBs, defined as two strand breaks in opposing strands of the DNA separated by 10 bp or less. Yields of single-strand breaks (SSBs), i.e., the remaining strand breaks not forming a DSBs, were also counted. If multiple DSBs and SSBs were concentrated in 10 bp or less, they were tallied individually. Therefore, the total number of strand breaks can be obtained as the number of SSBs plus twice the number of DSBs. We compared our results with those from simulations reported in recent studies by Kundrát et al. (64) using PARTRAC as well as Sakata et al. (23) using Geant4-DNA.
RESULTS
Direct and Quasi-Direct Damage
Figure 3a and b shows the yield (damage per Gy per Gbp) of direct base damages and strand breaks produced by protons (green shading and left to the vertical dashed line) and alpha particles (red shading and right to the vertical dashed line) without considering any quasi-direct damage, i.e., as if there was no hydration shell around the backbone volume. For these results a single threshold of 17.5 eV imparted in the corresponding volume was utilized to determine if a base damage or a strand break took place. Figure 3c shows the average number of strand breaks for proton simulations as a function of the volume of the backbone for each model, showing that this number of direct strand breaks (without quasi-direct effect) is proportional to the volume of the backbone.
FIG. 3
Panel a: Base damages and panel b: strand breaks per Gy per Gbp produced by protons (green shading and left to the vertical dashed line) and alpha particles (red shading and right to the vertical dashed line) by direct effect without considering the quasi-direct effect for the three geometry models considered. Panel c: Shows the average number of strand breaks produced by protons with any LET as a function of the volume of the backbone structure in each model.
The impact of the threshold for damage on the direct effect i.e., either a single energy threshold of 17.5 eV or an increasing linear probability from 5 eV to 37.5 eV, is shown in Fig. 4. As observed, the linear method tends to increase the recorded damage compared to the threshold method, but its impact is not uniform across different geometries. In particular, the half-cylinder model seems relatively insensitive to the method selected, while the cylinder sector and sphere models show larger differences for direct damage (without including quasi-direct damage).
FIG. 4
Yields of base damages (BDs) (top panels) and strand breaks (SBs) (bottom panels) from direct effects without considering quasi-direct effects for the three geometry models considered, i.e., cylinder sector (left), half-cylinder (middle) and sphere (right). Results when using a single threshold of 17.5 eV for the energy accumulated in a volume or an increasing linear probability from 5 eV to 37.5 eV are compared. Green shading (or the area left to the vertical dashed line) represents the LET range for protons and red shading (or the area right to the vertical dashed line) the LET range for alpha particles.
Damages attributable to the quasi-direct effect are included in Fig. 5, in which the total direct damage is also shown, that is, the accumulation of energy in bases and backbones plus ionizations in hydration shells. The quasi-direct effect is approximately twice as much in bases as in backbones as expected (see section representing DNA damage types in Monte Carlo volume models). The energy threshold of 17.5 eV to record base damages and strand breaks was used.
FIG. 5
Total direct damage and quasi-direct damage recorded for the three DNA models in FIG. 1 using the single threshold method to score base damages (BDs) (left panel) and strand breaks (SBs) ([right panel) for protons (green shading and left to the vertical dashed line) and alpha particles (red shading and right to the vertical dashed line). The legend is applicable to both panels.
Indirect Damage
We investigated the effects of varying the probability for a hydroxyl radical to react with a backbone Pint,backbone. As pointed out, this parameter should be specific for each geometry to meet the condition of distributing the reactions with bases and backbones as 80%/20% as well as 11% of the reactions with the DNA leading to strand breaks. We varied Pint,backbone for each geometry as shown in Fig. 6, where we also compare the results among geometries using the probabilities calculated in Table 1. Using these probabilities brings the predicted yields of strand breaks into agreement across the three models and within the 2.5-3 yield range suggested by Ward (57). Due to large computational requirements, only selected points for alpha particles were simulated to establish overall trends.
FIG. 6
Yield of indirect strand breaks (SBs) (top panels) and ratio of indirect base damages (BDs) and SBs (bottom panels) for the three DNA geometry models (from left-to-right cylinder sector, half-cylinder and sphere) for protons (green shading and left to the vertical dashed line) and alpha particles (red shading and right to the vertical dashed line) for various probabilities Pint,backbone for a •OH to react with a backbone. The right panel shows both quantities when an adjusted probability, as indicated in Table 1, is utilized for each geometry.
Yields of SSBs and DSBs
Combining strand breaks for direct damage (using the threshold method set to 17.5 eV), quasi-direct damage, and indirect damage (with Pint,backbone adjusted for each geometry according to Table 1), we investigate the induction of DSBs. Figure 7 shows the breakdown of the damage sites by their origin, i.e., direct, quasi-direct, or indirect. Hybrid DSBs, i.e., DSBs produced by the combination of a direct or quasi-direct strand break and an indirect strand break, are also shown. Major differences are due to the discrepancies in direct damage when using the same threshold for different geometries, as shown in Fig. 3. On the other hand, both indirect SSBs and DSBs are also consistent along with geometry models when using the adjusted probabilities suggested in Table 1. Ratios between base damages and strand breaks for all three geometries for direct and indirect effects are shown in Fig. 7.
FIG. 7
Top rows: Ratios between base damages (BDs) and strand breaks (SBs); middle rows: yields of SSBs; bottom rows: yields of DSBs for half-cylinder model (left column), cylinder sector model (middle column) and sphere model (right column). For BD/SB ratio and SSB yields: cyan points indicate quasi-direct damage; blue points indicate direct damage (including quasi-direct) with single threshold method; red points indicate indirect damage (with Pint,backbone adjusted according to Table 1); and black points indicate total damage (aggregated of direct and indirect). For DSB yields: blue points indicate fully direct DSBs (combination of direct/quasi-direct SB and direct/quasi-direct SB), red points indicate fully indirect DSBs, green points hybrid DSBs (combination of direct/quasi-direct SB and indirect SB) and black points indicate total yield of DSBs (aggregated of direct, indirect and hybrid).
Figure 8 shows our results using the cylinder sector DNA model, a single threshold of 17.5 eV for the energy deposited in a volume to score direct damage; and the probability of reactions with a backbone equal to 25%, compared to recent simulations of protons and alpha particles using PARTRAC by Kundrát et al. (64) and protons using the Geant4-DNA model developed by Sakata et al. (23). The latter only shows a breakdown of the strand breaks depending on their origin (direct or indirect) but not for either SSBs or DSBs. Neither of the two articles considers damages on base damages.
FIG. 8
Yields of strand breaks (SBs) (left panel), SSBs middle panel and DSBs (right panel) obtained in this work using our cylinder sector DNA model with single threshold of 17.5 eV for direct damage induction and Pint,backbone = 0.25; compared to results published by Kundrát et al. (57) using PARTRAC and Sakata et al. (20) using Geant4-DNA. Green shading (or the area left to the vertical dashed line) represents the LET range for protons and red shading (or the area right to the vertical dashed line) the LET range for alpha particles for all studies shown.
DISCUSSION
Different volume models can be employed to represent the DNA subcomponents, and the choice can be motivated by several reasons, such as historical, further compaction to represent larger structures, limitations of the Monte Carlo codes, or simplification for the sake of computational speed. However, one should expect consistency among the results produced by various Monte Carlo codes independently of the volume model selected, as each code aims to reproduce experimental results as much as possible. Yet, biology experiments also suffer from multiple factors leading to uncertainty, such as the quantification of strand breaks separately from the repair, the biological environment of the DNA, the presence of chemical scavengers, and the hydration level. Defining the ground truth that can be compared to results from Monte Carlo simulations is not easy. Nevertheless, experimental key values, e.g., the ratio of damages in bases and backbones, should be reproduced by all simulations.
In this study, we show that parameters to convert both direct and indirect damage are geometry dependent. For example, the threshold value of 17.5 eV was obtained by Charlton and Humm (32) to adjust results using their particular geometry and scoring method, and as we showed in Fig. 3, scored direct damage depends on the volume model employed. It seems reasonable, therefore, that the threshold value should be adapted for each Monte Carlo volume model. No rule for this adaptation seems trivial, as an experimental reference needs to be firstly selected, and then different simulations varying the energy values may be needed to obtain the best fit for the experiments. In fact, Thompson et al. (65) have recently revisited the analysis done by Martin and Haseltine to assess the convenience of this energy threshold using TOPAS-nBio, showing that higher values might be necessary to reproduce experimental data. In this work, we show the disparity of results when the same configuration (either the same energy threshold or the same linear probability energy range) is employed but with different volume models. In addition, it remains unclear whether using the same criterion to consider strand break and base damage for the energy deposited in backbones and bases, respectively, is justified. If the efficiency at producing damage can be considered similar for all DNA molecules, then a sub-component-wise energy threshold should also be adapted to account for volume differences. Regarding quasi-direct effects, the hydration shell in our geometries represents just an approximate model since the underlying mechanisms, i.e., charge transfers, cannot be simulated with volume-model-based MCTS simulations. Therefore, as happens with the base and backbone models, once a geometry model for the hydration shell is selected, additional parameters could be necessary to adjust the efficiency of the induction of base damages and strand breaks by charge transfers from the hydration shell to match experimental results.
To obtain the correct probability to include indirect stand breaks, Henthorn et al. (34) proposed the adjustment of the reaction probability to match the proportion of 35–65% between direct and indirect damage when simulating low-LET irradiations, which has been broadly reported (57). However, this adjustment depends on the criteria to include direct damage. Furthermore, reports of the proportion of 35-65% for direct-indirect damage mainly refer to cell killing and not necessarily strand break induction (66). By contrast, the criterion of assuming 11–13% of the •OH reacting with the DNA leading to strand breaks should be insensitive to irradiations with different LET and independent of how direct damage is scored. Nonetheless, this result was measured by Milligan et al. (55) in supercoiled plasmids transitioning to open circular plasmids through damage, and its repeatability in the cellular medium is not guaranteed. Indeed, these results are shown to considerably vary when plasmids are dissolved in scavenging solutions (55). Figure 9 shows the proportion of direct damage with respect to the total damage scored with each of our models in terms of strand breaks and DSBs. The cylinder sector model predicts the closest result to a 35% proportion of direct damage for both strand breaks and DSBs. Nonetheless, this is a manifestation of the same phenomenon described in the previous paragraph: while the indirect damage is consistent among geometries using a dedicated probability Pint,backbone, the direct damage differs because we used the same energy threshold to determine when strand breaks occur.
FIG. 9
Proportion of direct to total damage for strand breaks (SBs) (red) and DSBs (black) for the three DNA volume models considered in this work, using the same energy threshold for direct damage and adjusted probability for indirect damage.
In previous studies (31, 36), hydroxyl radicals were assumed to react with the backbone whenever they entered the volume, and therefore, their tracks were stopped and terminated. Each one of these events had a given probability assigned to lead to a strand break. This approach likely underestimates the number of base damages, particularly for the cylindrical geometries in which the bases are partially or fully shielded by the backbones since the vast majority of •OH are killed before they can reach the base volumes. With our approach, however, •OH tracks can reach the base volumes with a given probability, which, in turn, is adjusted for each geometry to keep consistency between the relative number of reactions to base and backbone moieties, as shown in Fig. 2. Also, note that we assumed a perfect efficiency for radicals to induce base damages, but this efficiency could be freed as a new parameter to control the base damage/strand break ratio, enlightening the possibility of using different geometry DNA models while conserving the same results. This is of interest since different geometry models may be needed to build increasingly complex geometries such as chromatin, chromosomes, and cell nuclei.
Monte Carlo simulations of damage to DNA typically do not report base damages. In this work, we consider not only the yields of base damages but also include a new geometry (sphere model) for which the overall ratio between base damages and strand breaks is shown to be within the 2.5–3 expected range reported in the literature (57). The relevance of a correct determination of base damages is not clear since their repair seems to be simpler for a cell than complex strand breaks (57). However, although DSBs have been proven to correlate with cell killing, the complexity of damage produced in a given DNA segment or site, including base damages, seems to be inversely correlated with the likelihood for this damage to be repaired (67). Despite the simpler repair, base damages may increase the complexity of DNA damage. For example, base excision repair (BER) is the main mechanism involved in the base damage repair, and although it is usually a slow process compared to strand break repair pathways, it has been suggested that parallel recruitment of components for different repair pathways may impair the repair efficiency (68). TOPAS-nBio allows for a standard record of not only strand breaks but also base damages using the standard DNA damage (SDD) format proposed by Schuemann et al. (69), which can be used as an input for systems modeling the DNA repair (70, 71). The methods and ideas presented in this work have been implemented in TOPAS-nBio. Therefore, TOPAS-nBio can be used to discriminate if strand breaks and base damages are in good agreement with experiments and radiation chemistry models (57).
Our results with the cylinder sector model, which for the selected energy threshold reproduces the 35–65% proportion of direct-indirect damage at low LET, approach those from PARTRAC (64) for both protons and alpha particles. Differences are attributable to using the single energy threshold for accounting for direct damage instead of a linearly increasing probability (see Fig. 4). The indirect number of strand breaks (and SSBs) in the Geant4 study from Sakata et al. (23) seems considerably higher than those calculated in this study and by PARTRAC. Potential differences between the study of Sakata et al. and our approach are that they used the geometry from Lampe et al. (19), in which sugar and phosphate are modeled as independent structures instead of a combined single backbone volume. They also employed an independent reaction time (IRT) method to propagate chemical tracks instead of the step-by-step method employed both in this work and in PARTRAC. Although Kundrát et al. used the same definition as in this work for the consideration of DSBs, Sakata et al. separated damages as occurring when a sequence of unbroken DNA greater than 100 bp was found. Then, each damage site was classified as SSB or DSB depending on whether a double-break in complementary strands was found within ten bp. Therefore, multiple DSBs concentrated separated by less than 100 bp were considered a single DSB. Nevertheless, the number of total DSBs in their study is in good agreement with the ones shown in this work and by PARTRAC. These differences illustrate how the scarcity of consistent experimental data with the same cell environment, including hydration level and scavenger concentration, hinders determining a ground-truth method for Monte Carlo simulations.
CONCLUSIONS
Volume models for Monte Carlo simulations of DNA damage can encompass different shapes and sizes: cylinder sector, half-cylinder, and sphere-spheroid. Recording the damage from direct, quasi-direct, and indirect damage in Monte Carlo simulations, therefore, requires the selection of characteristic parameters. Multiple values for Monte Carlo parameters important to determine DNA damage induction can be found in literature and are sometimes propagated between studies. However, this study shows that those parameters should be specific to the employed models, and they must be selected to meet experimental results. In particular, the threshold of energy imparted to a volume for direct damage consideration needs to be specific to the employed volume. The probability for hydroxyl radicals to induce indirect damage should further be chosen to meet basic results from radiation chemistry: backbones receive 20% of the hydroxyl induced damages to the DNA, and 11–13% of these damages induce strand breaks. We also present results of the yield of base damages as an additional checkpoint for validation of our methods. This yield is expected to be between 2.5 and 3 times the yield of strand breaks, which was achieved using the sphere-spheroid model. TOPAS-nBio integrates different volume models at the nucleotide scale. Here we present methods and parameters to produce consistent results among the models provided by TOPAS-nBio.
ACKNOWLEDGMENTS
This work was supported in part by the National Institutes of Health/ National Cancer Institute (NIH/NCI) grant no. R01 CA187003: “TOPAS-nBio: A Monte Carlo tool for radiation biology research” and the NIH/ NCI grant no. K99 CA267560: “Radiation dosimetry for alpha-particle radiopharmaceutical therapy and application to pediatric neuroblastoma.”
