Rapid MCNP simulation of DNA double strand break (DSB) relative biological effectiveness (RBE) for photons, neutrons, and light ions

Clusters of DNA lesions generated by ionizing radiation are typically composed of one or a few individual DNA lesions (damaged nucleotides) formed within one or two turns of the DNA. Because the initial direct and indirect physiochemical reactions responsible for the creation of a cluster of DNA lesions are localized within volumes not larger than about 10 nm, the MCDS simulates damage induction in a small segment of a DNA molecule uniformly irradiated by monoenergetic charged particles. Damage within that segment is then scaled in direct proportion to the total amount of DNA in a cell, or conversely, reported as the number of clusters per unit absorbed dose per unit length of DNA. In general, the MCDS simulates the induction of damage using a two-step algorithm. First, the expected number of individual lesions (damaged nucleotides) that would occur per gigabase pair (Gbp) of DNA from the delivery of 1 Gy of monoenergetic radiation in the presence of a defined oxygen concentration is randomly distributed across a small DNA segment. Second, lesions within the DNA segment are grouped into one of three mutually exclusive categories of cluster (single strand break (SSB), DSB, and clusters composed only of base damage). Within each major category of damage, the numbers of each type of cluster are also tallied using metrics indicative of cluster complexity, such as the number of lesions forming the cluster. The algorithm used within the MCDS has three adjustable (input) parameters (σ_Sb, f, n_min) that are, as a good first approximation, independent of the charged particle type and kinetic energy and one parameter (n_seg) that depends on particle type and kinetic energy. Here, n_seg is the length of a segment of DNA in base pairs (bp), σ_Sb is the number of individual strand breaks per unit dose per amount of DNA in the cell (not the DNA segment), f is the ratio of base damage to strand breaks ~3 (the number of nucleotides in the DNA segment with base damage is fσ_Sb), and n_min is the minimum length (in bp) of undamaged DNA (in bp) such that adjacent lesions in the same or opposed DNA strands are considered part of different clusters. The ratio of three damaged bases per strand break (f  =  3) is the same as the value used in recent track structure simulations of non-DSB clusters (Watanabe et al 2015). For charged particles with kinetic energy T, the n_seg parameter is computed as (Semenenko and Stewart 2006)

$$\begin{eqnarray}&&{{n}_{\text{seg}}}(x)=149\,200-\frac{123600(x)}{(x+267)}~ ~ ~\text{(base} ~\text{pairs})\\ &&~ ~ ~ ~ ~ ~ ~\text{where} ~ ~x~=~{{\left[1-\frac{1}{1+\frac{T}{{{m}_{0}}{{c}^{2}}}}\right]}^{-1}}.\end{eqnarray} \tag{ 1 }$$

Here, m₀c² is the rest mass of an electron. The specific mathematical form of equation (1) was determined from an empirical fit of the numbers of SSB and DSB from the MCDS (Semenenko and Stewart 2006) to the results from track structure simulations (Nikjoo et al 1994, 1997, 1999, 2001, 2002, Friedland et al 2003, 2005) for electrons, protons and alpha particles.

The conceptual model used to simulate DNA damage in the MCDS is premised on the conventional dogma that ionizing radiation damages the DNA through the direct ionization or excitation of the DNA and through the indirect action of hydroxyl or other radicals formed in close proximity to the DNA. However, there are many reports in the literature of cellular and genetic damage arising from extracellular signals, i.e. the so-called bystander and non-targeted effects of radiation (Mothersill and Seymour 1998, Morgan 2003a, 2003b, Wright and Coates 2006, Wright 2010). Although the MCDS does not explicitly distinguish DNA damage arising from extracellular signaling in adjacent or nearby cells from damage induced from smaller-scale (<5–10 nm) physiochemical mechanisms, the aggregate effects of direct and all indirect DNA damage mechanisms are ultimately embedded (implicit) in the values selected for the σ_Sb and f parameters. Because the default parameters used in the MCDS were primarily selected to reproduce measured data from cell culture and other experiments (Stewart et al 2011 and references therein), the MCDS implicitly accounts for DNA damage arising from intra- and extracellular signals in adjacent and nearby irradiated cells, e.g. the overall levels of DNA damage observed in a uniformly irradiated monolayer cell culture. Although the MCDS has the ability to separate DNA damage arising from the direct and indirect mechanisms through simulations of the effects of radical scavengers (Semenenko and Stewart 2006), the MCDS does not have any means to separate DNA damage arising from extracellular signaling from damage arising through indirect, subcellular mechanisms. DNA damage arising in non-irradiated ('bystander') cells from extracellular signals generated by (released from) irradiated cells are neglected in the MCDS.

Track structure simulations of the induction of DNA damage by ionizing charged particles (Holley and Chatterjee 1996, Nikjoo et al 1997, 1999, 2001, Ballarini et al 2000, Bernhardt et al 2002, Friedland et al 2002, 2003, Campa et al 2005, Alloni et al 2013) begins with a detailed physical model of the interactions of a particle passing through a liquid medium, usually water. Then, the locations of individual energy transfer events (ionization and excitation) along the trajectory of the ion, as well as any secondary charged particles formed by the primary ion, are typically superimposed on a model of the DNA and chromatin. Energy deposition events that directly involve (coincide with) a nucleotide, or the phosphodiester bond linking adjacent nucleotides, are damaged with a certain probability to form individual strand breaks or base damage, i.e. the direct DNA damage mechanism. The formation of strand breaks and base damage through the direct effect is often parameterized by introducing an energy threshold on the order of a few tens of eV (Nikjoo et al 2001). Energy deposition events in the DNA larger than the threshold produce a DNA lesion (strand break or base damage) whereas smaller energy deposition events do not. Alternatively, a linear or other simple function may be used to relate the size of an energy deposit to lesion induction through the direct mechanism (Holley and Chatterjee 1996, Friedland et al 2003).

