Chromosome aberrations (CAs) are one of the effects of radiation exposure and can have implications for human health in the space environment, since they are related to cancer risk. In radiation research, chromosome aberrations are a convenient biomarker for carcinogenesis. To shed light on the formation and quality of chromosome aberrations in the space environment, many experiments and simulations have been performed using chromosome aberrations in human cells, induced by heavy ions, which are present in galactic cosmic rays (GCRs). In this work, the new simulation program, radiation-induced tracks, chromosome aberrations, repair and damage (RITCARD), is presented. This software program is based on the algorithm used in the NASA Radiation Track Image (NASARTI) model with some improvements. NASARTI and RITCARD are both comprised of four parts: a random walk (RW) algorithm for simulating chromosomes in a nucleus; a deoxyribonucleic acid (DNA) damage algorithm; a break repair process; and a function to assess and count chromosome aberrations. Prior to running RITCARD, the code, relativistic ion tracks (RITRACKS), is used to simulate detailed radiation track structure and calculate time-dependent differential voxel dose maps in a parallelepiped centered on a cell nucleus. The RITCARD program reads the pre-calculated voxel dose and locates the intersections between the voxels and the chromosomes that were simulated by random walk. Radiation-induced breaks occur strictly at these intersections with a probability that is a function of the voxel dose. When a break occurs in the random walk, the corresponding chromosome piece is cut into two fragments where each has a free end at the position of the break. RITCARD generates a collection of all fragments, free ends, and enlists free end pairs. In the next step, the algorithm simulates the time-dependent rejoining of free end pairs, using different probabilities for pairs originating from a given break (proper) or from different breaks (improper), which results in the formation of fragment sequences. By grouping these sequences, the program determines the number and types of aberrations, based on the same criteria used in our experiment. The new program is used to assess the yields of various types of chromosome aberrations in human fibroblast cells for several ions (1H+, 4He2+, 12C6+, 16O8+, 20Ne10+, 28Si14+, 48Ti22+ and 56Fe26+) with energies varying from 10 to 1,000 MeV/n. The results show linear and linear-quadratic dose dependence for most chromosome aberrations types. The calculation results were compared with those obtained by fluorescence in situ hybridization (FISH) experiments that were performed by our group. The simulations and experiments are in better agreement at lower LET. Regarding the simulation results, the coefficient of the linear part of the dose-dependence curve also peaks at an LET value of approximately 100 keV/lm, which evokes a relative biological effectiveness (RBE) peak found by other researchers.
Radiation poses a hazard for human spaceflight and is of particular concern for future manned missions to Mars, during which astronauts may be exposed to doses of approximately 1 Sv (1). Radiation in the space environment comes from trapped particle belts, solar particle events and cosmic rays (2). Protons (87%), helium ions (12%) and heavier ions (1%) compose the galactic cosmic ray (GCR) spectra. The health risks from exposure to these radiation types are difficult to quantify, notably because of limited knowledge of the biological effects of high-energy and high-charge (HZE) particles. In particular, HZE particles (fully or partially ionized ions from deep space) typically have a higher relative biological effectiveness (RBE) than low-linear energy transfer (LET) X rays and gamma rays for many biological end points, notably deoxyribonucleic acid (DNA) and chromosome damage (3).
Chromosomal aberrations (CAs) have been used in mechanistic studies of radiation effects, in biodosimetry and as possible biomarkers for carcinogenesis (4). The frequency of CAs in circulating lymphocytes is considered by many investigators to be a reliable indicator of biological dose. Indeed, it has been used to measure radiation exposures in nuclear accidents (5), in survivors of the Hiroshima and Nagasaki atomic bombings (6) and in astronauts exposed during spaceflight (7). Therefore, many experimental and theoretical studies have been conducted on radiation-induced CAs over the past few years, including mechanistic studies, biodosimetry, genomic instability and investigations of possible biomarkers for carcinogenesis (8–14), and more experiments are still ongoing. One of the difficulties in assessing chromosome damage experimentally is the need for an available technique that can provide a high-resolution assessment of all types of chromosome rearrangements simultaneously (15).
In this work, a new simulation program was developed and is presented here; the radiation-induced tracks, chromosome aberrations, repair and damage (RITCARD) program is used to simulate radiation-induced CAs in cells. This code is inspired by the program, NASA radiation track image (NASARTI) (8, 9). RITCARD was developed using modern software engineering practices to enable current and future model advancements (i.e., mixed-field and dose-rate simulations) to be implemented more efficiently and reliably. Although RITCARD and NASARTI are based on the same fundamental theoretical framework, the numerical implementation is different, with the new code, RITCARD, being more reliable, portable, efficient and well documented. NASARTI and RITCARD are both composed of four parts: a random walk (RW) algorithm for simulating chromosomes in a nucleus; radiation-induced DNA damage calculations; a break repair process; and a function to assess and count the types of CAs generated. RITCARD was used to calculate the yield of various types of CAs occurring in human fibroblast cells for many ions (1H+, 4He2+, 12C6+, 16O8+, 20Ne10+, 28Si14+, 48Ti22+ and 56Fe26+) with energies varying from 10 to 1,000 MeV/n. The simulation results show linear and linear-quadratic dose dependence for most CA types, and the yields are relatively comparable to available experimental data. The coefficient of the linear part of the dose-dependence curve also peaks at a LET value of approximately 100 keV/µm, which evokes a RBE peak. Furthermore, the quadratic coefficient of the dose-response curve decreases with LET, suggesting that CAs caused by delta rays originating from several radiation tracks are more likely to occur at low LET. This program will be used to help analyze and understand results from CA experiments such as those performed recently at the NASA Space Radiation Laboratory, Brookhaven National Laboratories (NSRL).
MATERIALS AND METHODS
Code Programming and Maintenance
RITCARD is composed of several programs. The codes were made using the C++ programming language, version 98, using the Borland C++ environment. C++ classes are used extensively to model the various structures used in the program, such as chromosome fragments and reactive free ends. The classes are useful for modeling the objects, their attributes and their behaviors. The calculation part of the code is a console application that is platform-independent and also compiles and executes in Linux. It is also designed for execution in parallel on many central processing units (CPUs), both on a Linux cluster and on a Microsoft® Windows® (Redmond, WA) platform, as described elsewhere (16). In the Windows environment, a graphics user interface (GUI) is used to set up simulation parameters, whereas this is done manually in the Linux environment. Another GUI can be used to display the results in Windows. The tracks and the CAs can also be displayed with the RITRACKS OpenGL interface. That is, the Windows version is used mostly for development and testing of the code, whereas the Linux version is used for parallel execution of the code with large number of histories (production runs). The code is maintained and changes are tracked using the in-house GitLab environment accessible to NASA users, in compliance with NASA software engineering standards. A method toxicity metrics analysis was performed to assess the number of statements, the number of parameters, the depth of “if” functions and the cyclomatic complexity for each function of the code. For the calculation part, the maximum toxicity metrics score is 1.117, slightly exceeding the recommended threshold of 1.0. This score indicates that the code is well structured and in compliance with coding standards.
Part 1A: Simulation of chromosomes by random walk. The development of the chromosome simulation program that is a part of NASARTI spans several years (17, 18). Several updates have been made during this period, the most notable one in 2006 that included chromatin loops (19). Beginning in 2012, a CA simulation algorithm was implemented in NASARTI, but some coding errors were discovered during the course of this project and were resolved only recently. The RW model is based on the following principles:
1. The chromatin has physical and scaling properties analogous to those of long polymer chains;
2. A chromosome is roughly described by a random coil physical model;
3. Each chromosome is predominantly located in a part of the cell nucleus (domain);
4. Many large and intermediate-size DNA loops are present;
5. Telomeres are located at the extremities of the chromosomes, i.e., at the first and last ends of the RW.
Therefore, the simulation of chromosomes is done by a discrete RW, meaning that each chromosome is approximated by a coil comprising a series of monomers of 0.02 µm, each monomer containing 1 kbp of DNA. The nucleus can be modeled either as a sphere (e.g., for lymphocytes) or as an ellipsoid (e.g., for fibroblasts). In this work, since the simulations are done with fibroblast nuclei, the ellipsoid shape (height: 6 µm; large axis: 17 µm) is used. The extremities of the monomers occupy the points of a cubic lattice (X,Y,Z), where X, Y and Z are integers. The beginning of the chromosome (the first telomere) is determined by generating a position (X,Y,Z) randomly in the corresponding domain. The position of the next monomer is determined randomly and can be left (X-1,Y,Z), right (X+1,Y,Z), backwards (X,Y-1,Z), forward (X,Y+1,Z), down (X,Y,Z-1) and up (X,Y,Z+1). The process is repeated until the last monomer (the last telomere) is simulated.
The probability for each direction is equal at 1/6, but several constraints exist so that some directions may be forbidden during one step; therefore, these probabilities are different for some monomers. The new position of the monomer must be located inside the volume of the nucleus and in the domain of the chromosome. The option to use intermediate-size loops is also available. As described elsewhere (19), these loops are 60 monomeric units in length, separated by fragments of 60 monomer units. The program generates a subset of RW of 60 steps, until the last step of the RW returns to the original point of the loop.
The chromosome RW model included in the program NASARTI has been extensively validated, notably by pulse-field gel electrophoresis experiments (17). However, these calculations were made using large doses and amorphous track structures. For this work, since many new features have been added and stochastic track structures are used instead of amorphous tracks, the RW simulation routine was isolated from NASARTI, and is now a stand-alone program. Therefore, the RW simulation program pre-calculates nuclei that are used in RITCARD.
Part 1B: Simulation of volume irradiation by stochastic radiation tracks. Stochastic radiation tracks were simulated with the program RITRACKS, which has been described elsewhere (20). Briefly, RITRACKS calculates the energy deposition events, ionizations and excitations of water molecules by protons and HZE particles, as well as the energies, positions and directions of the secondary electrons and their full transport in an irradiated volume. The modeled cell nucleus was centered in a parallelepiped, which is considered the irradiated volume. Using calculated energy deposition, dose was calculated in voxels of 20 nm (21).
Previous versions of RITRACKS were used to simulate the irradiation of a parallelepiped or cylindrical volume with multiple tracks. However, the number of tracks was required to be fixed and identical for every simulation history. The type of ion and the energy were also required to be identical for each simulated ion. Additionally, all tracks are generated at time, t = 0. Recently, a new mixed-field functionality was added in RITRACKS. The new mixed-field functionality allows the simulation of the irradiation of a volume by multiple ion types and/or energies. The irradiation time should also be specified for each ion type and energy, which allows the program to simulate consecutive and/or overlapping beams.
For each simulation history, the number of tracks of each type impinging the volume is obtained by sampling from the Poisson distribution,
where λ = ϕA is the average number of tracks, ϕ is the fluence and A is the area of the irradiated surface. The value of ϕ is obtained from a given dose using the well-known equation:
The LET of a particle is calculated using the Bethe equation (22).
The user also needs to specify the time interval for each particle type. It is therefore possible to simulate consecutive or overlapping beams using the code. The GUI to enter mixed fields data is shown in Fig. 1.
The grid shows the summary of each ion type and energy. The contribution of each beam type can be modified. In principle, there are no limits on the number of particle types and energies to use in the simulation. For this work, the calculations are limited to one ion type of one energy, with acute dose rate (i.e., all tracks are simulated at time, t = 0).
Calculations performed using RITRACKS led to the discovery that the dose calculated in the nucleus was significantly lower than the requested dose, by up to 25%, especially for ions with high energy per nucleon. An investigation of this problem revealed that a significant amount of the energy is not deposited in the irradiated volume due to secondary electrons escaping. To improve the calculated dose, we first attempted to increase the size of the irradiated area. However, for particles of moderate to high energies, many secondary electrons have sufficient energy to leave the irradiated volume. For instance, electrons with energy in the MeV range can travel up to 100 mm (23); thus, increasing the size of the irradiated area did not help much. In actual experiments, we would also expect secondary electrons from neighboring volumes to contribute to energy deposition in the volume under consideration. To consider this, we simulated the contributions from the electrons originating from outside the irradiated volume by using periodic boundary conditions (PBC). We have found that using PBC significantly improves the dose calculations, with values very close to the predicted values given by using the Bethe equation. However, PBC is an imperfect approach. For example, at low doses the electron tracks that return to the irradiated volume after crossing a boundary represent secondary electron tracks generated by ion tracks impinging neighboring volumes. At low doses these tracks occur at a time that is different from the ion track actually impinging the irradiated volume. For this work, since acute dose rates were used for the simulation, this is not expected to influence the results.
Part 1C: Simulation of chromosome breaks. Previously, NASARTI was used to calculate the number and distribution of double-strand breaks (DSBs) within nuclei, using amorphous track calculations of the radial dose combined with chromosome models simulated by random walk (8, 9). In this code, which is used essentially for comparison in this work, the probability (Ψ) to find at least one DSB at a monomer located in a voxel is given by Eq. (3),
where D(i,j,k) is the local (voxel) dose, i, j and k are the spatial coordinates of the voxel (in lattice units), and Q is an adjustable parameter representative of the intensity of the stochastic process of DSB formation. In this work, Q = 1.14 × 10–5 Gy–1. Since several tracks (ntracks) may be present, D(i,j,k) is calculated as
where rl is the radial distance to track l. In RITCARD, the simulation of chromosome breaks is done the same way, but stochastic tracks are used instead of amorphous tracks. Equation (3) is also used for calculating the breaks at specific monomers of the chromosomes. However, since stochastic tracks are used, the voxel dose values are obtained from the three-dimensional (3D) voxel dose map (21). To include possible dose-rate effects, the creation time of the break is also obtained by using the corresponding time in the voxel. A schema summarizing the first part of the program is shown in Fig. 2.
Part 2: Repair of breaks. In Part 1, the program calculates the break locations in the chromosomes. Similar to a rope, each break results in the formation of two fragments. An example of fragments generated by the program is shown in Fig. 3.
The newly created extremities are reactive free ends. The free end pairs are organized into a list, and each reactive free end pair can be repaired. The repair of reactive free ends that originates from the same break (e.g., the free ends no. 2 and 3 in Fig. 3) is a proper rejoining; otherwise, it is an improper repair. The relevant files generated by the program are described in detail in Appendix 1.
In this part of the program, the breaks are repaired using an algorithm that follows the structure of NASARTI. The program uses time steps of 5.18 s, for 24 h. The value of this time step was obtained from prior calibration of the model, in which the DSB repair kinetics are fitted as a bimodal exponential curve (9). At each time step, the program picks a free end pair randomly. If the free ends originate from the same break, it is a proper rejoining, which is repaired with probability P = 1. Otherwise, it is an improper repair. In this case, the probability of repair depends on the 3D Euclidian distance between broken ends as in Eq. (5),
where I, is the probability of an improper rejoining (misrejoining) in an elementary step of the algorithm; W is an empirically calibrated parameter; r is a 3D Euclidian distance between reacting DSB free end (µm); and σ is an adjustable parameter with the dimension of microns. Another condition that needs to be fulfilled for a break to repair is that the break must have been already created at the time of repair. This condition is added to consider possible dose-rate effects. When a free end pair is repaired, all free end pairs containing these free ends are removed from the list of free end pairs, so they do not participate in the repair process any further.
After the repair algorithm is complete, RITCARD analyzes the repaired chromosomes. At this point, RITCARD diverges from NASARTI; that is, it is using a different, but equivalent, algorithm. Briefly, the program is building new objects, named fragment sequences, which are composed by chromosome fragments. A new sequence begins with the first fragment comprising a free end. Then, if a free end is present on the right side, the program inspects the list of repaired free ends to find if it has been repaired. If it is the case, the fragment to which the other free end of the pair belongs is added to the right of the sequence. If the last free end of the fragment is selected, the added fragment is in reverse order. The other end of the fragment is then examined. If it is a telomere, the repair process is over from this side; otherwise, the process is repeated until the last free end is a telomere or not repaired. The repair process is complete on the right side. If the initial fragment does not begin with a telomere, the repair process is repeated from the left side. The repair algorithm is illustrated in Fig. 4.
Part 3: Classification of chromosome aberrations. During the repair process, the fragments are rejoined, which leads to the formation of fragment sequences. Figure 5 illustrates one possibility of fragment sequences for the fragments shown in Fig. 3 (a more detailed file is shown in the Appendix).
The fragment sequences generated are not only a clear visualization of the formation of CAs but also highly useful for verifying that they are correctly identified by the algorithm.
A detailed review of CA types is given in (24). The code detects intact chromosomes (properly repaired), terminal deletions, deletions, inversions, intrachromosomal translocations, simple exchanges, rings, dicentrics, open fragments and complex aberrations. These aberration classes are not necessarily exclusive; for example, a ring can also be a dicentric.
Group Sequences with Material from the Same Chromosomes (Fig. 5)
The classification algorithms (Figs. 6–8) consider each group of fragment sequences separately. In the beginning of this part, the program determines the number of chromosomes present in the group and analyzes the sequences accordingly. The groups containing fragments with material from one chromosome are analyzed, as shown in Fig. 7. Similarly, groups containing fragments with material from two chromosomes are analyzed, as shown in Fig. 8. The groups containing fragments from three or more chromosomes are necessarily complex aberrations.
hTERT-immortalized normal human skin fibroblast cells 82-6 (25), kindly provided by J. Campisi (Lawrence Berkeley National Laboratory, Berkeley, CA), were grown as monolayers at 37°C under 5% CO2 in DMEM (Gibco®/Invitrogen™, Carlsbad, CA) supplemented with 10% fetal bovine serum, antibiotic-antimycotic (1×) and L-glutamine (2 mM). Cells were irradiated in confluent state with accelerated ions, as described in the next section.
Cells were grown in T25 flasks, and irradiated with proton (250 MeV), helium (250 MeV/n), oxygen (350 MeV/n) and titanium (300 MeV/n) ions at NSRL. Dose rates were between 0.5 and 40 cGy/min, depending on the dose delivered. Flasks were filled with media during irradiations.
Premature Chromosome Condensation (PCC)
The PCC technique was used to collect G2/M-phase chromosomes as previously described elsewhere (12, 13, 26). After irradiation, fibroblast cells were allowed to repair for 16 h, and then subcultured at low-cell-population density. After 10-h incubation, cells were arrested in mitosis by adding colcemid to a final concentration of 0.1 µg/ml in the culture media, then cells were incubated for an additional 6 h. Approximately 30 min before collection, 50 nM of calyculin A (Wako Pure Chemical Industries Ltd., Osaka, Japan) was added to the culture media to condense the chromosomes in the G2 phase of the cell cycle.
Fluorescence In Situ Hybridization (FISH)
Chromosomes were hybridized in situ with a combination of three fluorescence whole chromosome human DNA probes for chromosomes 1 (red), 2 (green) and 4 (yellow) (Aquarius®, Cytocell, Oxford Gene Technology, Oxfordshire, UK), using the protocol recommended by the manufacturer. Background chromosomes were stained with 4′,6-diamidino-2-phenylindole (DAPI). Chromosomes were analyzed with Leica CytoVision® FISH system that includes Leica fluorescent microscope with charge-coupled device (CCD) camera and karyotyp-ing software (Nussloch, Germany). The images of all damaged cells were captured electronically. Complex exchanges were scored when it was determined that an exchange involved a minimum of three breaks in two or more chromosomes (27). An exchange was defined as simple if it appeared to involve two breaks in two chromosomes, that is, dicentrics and translocations. Incomplete translocations and incomplete dicentrics were included in the category of simple exchanges, assuming that in most cases the reciprocal fragments were below the level of detection. Each type of exchange (dicentrics, apparently simple reciprocal exchanges, incompletes or complex) was counted as one exchange, and values for total exchanges were derived by adding the yields. When two or more painted chromosomes were damaged, each was scored separately. For each experiment consisting of a single beam at multiple doses, at least 1,000 cells were scored for each data point.
The frequencies of CAs in the painted chromosome(s) were evaluated as the ratio between aberrations scored and total cells analyzed, and extrapolated to whole genome equivalent using a modified version of the Lucas et al. formula (14, 28, 29). Standard errors for aberration frequencies were calculated assuming Poisson statistics. Error bars in each figure represent standard errors of the mean values. The area of the cell nucleus of the 82-6 cells was determined by confocal microscopy with fixed cells at 162.0 µm2.
For comparison of experimental data with modeling results, the background level of CAs present in the cells needed to be subtracted from the results measured for irradiated cells. Consistent with the usual assumption of Poisson statistics, this subtraction was performed by direct Monte Carlo sampling. The average number of CAs from nonirradiated control cells and irradiated cells were obtained in the experiment and assumed to follow Poisson distributions. These distributions were sampled numerically using Monte Carlo methods in an uncorrelated manner, and subtraction was performed on a trial-by-trial basis (i.e., number of CAs sampled from background distribution was subtracted from number of CAs sampled from irradiated distribution). This process formed a new distribution representing the radiation-induced CAs in irradiated cells from which a mean and standard error could be computed. In general, the means obtained from this process were nearly identical to what would be obtained if the average measured values were subtracted directly. The errors were similar to what would be obtained if the measured errors were combined in quadrature. However, this process also allowed us to identify background-subtracted experimental values that were statistically greater than zero (i.e., irradiated cells with CA levels significantly above background levels) at a 95% confidence level. These experimental values are noted on the plots presented later in this article.
RESULTS AND DISCUSSION
Simulation of Radiation Tracks and Nuclei
The first part of the program is the calculation of nuclei chromosomes by RW and stochastic tracks. An example of a fibroblast nucleus irradiated by proton tracks is given in Fig. 9. The nucleus is an ellipsoid of dimensions 17 µm × 17 µm × 6 µm. The chromosomes are different colors, occupying clearly visible domains. The tracks and voxels are confined to the irradiation box encompassing the nucleus, due to periodic boundary conditions.
Fragment Size Distribution
Simulation of chromosome breaks results in the formation of fragments of different length. The distributions of fragment sizes were obtained using RITCARD and NASARTI for a simulation of 16O8+, 350 MeV/n, with a dose of 0.4 Gy and 100,000 histories. The bin size is 1 monomer. Results are shown in Fig. 10.
The vertical lines observed at over 10,000 monomer units correspond to intact chromosomes. The smaller fragments are generated by the chromosome breaks. The fragment size distributions obtained from NASARTI and RITCARD are essentially indistinguishable; therefore, these two codes yield mostly equivalent results at this stage of the simulation.
Repair kinetics curves are obtained by recording the time at which breaks are repaired, and averaging over all simulation histories (Fig. 11). The proper and improper rejoined breaks are recorded separately.
We have performed several tests on the restitution kinetics. As discussed elsewhere (9), the number of remaining breaks as a function of time is often approximated by a bi-exponential decay curve. We have observed that in general, the CA yields are not very sensitive to the shape of the restitution kinetics curve, and that parameters such as the decay times do not influence the results by much. The reason is that the repair is largely dominated by proper repair, which does not lead to aberrations. As the irradiation is delivered at time, t = 0, for this work, most breaks are repaired at 24 h. The most important thing that influences the CA yields appears to be the ratio between the proper and improper repair at 24 h.
The RITRACKS OpenGL interface allows for visualization of CAs, as shown in Fig. 12. The visualization interface clearly illustrates the interaction of tracks leading to a CA, which is a powerful verification and validation tool for the program.
It is also possible to visualize other types of aberrations. However, in many cases, especially for high-LET radiation, many breaks occur in close proximity and the fragments are rather small. For example, as shown in Fig. 5, we observe that most fragments created by the breaks are smaller than 10 monomer units. Since most aberrations involve small fragments, they are difficult to see even under magnification.
The program RITCARD was used to calculate the yield of various types of CAs for a number of ions (1H+, 4He2+, 12C6+, 16O8+, 20Ne10+, 28Si14+, 48Ti22+ and 56Fe26+), with energies varying from 10 to 1,000 MeV/u. The error on the results was estimated using 1.96 times the standard error on the mean, so that they should represent a confidence interval of 95%. Each data point was calculated using 100,000 simulation histories. Selected dose-response curves are shown in Fig. 13. The corresponding experimental data are included in Table 2. In Fig. 13, we observe that the simulation results can be well approximated by a linear-quadratic curve:
Values of a and b Calculated from Fitting of the Simulation Results
In general, experimental results are in good agreement with the simulation results at low dose. At higher dose, the experimental results are lower than the simulation results. Furthermore, it is difficult to infer any conclusion on a trend, notably due to large error bars and lack of points at higher doses. At this time, the cause of this discrepancy is not clear. In the repair process implemented in the code, the rejoining probability is exclusively function of the Euclidian distance, and no restrictions are imposed on the fragment length. It is possible that experiments are unable to detect CAs containing very small fragments, or that allowing CAs to occur with fragments of any length in the simulation model overestimates the number of CAs.
Notably, experimental techniques usually cannot detect all aberrations at once. For example, in the FISH technique, probes are used to paint chromosomes 1 in red, 2 in green and 4 in yellow. All other chromosomes are stained with DAPI. Using this technique, a subset of the exchanges in the sample can be detected, so that the number counted by the experimentalist must be multiplied by a factor to extrapolate the number of exchanges occurring in the entire genome.
Calculation of α and β coefficients
The dose-response curves for all ions simulated were fitted by a linear-quadratic curve [Eq. (6)] in Origin software using the built-in functions. The coefficients α and β and the standard errors on their values were calculated. The coefficients α and β of the fitted curves are shown in Fig. 14. The program clearly shows that there is a peak in the α component at approximately 100 keV/µm. This is an important finding, since such a peak is observed for many measurable end points, including in chromosome aberration simulations (11) and experiments (7, 29, 30). Notably, this peak appears naturally from the calculations without doing specific assumptions in the model.
The results in Fig. 14 (right) show that the quadratic component appears to be less important for high-LET ions. This trend is in line with the analysis made by Loucas et al. (14), as they used purely linear models for high-LET ions. As reported by Loucas et al. (14), for photon irradiation the linear portion is thought to be the consequence of the misrejoining of break pairs formed by ionizations that occur along a single electron track that is set in motion by photo-absorptive events. These events dominate at low doses and low dose rates and are independent of dose rate. At higher doses and dose rates, contribution from the β coefficient begins to dominate and the interaction of breaks caused by independent electron tracks becomes more likely. This process also leads to misrejoining involving two or more breaks, and results in a dose-dependent curvature in the dose response, which is profoundly affected by dose rate. In stark contrast, low-energy charged particles, such as those from naturally emitted alpha particles, have a much higher LET and produce densely spaced ionizations distributed among far fewer tracks to achieve the same dose. For low-LET tracks, the same basic concepts involving the interaction of multiple proximate breaks apply. However, since the average distance between high-LET tracks for a given dose is greater compared to that for low-LET tracks, high-LET track interaction is effectively eliminated for any dose that is likely to be used in an experimental setting. Thus, chromosome exchange formation for high-LET tracks is effectively confined to interactions of breaks caused by single-track action, which is in favor of a linear dose response.
The experimental dose-response curves shown in Fig. 13 do not align well with a linear-quadratic model. However, many experiments (14, 29, 30) and theories (31) support the assumption of a linear-quadratic model for CA. The main difference is that the doses used for these experiments are often an order of magnitude higher than those used in the current work. At low doses, the calculated CA yield is more difficult to obtain due to poor statistics. Furthermore, some CAs are also present in sham-irradiated cells (referred to as background). These background CAs have much more influence on the total CA yield at low doses than at high doses, because the CA yield increases with dose. Furthermore, there is a possibility of non-targeted effects at low dose, especially for high-LET ions. Overall, these contributing factors cause the results to deviate from the linear-quadratic model.
This code is largely based on the previous NASARTI work, which dates back several years. It is possible that a recalibration of the model using newer experimental data is necessary. It could also be useful to review some details included in the model such as DNA break complexity and different repair pathways.
In this work, the new simulation program RITCARD, which is inspired by the code NASARTI, was presented. RITCARD was developed using modern software engineering practices to enable future model advancements (i.e., mixed-field and dose-rate simulations) to be implemented more efficiently and reliably. Although RITCARD and NASARTI are based on the same fundamental theoretical framework, the numerical implementation is vastly different, with the new code (RITCARD) being more reliable, portable, efficient and well documented. By building this code, the structure and algorithms of NASARTI were followed as closely as possible. However, this program and RITRACKS include many improvements as well as new simulation capabilities. For instance, the RITRACKS program now allows multiple ion types to be simulated during irradiation. Amorphous tracks were replaced by stochastic tracks in our CA simulation models. RITCARD is also built to include time dependence in the model, which will allow the study of dose-rate effects on CA in future work. RITCARD uses extensively C++ classes, which are well suited to model and manipulate real-life objects. Due to numerous changes, it was appropriate that a different name be used for this new code. The code is well structured, as shown by a low methods toxicity score. The changes are tracked using the GitLab environment, in accordance with modern software engineering standards. Thanks to its modular structure, the code can be easily modified to include new elements and improvements in the model, notably regarding breaks and repair mechanisms.
The programs were used to perform a systematic study of radiation-induced CA. The dose-response curves of the yields of CA were obtained for several ions of energies varying from 10 to 1,000 MeV/n. In general, the simulated dose-response curves are well approximated by a linear-quadratic model. The coefficient α of this curve appears to peak at approximately 100 keV/µm, suggesting relative biological effectiveness. This peak appears naturally in the calculations, without introducing additional assumptions or parameters in the model. This is very interesting because the model is essentially based on physical calculations. For instance, basically, the radiation-induced chromatin fiber damage is included in Eq. (3), with the voxel dose obtained from the track structure. The experimental data are in relatively good agreement with the simulation results at very low doses (0–0.1 Gy), but the simulation and experimental results diverge as dose increases beyond 0.1 Gy. Several possible explanations have been suggested in the Discussion. A recalibration and possibly a review of the repair kinetics, in which newer experimental data are considered, might be warranted. Thanks to the highly modular structure of RITCARD, it will be relatively easy to implement and test new repair kinetics models without interfering with other parts of the code.
This program will be very useful for the NASA Space Radiation Program. Further use of the code will include calculation of radiation-induced CA for mixed-field and chronic exposures, which are more representative of the GCR environment.
This work was supported by NASA Space Radiation Program (grant no. NNX16AR97G to MH), NASA Human Health and Performance Contract (HHPC; no. NNJ15HK11B), and TAMU Chancellor's Research Initiative, Texas A&M University.
Appendix 1: Files Generated
This appendix provides a description of the files generated by the program RITCARD, with examples (Figs. A1–A9). As counters in C++ usually start with 0 (not 1), many first elements in files are numbered with index 0.
RWIntersectionsBreaks.dat. This file (shown in Fig. A1) comprises information about the break locations.
FreeEnds.dat. This file (shown in Fig. A2) comprises information on the free ends.
FragmentsFE.dat. This file (shown in Fig. A3) contains information on the fragments comprising free ends.
FragmentsFESchema.dat. This file (shown in Fig. A4) schematically displays the fragments.
FreeEndPairs.dat. This long file (shown in Fig. A5) comprises all possible free end pairs in the system.
FreeEndPairsRepairs.dat. This file (shown in Fig. A6) is essentially identical to the file FreeEndPairs.dat, but comprising only the free end pairs that were repaired. All other free end pairs were discarded.
FragmentsSequence.dat. This important file file (shown in Fig. A7) contains the information about the fragment restitution after repair.
FragmentsSequenceSchema.dat. This file (shown in Fig. A8) displays the fragment sequence as well as the type of aberrations calculated. This is essentially the file FragmentsSequence.dat in schematic form.
FragmentsGroups.dat. The groups in this file (shown in Fig. A9) are used to classify aberrations.