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:

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.

FIG. 1

Graphics user interface for entering multiple beam information in RITRACKS. For each contribution, the user must specify the ion type, the energy, the requested dose and the initial and final times of irradiation. Quantities such as the dose, the fluence, the flux, the LET and the dose rate are calculated by the program.

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 r_l 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.

FIG. 2

Schema showing the first part of the program, leading to the formation of breaks, free ends, free end pairs and chromosome fragments.

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.

FIG. 3

Fragments. The numbers on the first line indicate the index of the reactive free end, or if it is a telomere (Tel). In the second line, “o” is a telomere, “x” is a free end, “C” indicates a centromere and “>” represents the direction of the DNA (3′–5′). The numbers in the second line indicate the chromosome number (1–45). Chromosome 0 has three breaks, resulting in four fragments. Chromosomes 18, 24 and 28 have two breaks, resulting in two fragments each. Chromosome 23 has two breaks, resulting in three fragments. The various counts in this file and others generated by RITCARD start at 0, since the first elements of arrays and vectors in C++ are indexed as 0.

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.

FIG. 4

Repair algorithm in RITCARD, illustrating the repair of broken chromosomes (left part) and the formation of fragment sequences (right part).

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).

FIG. 5

Fragment sequences. As in Fig. 3, the numbers on the first line of each sequence indicate the index of the free end, or if it is a telomere (Tel) “o” is a telomere, “x” is a free end, “C” indicates a centromere and “>” represents the direction of the DNA (3–5′). The numbers in the second line indicate the chromosome number.

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.

To determine the CA types it is convenient to group sequences containing material originating from the same chromosomes. The groups for the case illustrated in Fig. 5 are listed in Table 1.

TABLE 1.

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.

FIG. 6

Schema showing an overview of the analysis of the groups of chromosome fragments. As mentioned in the text, the program analyzes the fragment sequences group by group.

FIG. 7

Classification of aberration types in groups containing fragments from one chromosome. In the schema, a break is a type of aberration not considered in the program for now.

FIG. 8

Classification of aberration types in groups containing fragments from two chromosomes.

Cell Culture

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% CO₂ 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.

Irradiation

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 G₂/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 G₂ 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.

Statistical Analysis

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 µm².

Background Subtraction

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.

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.

FIG. 9

Simulation of the irradiation of a fibroblast nucleus by 10-MeV protons, at a dose of 10 cGy. The tracks are shown on the top panel, with each dot representing a radiolytic species created, and the corresponding dose voxels (20 nm in size) are shown on the bottom panel.

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 ¹⁶O⁸⁺, 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.

FIG. 10

Fragment size distribution for a nucleus irradiated by ¹⁶O⁸⁺, 350 MeV/n, at a dose of 0.4 Gy, as calculated by the programs RITCARD (symbols) and NASARTI (line).

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

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.

FIG. 11

Repair kinetics curves for a 1-GeV/n iron ion, 1.2 Gy. The calculations by NASARTI and RITCARD are almost identical.

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.

Chromosome Aberrations

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.

FIG. 12

Illustration of a radiation-induced simple exchange between chromosomes 18 and 23. The nucleus is irradiated by protons, 0.1 Gy. For clarity, only one track is shown. Panels A and B: The intact chromosomes 18 and 23 are shown, with a telomere shown in blue. They are both intersected by a radiation track, which damages the chromosomes (shown in red). Chromosome 18 has one break, whereas chromosome 23 has two. Panels C and D: The fragments are misrejoined, which leads to an exchange.

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 (¹H⁺, ⁴He²⁺, ¹²C⁶⁺, ¹⁶O⁸⁺, ²⁰Ne¹⁰⁺, ²⁸Si¹⁴⁺, ⁴⁸Ti²²⁺ and ⁵⁶Fe²⁶⁺), 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:

FIG. 13

Dose-response curves for simple exchanges for (panels A–D) 250 MeV protons, 250 MeV/n helium, 350 MeV/n oxygen and 300 MeV/n titanium, respectively.

TABLE 2.

Experimental Data

TABLE 3.

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.

FIG. 14

Coefficients alpha (panel A) and beta (panel B) of the dose response curves, as a function of LET, for chromosome exchanges, as calculated by RITCARD.

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.