To model the formation of individual lesions from the indirect damage mechanism, analog Monte Carlo simulations (Ballarini et al 2000, Kreipl et al 2009, Alloni et al 2013) or simple probability models are used to mimic the diffusion and interactions of hydroxyl (⋅OH) or other radicals with the DNA. Models of the indirect action of radiation typically require the introduction of at least one or two additional parameters, such as a parameter to relate indirect (radical-mediated) hits on a nucleotide to the formation of a strand break and another parameter to relate 'hits' to base damage. For example, a  ⋅OH formed within about 4–6 nm of a nucleotide has a 0.13 chance of creating a strand break (Nikjoo et al 1999, 2001). Similar probability models are used to determine the chance a base is damaged by an  ⋅OH.

An appealing aspect of the track structure simulation approach to modeling DNA damage is that the clustering of two or more individual DNA lesions to form a DSB, or other complex cluster of DNA lesions, arises as a direct consequence of the ionization density (spatial distribution of energy deposits) on a nanometer scale. The parameters and probability models used to relate the spatial distribution of energy deposits to the formation of individual and clusters of DNA lesions are presumably the same for all types of low and high LET radiations. The structure of the track is determined from the physical interactions of the charged particle with the medium and can be considered free of ad hoc adjustable parameters for most purposes. Any effects on the clustering of lesions associated with the higher order structure of the DNA and chromatin arise from the larger-scale features of the track as an ion passes through the cell nucleus.

In the MCDS, the σ_Sb and f parameters quantify the initial numbers of individual lesions formed per unit dose of radiation and per unit length of DNA. These parameters serve the same role in the MCDS as the energy thresholds and probability models used in track structure simulations to relate the magnitude of an energy deposition event in or near the DNA to the formation of individual DNA lesions. As in track structure simulations, the MCDS treats the σ_Sb and f parameters as independent of particle type and LET (kinetic energy). The most significant difference between the MCDS and analog track structure simulations relates to the conceptual model/algorithm used in the placement of individual lesions within the DNA to form clusters of DNA lesions, such as a DSB or complex SSB (e.g. a strand break with base damage in opposed DNA strand). In the MCDS, the placement of individual lesions to form a cluster is indirectly controlled in the MCDS by inserting a fixed number of individual lesions at random into a DNA segment of varying length. The length of the DNA segment depends on the particle type and energy, i.e. the n_seg parameter as defined by equation (1). Large DNA segments correspond to low LET radiations (high-energy charged particles) and small DNA segments correspond to high LET (low-energy charged particles).

The ad hoc seeming method used in the MCDS to control the cluster of lesions within the DNA arose as a fortuitous byproduct of an independent effort to simulate the base and nucleotide excision repair of clusters of DNA lesions other than the DSB (Semenenko and Stewart 2005, Semenenko et al 2005). The accuracy and ultimate outcome from the attempted excision repair of a cluster and DNA lesions (other than a DSB) is critically dependent on the spatial distribution, orientation (same or opposed DNA strand) and the types of lesions (base damage or strand break) within one or two turns of the DNA (~10 to 20 base pairs). The algorithm used for lesion clustering in the MCDS (Semenenko and Stewart 2004) was initially considered no more than a computationally efficient way to generate representative configurations of clusters of DNA lesions to serve as an input into the Monte Carlo excision repair (MCER) model (Semenenko and Stewart 2005, Semenenko et al 2005). However, the simple lesion clustering algorithm used in the MCDS reproduces the numbers of SSB and DSB as well as the relative numbers of eight other categories of simple and complex DNA clusters (SSB+, 2 SSB, DSB+, DSB++, SSB_c, SSB_cb, DSB_c, DSB_cb) with surprising accuracy (Semenenko and Stewart 2004, 2006, Stewart et al 2011). The algorithm used in the MCDS has also been adapted to simulate the effects of radical scavengers (Semenenko and Stewart 2006), the interplay between particle LET and the concentration of oxygen present at the time of irradiation, and the formation of Fpg and Endo III clusters (Stewart et al 2011). We have also shown through combined MCDS  +  MCER simulations that trends in the numbers of HPRT mutations and the numbers of DSB formed through the aborted (incomplete) excision repair of non-DSB clusters are in good agreement with measurements (Semenenko and Stewart 2005, Semenenko et al 2005).

The successes of the MCDS at reproducing measurements and the results of analog track structure simulations may point to a simplified, but perhaps equally fundamental, conceptual model for the way clusters of DNA lesions are formed on a nanometer scale (few turns of the DNA). Namely, the distribution of individual lesions formed within one or two turns of the DNA (i.e. on the spatial scale in which clusters are created) is governed by a stochastic process that is largely determined by the number of DNA lesions formed per unit length of DNA. In the MCDS, the number of lesions per unit length of DNA is determined by the σ_Sb, f and n_seg parameters. The tendency towards increasing levels of lesion clustering as particle LET increases arises because the number of lesions per unit length of DNA increases with increasing LET (σ_Sb and f are constant and n_seg tends to decrease with increasing LET).

For the scheme used in the MCDS to be considered more than an empirical tool to reproduce the results of track structure simulations, all nucleotides within about one or two turns of the DNA (~10 to 20 bp) must be about equally likely to be damaged through the direct and indirect action of ionizing radiation. As conceptual motivation for this approach, consider a linear section of DNA 10 bp in length (i.e. ~length in bp of a DSB composed of exactly two opposed strand breaks) occupying a cylindrical volume not larger than about than 2 nm in diameter and 3 nm in length (~0.3 nm per bp), i.e. ~10 nm³. Because the orientation and packaging of this small segment of DNA is independent of the (random) trajectory of an ion (and the secondary electrons created by the primary), it is reasonable to assume that all 20 nucleotides (=10 bp) in a 10 nm³ are about equally likely to sustain damage from the direct mechanism. As a first approximation, it also seems quite reasonable to assume that nucleotides within a 10 nm³ region of interest are about equally likely to be indirectly damaged because the average distance an  ⋅OH can diffuse in a cellular environment is on the order of 4–6 nm (Roots and Okada 1975, Goodhead 1994, Nikjoo et al 1997, Georgakilas et al 2013), i.e. the same order of magnitude as the dimensions of the molecular domain of interest for cluster formation (~10 nm³). To the extent that all nucleotides within a region of interest associated with the formation of a cluster of DNA lesions (~10 nm³) are about equally likely to be damaged, the MCDS can be considered a simple but still mechanistic alternative to track structure simulations. The σ_Sb and f parameters relate the deposition of energy (dose) to the numbers of individual lesions, and n_seg parameterizes the effects on lesion clustering of differences in track structure and ionization density.

We used the MCDS Version 3.10A (Stewart et al 2011) to generate estimates of the number of DSB Gy⁻¹ Gbp⁻¹ for selected monoenergetic electrons and light ions (i.e. e⁻, ¹H⁺, ⁴He²⁺, ¹²C⁶⁺, and ⁵⁶Fe²⁶⁺) with ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ as large as 2  ×  10⁵. Large values of ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ correspond to low energy, high linear energy transfer (LET) particles. For very high energy ions, the speed $\beta$ of the ion relative to the speed of light approaches unity, and the effective charge approaches Z (Barkas and Evans 1963), i.e.

$$\begin{eqnarray}&&{{z}_{\text{eff}}}=Z\left[1-\exp \left(-125\cdot \beta \cdot {{Z}^{-2/3}}\right)\right],\end{eqnarray} \tag{ 2 }$$

and

$$\begin{eqnarray}&&\beta =\sqrt{1-\frac{1}{{{\left(1+T/{{m}_{0}}{{c}^{2}}\right)}^{2}}}}.\end{eqnarray} \tag{ 3 }$$

Here, T is the kinetic energy of the charged particle (in MeV), and m₀c² is the rest mass energy of the charged particle (in MeV).

All of the reported results, as well as results reported in numerous other publications (Semenenko and Stewart 2006, Reniers et al 2008, Stewart et al 2011, Neshasteh-Riz et al 2012, Wang et al 2012, Hsiao et al 2014, Pater et al 2014, Van den Heuvel 2014, Polster et al 2015), are based on the default MCDS parameter values (σ_Sb  =  217 Gy⁻¹ Gbp⁻¹, f  =  3, n_min  =  9 bp), which have been shown to adequately reproduce cluster yields for a wide range of particle types, kinetic energies and oxygen concentrations (Stewart et al 2011 and references therein). All MCDS simulations are for a minimum of 10 000 cells (standard error of the mean  ⩽0.2%). MCDS simulations were performed for cells irradiated under normoxic conditions (O₂ concentration  ⩾8%) and anoxic conditions (0% O₂). The effects of O₂ concentration on DSB induction are negligible beyond about 8% (figure 2 in Stewart et al 2011).

We computed the RBE for DSB induction relative to ⁶⁰Co γ-rays by dividing the number of DSB Gy⁻¹ Gbp⁻¹ for the ion of interest by the number of DSB Gy⁻¹ Gbp⁻¹ for cells irradiated by ⁶⁰Co γ-rays in vitro (8.32 DSB Gy⁻¹ Gbp⁻¹ normoxic cells, 2.86 DSB Gy⁻¹ Gbp⁻¹ anoxic cells). For comparison to measurements or simulations in which a radiation other than ⁶⁰Co is used as the reference radiation, the estimates of RBE_DSB relative to ⁶⁰Co reported in this work should be multiplied by the DSB yield of ⁶⁰Co and divided by the DSB yield of the desired reference radiation. For example, 125 kVp x-rays produce about 1.11 times as many DSB as ⁶⁰Co γ-rays (Kirkby et al 2013). Therefore, estimates of the RBE reported in this work would need to be divided by 1.11 to use 125 kVp x-rays as the reference radiation. See also table 2 in this work for the RBE of other selected, low LET reference radiations.

The TableCurve 2D Version 5.01 (SYSTAT Software Inc. 2002) software was used to perform an automated regression analysis of the tabulated MCDS estimates of RBE_DSB for the selected ions. Fits to the data were generated for all 3,667 built-in functions and then sorted by the quality of the fit (r²). For normoxic cells, the empirical formula that provides the best fit to the MCDS estimates of RBE_DSB is

$$\begin{eqnarray}&&\text{normoxic} ~\text{cells}:~\text{RB}{{\text{E}}_{\text{DSB}}}=a+b-{{\left\{{{b}^{(1-d)}}+cx(d-1)\right\}}^{\frac{1}{1-d}}},\,2\leqslant x\leqslant {{10}^{5}},\end{eqnarray} \tag{ 4 }$$

where $x\equiv {{\left({{z}_{eff}}/\beta \right)}^{2}}.$ For anoxic cells, the best-fit empirical formula is

$$\begin{eqnarray}&&\text{anoxic} ~\text{cells}:~\text{RB}{{\text{E}}_{\text{DSB}}}=\frac{a+\sqrt{x}\left\{c+\sqrt{x}\left[e+\sqrt{x}\left(g+i\sqrt{x}\right)\right]\right\}}{1+\sqrt{x}\left\{b+\sqrt{x}\left[d+\sqrt{x}\left(\,f+h\sqrt{x}\right)\right]\right\}},\,2\leqslant x\leqslant {{10}^{5}}.\end{eqnarray} \tag{ 5 }$$

For aerobic and anoxic cells when x  <  2, RBE_DSB approaches the asymptotic (low LET, high energy) limit 0.995  ±  0.005. In the low energy (high LET) limit (x  >  10⁵), RBE_DSB  =  3.41  ±  0.01 for normoxic cells and RBE_DSB  =  9.93  ±  0.01 for anoxic cells. Table 1 lists the fitted parameters for equations (4) and (5).

To estimate the overall (average) level of DNA damage at the multi-cellular and tissue levels, Monte Carlo simulations need to correct for the effects of spatial variations in the delivered dose as well as spatial variations in the nature of the radiation field (numbers and types of particles). Formally, the dose-averaged RBE in a volume element of interest (VOI) can be computed by integrating the product ${{D}_{i}}(E)\cdot \text{RB}{{\text{E}}_{i}}(E)$ over all i particle types and kinetic energies (E) and then dividing by the total absorbed dose D, i.e.

$$\begin{eqnarray}&&\text{RBE}\equiv \frac{1}{D}\underset{i}{\mathop \sum}\,\int_{0}^{\infty}\text{d}E{{D}_{i}}(E)\text{RB}{{\text{E}}_{i}}(E).\end{eqnarray} \tag{ 6 }$$

From equation (6), it follows that the RBE-weighted dose computed by the product D  ×  RBE in the VOI is

$$\begin{eqnarray}&&D\cdot \text{RBE}=\underset{i}{\mathop \sum}\,\int_{0}^{\infty}\text{d}E{{D}_{i}}(E)\text{RB}{{\text{E}}_{i}}(E).\end{eqnarray} \tag{ 7 }$$

Although equations (6) and (7) account for DNA damage arising from direct and indirect mechanisms within a VOI, DNA damage arising from extracellular signaling transmitted from one VOI to another are neglected.

MCNP has built-in tallies to record the dose from specific types of particles (e.g. F6:<pl> tallies, where <pl> denotes particle type) as well as from all particles (e.g.  +F6 tally). The standard dose tallies provided by MCNP can also be modified by a user-supplied dose–response function in the form of a data table consisting of an energy-specific data entry (DE) card with the corresponding dose–response (DF) function specified on a second data entry card. RBE_DSB dose–response functions for selected particles (n, e⁻, ¹H⁺, ²H⁺, ³H⁺, ³He²⁺, ⁴He²⁺, ¹²C⁶⁺) in a format suitable for use in MCNP are available for download on the MCDS website (http://faculty.washington.edu/trawets/mcds/). For charged particle transport in MCNP, we used the default Vavilov energy straggling and the finest-allowed energy resolution in stopping powers (efac  =  0.99). Except where explicitly noted otherwise, the default MCNP6 (CEM) physics models were used for neutron interactions.

Figure 1 shows a comparison of RBE_DSB estimates from equations (4) and (5) to estimates obtained directly from the MCDS. The analytic formulas with the parameter estimates listed in table 1 reproduce the MCDS estimates for RBE_DSB within the precision of the Monte Carlo simulations (mean  ±  0.01) for electrons and light ions with Z  ⩽  26 and ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ from 2 (high energy) to 10⁵ (low energy). For ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$  <  2 (low LET charged particles), the RBE_DSB is approximately equal to 0.995  ±  0.005 for cells irradiated under aerobic and anoxic conditions. Values of ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ less than unity are not meaningful because charged particles cannot exceed the speed of light in a vacuum ($\beta$  <  1). In the low energy, high LET limit (i.e. ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$  >  10⁴), estimates of RBE_DSB approach the asymptotic values of 3.41 for normoxic cells and 9.93 for anoxic cells. As illustrated in figure 1, electrons and light ions with the same impact parameter ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ produce the same overall number of DSB Gy⁻¹ Gbp⁻¹ (or per cell). Particles with the same ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ produce about the same number of DSB Gy⁻¹ Gbp⁻¹ because DSB are primarily created in close spatial proximity (~5 nm) to the track core where numerous very low energy (high LET) δ-rays are created. The distance of 5 nm is comparable to the approximate 4 to 6 nm diffusion distance of the hydroxyl radicals usually considered responsible for the indirect action of ionizing radiation on the DNA (Roots and Okada 1975, Goodhead 1994, Nikjoo et al 1997, Georgakilas et al 2013). Although it is tempting to consider extrapolating the trends seen in figure 1 to particles with even larger values of ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$, there is considerable uncertainty in estimates of the relative numbers of DSB induced by massive, high LET (${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$  ⩾  10⁴) particles (e.g. figure 3 in Stewart et al 2011).

Zoom In Zoom Out Reset image size

For large ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$, the asymptotic trends in RBE_DSB apparent in figure 1 can also become inaccurate when the continuous slowing down approximation (CSDA) range of a charged particle is comparable to the dimensions of the cell nucleus (~4–6 μm). The MCDS simulates damage in small (⩽149.2 kbp) segments of DNA given a uniform absorbed dose of the same particle type and energy (Semenenko and Stewart 2006) or a mixture of charged particle types and energies (Stewart et al 2011). The estimate of the number of DNA clusters Gy⁻¹ Gbp⁻¹ formed in the small DNA segment is then scaled in linear fashion to the (optional) user-input DNA content of the cell nucleus. Therefore, the conceptual model of DNA damage induction implemented in the MCDS is most applicable to irradiation geometries in which the entire cell nucleus receives a uniform (expected or average) absorbed dose. Alternatively, the model only applies to the irradiated portion of the DNA within the nucleus that receives a uniform dose. For non-uniform irradiation, the cell nucleus needs to be sub-divided into smaller regions that receive a uniform absorbed dose. Then, damage to each uniformly irradiated sub-region can be simulated in separate runs of the MCDS. Hsiao et al (2014) took this approach when modeling DNA damage induction from Auger-electron emitting radioisotopes distributed within the nucleus.

The interaction of break-ends associated with two different DSB can give rise to lethal and non-lethal chromosome aberrations (Sachs et al 1997b, Hlatky et al 2002, Carlson et al 2008). Also, pairs of DSB in close spatial and temporal proximity are more likely to interact and form a chromosome exchange (Sachs et al 1997a, Carlson et al 2008). The DSB formed by a single ion (or the secondary charged particles produced by this ion) tend to be in closer spatial and temporal proximity than DSB formed by different (independent) light ions. Therefore, the relative numbers of DSB formed per primary particle track through the cell nucleus is a potentially useful, biological metric of the RBE for chromosome damage. Because DSB induction is a linear function of dose up to several hundred or thousands of Gy (Stewart et al 2011 and references therein), the average number of DSB Gbp⁻¹ track⁻¹ can be computed as the product ${{N}_{\text{bp}}}{{\bar{z}}_{\text{F}}}{{\Sigma}_{\gamma}}\text{RB}{{\text{E}}_{\text{DSB}}}$. Here, N_bp is the DNA content of the cell nucleus (e.g. 6 Gbp for a normal diploid human cell), ${{\bar{z}}_{F}}$ is the frequency-mean specific energy (in Gy), and Σ_γ is the number of DSB Gy⁻¹ Gbp⁻¹ produced by the reference radiation (e.g. ⁶⁰Co γ-rays). For particles with a CSDA range that is large compared to the diameter d of a target nucleus, the frequency mean specific energy is well approximated by ${{\bar{z}}_{\text{F}}}=LET/\rho {{d}^{2}}$ (ICRU 1983), where ρ is the mass density of the cell nucleus. When the range of the charged particle becomes comparable to the dimensions of the cell nucleus, ${{\bar{z}}_{\text{F}}}$ is usually less than $LET/\rho {{d}^{2}}$ because of pathlength straggling and because the stopping power is not constant as the particle passes through the cell. A very short-range particle may also lose all of its kinetic energy before passing through the nucleus. For these so-called stoppers, the frequency-mean specific energy is well approximated by ${{\bar{z}}_{\text{F}}}\cong KE/\rho {{d}^{3}}$, where KE is the kinetic energy of the particle incident on the cell nucleus.

The left panel of figure 2 shows the predicted trends in ${{\bar{z}}_{\text{F}}}\cdot \text{RB}{{\text{E}}_{\text{DSB}}}$ (=per track relative number of DSB) as a function of ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ for selected particles (e⁻, ¹H⁺, ⁴He²⁺, ¹²C⁶⁺, ⁵⁶Fe).⁸ The right panel of figure 2 shows an estimate of the average distance between DSB formed by an ion passing through the nucleus of a representative human cell with a cell nucleus 5 μm in diameter (horizontal dashed line in the right panel of figure 2). The average distance between DSB can be computed, under the assumption of random straight line transversals of a particle through the nucleus, by dividing the mean chord length $\bar{l}$ by the number of DSB per track ${{\bar{z}}_{\text{F}}}\text{RB}{{\text{E}}_{\text{DSB}}}{{\Sigma}_{\gamma}}{{N}_{\text{bp}}}$, i.e. $\mu \text{m} ~\text{per} ~\text{DSB} ~\text{=} ~\bar{l}/{{\bar{z}}_{\text{F}}}\text{RB}{{\text{E}}_{\text{DSB}}}{{\Sigma}_{\gamma}}{{N}_{\text{bp}}}.$ The mean chord length for random particles passing through a sphere with diameter d is 2d/3.

Zoom In Zoom Out Reset image size

Although trends in RBE_DSB for electrons and light ions (Z  ⩽  26) are the same when the nucleus is uniformly irradiated (i.e. figure 1), substantial deviations in the per track relative number of DSB (=${{\bar{z}}_{\text{F}}}\cdot \text{RB}{{\text{E}}_{\text{DSB}}}$) are evident as ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ increases. The differences in the per track number of DSB seen in the left panel of figure 2 are primarily due to the finite range of the charged particle (dotted versus solid lines), although some small effects of changes in the stopping power as the particle passes through the target can be seen in the left panel of figure 2 near the peaks in the ${{\bar{z}}_{F}}\cdot \text{RB}{{\text{E}}_{\text{DSB}}}$ curves. In the right panel of figure 2, the average distance between DSB formed by the same track decreases with increasing ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ (and LET) until the CSDA range of the particle is approximately the same as the diameter of the cell nucleus, which is shown as a horizontal dashed line, and then increases with increasing ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$. For electrons with kinetic energies above ~100 eV, the average distance between DSB formed by the same track is larger than a representative cell nucleus 5 μm in diameter. This observation suggests that individual electrons seldom produce more than one DSB while passing through a cell, regardless of kinetic energy. In stark contrast, protons and other light ions up to ⁵⁶Fe²⁶⁺ have the potential to produce multiple DSB in close spatial proximity, which enhances the chance pairs of DSB will be incorrectly rejoined to form a chromosome aberration. Formation of an incomplete exchange-type chromosome aberration may also manifest in assays for DSB as an unrepairable DSB. Because DSBs formed by ionizing radiation need to be correctly repaired, or at least correctly or incorrectly rejoined, for a cell to survive, the downward trends in ${{\bar{z}}_{\text{F}}}\cdot \text{RB}{{\text{E}}_{\text{DSB}}}$ and the upward trend in the average distance between DSB seen in figure 2 for large ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ are a potential explanation for some or all of the trends in particle RBE observed in cell survival experiments (Furusawa et al 2000, Mehnati et al 2005, Czub et al 2008). That is, cell lethality increases with increasing ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ (and LET) up to a particle-specific peak (e.g. 2.3 keV μm for electrons, 57.4 keV μm⁻¹ for protons, and 199.3 keV μm⁻¹ for ⁴He²⁺ ions) because increasing numbers of DSBs (figure 1) are created in closer and closer spatial proximity to each other (figure 2). Beyond this peak, which occurs when the CSDA range of the ion is approximately the same as the diameter of the cell nucleus, cell lethality decreases because the number of DSB Gy⁻¹ Gbp⁻¹ approaches an asymptotic value (figure 1) and the average distance between adjacent DSB formed by the same track increases with increasing ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ (figure 2). The former implies no change in cell lethality per induced DSB with increasing LET, and the latter suggests decreasing lethality per induced DSB with increasing LET. In addition to characteristic changes in the number and average distance between adjacent DSB formed by the same track, the local complexity of the individual DSB, which tends to increase with increasing particle LET (Nikjoo et al 2001, Semenenko and Stewart 2004, 2006, Georgakilas et al 2013), are likely to have a significant impact on the regulation of alternate DSB repair pathways and on the overall kinetics and fidelity of the DSB rejoining process (Jeggo 2002, Pastwa et al 2003, Frankenberg-Schwager et al 2009, Yajima et al 2013, Averbeck et al 2014).

Figure 3 shows an idealized model of a representative (generic) diagnostic x-ray machine. In the model, a mono-energetic, mono-directional beam of electrons is incident on a W target (19.25 g cm⁻³) at an angle of 45° (figure 3, left panel). The right panel of figure 3 shows an example of a 200 and 250 kVp x-ray energy spectra with 1.5 mm and 1.3 mm of Cu filtration, respectively. To compute the RBE-weighted dose in the cell layer, we used a F6:e tally in MCNP 6.1 modified by a DE DF (RBE_DSB) function for electrons. As an initial test of the MCNP RBE tally, we first used the monolayer cell culture model shown in figure 3 to compute the dose-averaged RBE of ¹³⁷Cs γ-rays and ⁶⁰Co γ-rays incident. As expected, the MCNP 6.1 simulation of ⁶⁰Co γ-rays incident on the bottom of the cell culture dish gives an RBE_DSB of 1.003  ±  0.001 at the depth of maximum dose (~4.5 mm) where electronic equilibrium has been established. However for depths less than the 4.5 mm, RBE_DSB increases with decreasing depth to a maximum of 1.015 at a depth of 0 mm (e.g. skin surface). For ¹³⁷Cs γ-rays, the dose-averaged RBE_DSB (relative to ⁶⁰Co γ-rays) at the depth of maximum dose (~1.8 mm) is 1.0170  ±  0.0001 and increases with decreasing depth to 1.040 at a depth of 0 mm.

Zoom In Zoom Out Reset image size

For ⁶⁰Co and ¹³⁷Cs γ-rays, the dose-averaged value of RBE_DSB remains constant from the depth of maximum dose (1.8 mm ¹³⁷Cs, 4.5 mm ⁶⁰Co) for depths up to at least 2 cm. For shallower depths, RBE_DSB may increase by 2–4% with decreasing depth. To avoid introducing systematic or seemingly random errors into dose–response measurements, experiments should be designed to expose relevant biological targets at depths beyond the depth of maximum dose in water-equivalent media, e.g.   >  1.8 mm ¹³⁷Cs,⩾4.5 mm ⁶⁰Co,⩾1.5 cm for 6 MV x-rays,⩾2 cm for 10 MV x-rays and  ⩾3.5 cm for 18 MV x-rays. To ensure reproducibility, experiments to determine neutron RBE should also be designed with the relevant biological targets at the depth of maximum dose, which increases with increasing neutron energy.

Table 2 lists estimates of the RBE_DSB for cells irradiated by 60 to 250 kV x-rays under normoxic and anoxic conditions. For 200 and 250 kV x-rays with Cu filtration, the dose-averaged RBE_DSB in a 4 μm cell layer is 1.178 and 1.163, respectively. For 200 and 250 kV x-rays, changes in RBE_DSB with depth out to at least 2 cm are negligible (data not shown). As illustrated in table 2, estimates of RBE_DSB are far more sensitive to the amount and type of filtration than to the composition of the anode (Mo, Rh, or W) or applied voltage. For example, RBE_DSB for normoxic cells ranges from 1.097 for 250 kV x-rays with 1.3 mm of W filtration to a maximum of 1.240 without W filtration. For cells irradiated under anoxic conditions, RBE_DSB ranges from 1.157 for 250 kV x-rays with 1.3 mm of W filtration to a maximum of 1.441 without W filtration. Estimates of RBE_DSB are sensitive to the amount of filtration because differences in estimates of RBE_DSB primarily arise from the effects of secondary electrons with kinetic energies below 50 keV. That is, RBE_DSB increases with decreasing electron energy below ~50 keV. For electrons with kinetic energies above 50 keV, RBE_DSB  ≅  1.

As protons interact with tissue or other media and lose kinetic energy, proton stopping power first increases with decreasing energy and then decreases with decreasing energy for very low energy protons. Near the end of expected (average) proton range, fluence also rapidly decreases towards zero. These characteristic changes in proton fluence and stopping power give rise to the well known peak in physical dose with depth, i.e. the Bragg peak. Figure 4 shows the results of a MCNP simulation of integral depth dose (IDD) curve for a 163 MeV pencil beam normally incident upon a 40 cm  ×  40 cm  ×  40 cm water phantom. As illustrated in the left panel of figure 4, the magnitude of the Bragg peak tends to decrease as the cutoff energy for proton transport decreases. The location of the Bragg peak also shifts to slightly larger depths (⩽0.1 mm) as the cutoff energy decreases. Smaller cutoff energies substantially increase the CPU requirements of MCNP simulations but provide improved modeling of the spatial and RBE effects of high LET track ends. For cutoff energies  ⩽5 MeV, the choice of cutoff energy has a negligible impact on the shape of proton depth-dose (physical) curve. However as shown in the right panel of figure 4, estimates of the dose-averaged RBE_DSB increase with decreasing proton cutoff energies down to ~0.5 MeV. That is, the estimate of RBE_DSB at the Bragg peak (vertical dashed lines in figure 4) is ~1.1 for E_cut  =  10 MeV, 1.14 for E_cut  =  5 MeV and 1.2 for E_cut  ⩽  0.5 MeV.

Zoom In Zoom Out Reset image size

Beyond the Bragg peak, estimates of RBE_DSB continue to increase with increasing depth up to about 193  ±  1 mm (~8 mm beyond the Bragg peak) and then begin to decrease. The downward trend in RBE_DSB beyond the Bragg peak arises because of (1) the changes in the relative numbers of low and intermediate energy primary-beam protons with depth and (2) because indirectly ionizing photons and neutrons created proximal ('upstream') to the Bragg peak subsequently create low numbers of secondary charged particles distal ('downstream') to the Bragg peak. For a monoenergetic 163 MeV pencil beam, the maximum value of the RBE_DSB is ~1.55. Because of energy and pathlength straggling, the peak value of RBE_DSB increases as the kinetic energy of the incident proton beam decreases. In general, estimates of the peak value of RBE_DSB are larger for monoenergetic beams than for a beam with a distribution of incident proton energies (average energy  =  energy of monoenergetic beam). For a nominal 163 MeV proton beam, the peak value of RBE_DSB decreases from 1.55 for monoenergetic protons incident on the water phantom to 1.45 for incident protons with a normal distribution (E_avg  =  163) and full-width at half maximum (FWHM) of 3 MeV. As the average beam energy goes from 96 MeV to 226 MeV (FWHM  =  3 MeV), the maximum value of RBE_DSB decreases from 1.56 (E_avg  =  96 MeV) to 1.39 (E_avg  =  226 MeV).

Figure 5 shows a 2D depth-dose profile for a 163 MeV proton pencil beam (5.3 mm spot size) incident on a water phantom. The bottom panel in figure 5 shows the absorbed dose profile, and the top panels illustrate the effects of cutoff energy and oxygen concentration on the size and intensity of the biological hotspot (dose  ×  RBE_DSB) near the Bragg peak. The dose in all panels is normalized to the maximum physical dose for a simulation with a 1 keV proton cutoff energy. The size and intensity of the biological hotspot in the Bragg peak increases with decreasing proton cutoff energy. Increases in the width of the biological hotspot indicate the need to account for changes in proton RBE in the lateral beam direction as well as with depth. The magnitude of the biological hotspot is also larger for cells irradiated under anoxic conditions than under normoxic conditions, which suggests proton dose painting to address hypoxia (Flynn et al 2008) may be more advantageous than a plan (outcome) assessment based on physical dose alone (×constant RBE  =  1.1).

Zoom In Zoom Out Reset image size

Figure 6 shows the estimates of RBE_DSB for monolayer cell cultures irradiated under normoxic conditions by monoenergetic neutrons from 1 keV to 100 MeV. To estimate RBE_DSB, we first used the MCDS to tabulate values of RBE_DSB for light ions (¹H⁺, ²H⁺, ³H⁺, ³He²⁺, and ⁴He²⁺) from 1 keV to 1 GeV. Then we used the standard MCNP F6 heating tally, modified by an ion-specific RBE_DSB dose–response function (DE DF card in MCNP), to tally RBE  ×  dose in a 10 μm cell layer (e.g. as in the left panel of figure 3). The dose-averaged value of the RBE_DSB is then computed by summing (RBE  ×  dose)_i over all i light ions (i  =  ¹H⁺, ²H⁺, ³H⁺, ³He²⁺, and ⁴He²⁺) and then dividing by the total dose from the same ions (dashed black line in figure 6). Neutrons and light ions were transported down to the lowest allowed energy in MCNP 6.1 (10⁻¹³ MeV for neutrons and 1 keV for light ions). For neutrons from 10⁻¹³ MeV to 20 MeV, continuous energy cross sections from ENDF/B VIII were used for H (1001.80c and 1002.80c) and O (8016.80c and 8017.80c). From 20 to 100 MeV, we used the default (CEM) MCNP6 neutron physics models for interactions with ²H and ¹⁷O, and the ENDF/B VI continuous energy cross section library for ¹H (1001.66c), and the ENDF/B VIII cross sections for ¹⁶O (1001.80c). Because ions with kinetic energies less than 100 MeV and Z  >  2 have a large ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$, the contributions to the overall neutron RBE of more massive ions (Z  >  2) were estimated by multiplying the dose deposited in the cell layer by these ions by the maximum RBE_DSB  =  3.41 (horizontal dotted line in figure 6). For each energy, simulations are based on 10⁹ source neutrons passing through the cell culture dish.

Zoom In Zoom Out Reset image size

For comparison to the estimates of RBE_DSB from the MCNP simulations, we also estimated RBE_DSB for the special case when neutrons only interact with ¹H through an isotropic (n, p) elastic scattering mechanism (figure 6, blue dashed line). For this special case, the distribution of secondary proton energies produced by neutron elastic scattering is well approximated by a uniform probability distribution with a mean equal to ½ the kinetic energy of the incident neutron (Caswell 1966). As illustrated in figure 6, (n, p) elastic scattering with H is the dominant neutron interaction mechanism in water between about 1 keV and 1 MeV. Below a few hundred keV, the neutron RBE for DSB induction approaches the asymptotic, high ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ limit, of RBE_DSB  =  3.41 (dotted horizontal line in figure 6). Above about 4 MeV, (n, α) and other nuclear reactions, many of which have resonances in the cross sections up to about 30 MeV, have a significant impact on the estimates of neutron RBE.

Table 3 lists dose-averaged estimates of RBE_DSB for selected neutron fission spectra and for the fast neutron Clinical Neutron Therapy System (CNTS) at the University of Washington (Stelzer et al 1994, Kalet et al 2013). For fission neutrons without additional moderation, as well as the higher energy CNTS neutron energy spectrum, 65% to 80% of the dose deposited in cell layer (composed of water) is from light ions (Z  ⩽  2) and the remainder is (primarily) from recoil ¹⁶N⁷⁺ and ¹²C⁶⁺ ions. A 100 MeV ¹²C⁶⁺ ion has a   =  2012 and RBE_DSB  =  2.875, and a 50 MeV ¹²C⁶⁺ ion has a ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$  =  3825 and RBE_DSB  =  3.158. Recoil ¹⁶N⁷⁺ ions have larger values of ${{\left({{z}_{\text{eff}}}/\beta \right)}^{2}}$ and RBE_DSB than ¹²C⁶⁺ ions with the same (corresponding) kinetic energy. As illustrated in table 3, estimates of RBE_DSB for fission, moderated fission, PuBe and the CNTS neutrons differ by 3–7% depending on the method used to account for the effects of ions with Z  >  2. However as a first approximation, the data reported in table 3 suggest that an RBE_DSB  =  2.7  ±  0.1 is appropriate for fission neutrons, CNTS neutrons and PuBe neutrons. For ²⁵²Cf neutrons moderated by 30 cm of heavy water, a slightly higher RBE_DSB of 2.95  ±  0.05 is more appropriate.