A closed parameterization of DNA-damage by charged particles as a function of energy - A geometrical approach

Purpose: To present a closed formalism calculating charged particle radiation damage induced in DNA. The formalism is valid for all types of charged particles and due to its closed nature is suited to provide fast conversion of dose to DNA-damage. Methods: The induction of double strand breaks in DNA--strings residing in irradiated cells is quantified using a single particle model. This leads to a proposal to use the cumulative Cauchy distribution to express the mix of high and low LET type damage probability generated by a single particle. A microscopic phenomenological Monte Carlo code is used to fit the parameters of the model as a function of kinetic energy related to the damage to a DNA molecule embedded in a cell. The model is applied for four particles: electrons, protons, alpha--particles, and carbon ions. A geometric interpretation of this observation using the impact ionization mean free path as a quantifier, allows extension of the model to very low energies. Results: The mathematical expression describes the model adequately using a chi--square test ($\chi^2/NDF<1$). This applies to all particle types with an almost perfect fit for protons, while the other particles seem to result in some discrepancies at very low energies. The implementation calculating a strict version of the RBE based on complex damage alone is corroborated by experimental data from the measured RBE. The geometric interpretation generates a unique dimensionless parameter $k$ for each type of charged particle. In addition, it predicts a distribution of DNA damage which is different from the current models.


I. INTRODUCTION
An extensively used model to describe the biological impact of charged particles on cells is the categorization into simple and complex damage introduced by the impacting particles.
The Bethe-formalism describes this impact in terms of energy loss of the charged particles in inelastic collisions with the electrons of the medium, through the notion of mass stopping power (dE/ρdx)[1]. Linear energy transfer (LET) is a measure for the density of ionization taking place along the track of a charged particle through a medium. Hans Bethe showed that stopping power and therefore the closely related LET approximately depends inversely on the speed squared with (i.e. ∼ 1/β 2 with β = v/c). This is corrected by incorporating the effective charge of the charged particle (z ef f ) [2]. It follows that there is a direct relationship of LET with the kinetic energy of the depositing particle. From observation the simple damage was shown to be related to low LET events, while more complex damage is attributed to high LET interactions [3]. Brenner and Ward [4] argued that complex damage was related to single particle interactions with a high LET rather than the combination of simple damage events generated by single particles. In this paper we will use this observation to relate the single interaction probabilities to simple or complex damage.
On the other hand, any molecule, and specifically the DNA molecule, can be described using physical models governed by quantum theory. Therefore, single particle interactions creating complex damage need to follow both classical and quantum models.
To describe the impact of energetic charged particles on the DNA-structure, the science community has taken its recourse to using Monte Carlo simulations to quantify the damage introduced [5,6]. These codes are based on the knowledge provided by the Bethe formalism and are phenomenological in that they are combined with the current experimental knowledge. A more fundamental analytical approach is currently lacking, not in the least due to the underlying complexity of the DNA molecule, as well as the available experimental data which is mainly provided in terms of radiobiologically equivalent dose (RBE), a quantity combining physical, spectral, chemical, and biological factors. All of this hampering ab-initio calculations. Other efforts using more mechanistic models are also available and make use of a classical geometrical view of the DNA-molecule, not taking into account the quantum mechanical nature of the interactions. The Monte Carlo calculations are able to predict the induction of simple damage as well as more complex damage to the DNA-molecules. These findings are interpreted using the Bethe formalism in terms of LET and show that high LET particles introduce more complex damage.
In this paper we will develop an approach that leads to an adequate and universal description of damage induction in DNA-molecules. The approach is based on an overall interaction of the charged particle with the DNA-molecule as a whole. Describing the interaction in terms of quantum mechanical scatter theory. In this approach no assumptions are made regarding the specific structure, nor chemical or biological effects that affect the generation of damage. We also show that this formalism describes the current knowledge well.

A. Molecular Model
In this work we reduce the DNA-molecule to a free-electron molecular orbital model (FEMO), further simplified by assuming a constant potential of the particle in a box, a practice proposed independently by Simpson [7], Bayliss [8], and Kuhn [9]. The potential for all electrons is an effective potential. The model hamiltonian for such a system is thought to be the sum of all individual hamiltonians, We then take this one step further and provide an effective potential V ef f , to serve as single potential interacting with the incoming particle. The approach of replacing a complex system with a many body Schrödinger equation involving an effective potential has also been succesfully used to describe elastic and inelastic interactions if neutrons with nucleons [10,11]. The problem is thereby reduced to the problem of scattering particles by a square well, the size of which is comparable to the size of the DNA-molecule. This is a well known problem in scattering theory and can be solved using the partial wave approach or using the formalism based on the Green's function, both of which result in resonant behavior of the particle with the potential well. Note that there are no additional assumptions regarding the internal structure of the DNA-molecule, nor of any mechanistic process.

B. Resonance formalism
In nuclear physics, a resonance interaction is described by the Breit-Wigner expression (Eq 1) which is completely determined by a position E 0 and a spread Γ. There the parameter and form are well-defined by a theory solving the equations governing the scattering of particles [12,13]. The expression then provides the probability of interaction (or cross section) as a function of energy. In all cases the interaction probability of inelastic scattering is governed by the following equation: (1) This equation provides the probability of inelastic interaction with the DNA-molecule. The physical interpretation of such a resonance is that it describes a phase transition of the interaction of the scattering beam with the target.
In this work we are interested in the case where single particles generate more than one inelastic interaction in the DNA-molecule. We also know that the average distance between interactions in a medium changes with the energy. To model this we propose to use a hybrid model where we combine the quantum mechanical nature of the interaction with a classical average distance between interactions by introducing a forbidden zone in the target potential.
The size of this zone is related to the average distance. Figure 1 illustrates this concept by showing that only interactions (and their probability) confined to that volume can contribute to complex interactions. As the size of the forbidden zone diminishes with the energy of the incoming particle, the contribution of complex damage increases with decreasing energy, thereby sampling an increasingly larger part of the Breit-Wigner distribution. This implies that the probability of complex interaction varies as with energy as the cumulative function of Eq 1. This therefore leads to the ansatz to write the change in probability of the an interaction generating complex damage as the same expression.
The interpretation of this equation is that it provides an estimate of the change in contribution of complex damage to the total damage. Note that the subscript is changed to cd to reflect that we now speak of complex damage instead of multiple interactions, making it clear that this is a mathematical model to represent this rather than an actual mechanism, can be inflicted by a single particle.
It is straightforward that the damage behavior as a function of energy has the following form which is the cumulative expression of the Breit-Wigner function: With the parameters a and b related to the levels 1 and 2 as outlined above. From boundary conditions we find that at very large energies (i.e. E >> E 0 ) the expression is reduced to minimal number of double interactions (D min ) which is equal to a. The value of b is related to the maximal number of double interactions (D max ) as follows:

Validation using Monte Carlo Simulations
The use of microdosimetric calculations has provided important insight in the mechanisms and effects of radiation deposition. In the past Monte Carlo simulations of charged particle deposition by various modalities were used to quantify and typify the kinds of damage introduced by the different modalities [5].
The Monte Carlo Damage Simulation code (MCDS) developed by Semenenko and Stewart, generates spatial maps of the damaged nucleotides forming many types of clustered DNA lesion, including single-strand breaks (SSB), double strand breaks (DSB) as well as individual and clustered base damages [6]. This approach has been shown to be linear with dose up to very high doses, it follows that this parameterization is likewise. In this paper MCDS version 3.0 was used with the parameters described below. The DNA length which was chosen to be 1Gbp (Giga base pairs) and a nucleus diameter of 5µm. In the MCDS software the geometry of the DNA-molecule is not an explicit parameter. Here 4 parameters are used: 1) The DNA-segment length n seg , which is an ad hoc parameter expressed as base pairs Gy −1 cell −1 , 2) The number of strand breaks generated σ Sb , 3) The number of base pair damages generated σ Bb by defining f = σ Bb /σ Sb , and 4) A parameter N min (bp) describing the minimal separation for damage to be apart no to be counted as being in the same cluster. The values of these parameters is determined on the basis of other simulations and measurements. For a more in-depth treatment of these parameters we refer to the work by Semenenko and Stewart [14] Variable input parameters MCDS were: the modality (i.e. energy depositing particle (electron, proton,. . . )), the energy (in MeV), and the oxygen concentration in %. In the implementation described here we chose to omit any oxygen enhancement as this could be a confounding factor and is subject of another study.
Therefore, a concentration of 0% oxygen was used. For every particle type all complex damage was noted per Gy, per cell and per kinetic energy. To obtain differential information a normalized subtraction curve was generated. The fitting procedures were performed in the gnuplot[15]-software using a Levenberg-Marquardt minimization routine.

III. RESULTS
In Figure 2 the Breit-Wigner expression as outlined in Equation 1 together with a normalization factor is used to fit the energy differential cross-section for complex damage.
The fit is performed to minimize the χ 2 -value. In all cases the resulting χ 2 /NDF (NDF = Number of Degrees of Freedom) are lower than 1. The values of the parameters is provided in Table I. All fits are completely satifactory at energies higher than the Breit-Wigner peak.  On the lower energy side different discrepancies can be observed depending on the incoming particles, only for protons we see a satisfactory fit over the full energy range. We attribute this to optical effects with the electron interactions and the change of effective charge in the α-particles and C 6+ -ions as they tend to pick up electrons at lower energies, reducing the charge of the particle and therefore change the interaction parameters. To verify the model we calculate the different parameters starting from the results listed in Table I. Here we convert the values of E 0 and Γ to zp/E which is nothing more than the charge divided by the speed of the particle, which is closely related to the stopping power.
F cd (E) then denotes a response function converting dose to damage.

Energy spectrum -photon treatments
Dose deposition spectra rarely consist of a field of mono-energetic electrons. For a photon source with a given photon spectrum, an energy depositing electron fields exists, which is roughly constant throughout the target volume. Using Monte Carlo simulations it is possible to calculate this field and its spectrum Ψ(E). It then becomes possible to include the spectrum in the calculation of the damage matrices. This approach has been used already by different authors [16,17].
In the case of charged particle treatment, the particles are moderated and the energy spectrum changes depending on the position of the point where the dose is being deposited. It is therefore necessary to apply Equation 6 to each point separately with knowledge of the depositing energy spectrum in that point. Due to the closed nature of this approach this becomes feasible using off the shelf computing equipment.

D. Application: Proton treatment
Recently, the coupling of Monte Carlo simulations in dose deposition to micro-dosimetric code has been proposed and applied by several groups [16,17]. Here a two step approach is followed: 1) A general purpose Monte Carlo code (MCNPX 2.7b) [19] is used to estimate the spectrum of all different dose contributing particles,. 2) A micro dosimetric code [14] is used to determine the biological damage.
The framework for conversion of dose to biological effect is implemented on a simulation of a pristine 200MeV proton beam, taking into account the changing proton spectrum. The proton simulation is performed using MCNPX. Figure 4(a) shows the variation of the number of complex damage evenst as a function of energy of the proton. In addition the spectrum of depositing protons is shown at a position before the Bragg peak and at the Bragg peak. In Figure 4(b) the effect on the dose deposition is shown together with the RBE cd calculated as the complex damage yield generated by the protons at that particular position divided by the complex damage induced by a 6MV photon beam with the same spatial characteristics.
Note that the RBE cd is of the order of 1.1 with larger value of 2 a few mm distal from the Bragg peak. This is commensurate with the cell data reported by Paganetti et al. [20].

IV. DISCUSSION
It is clear that the model used here is over-simplified and is a tool rather than a description of reality. The choice of the Breit-Wigner function shape is a convenient guess based on some knowledge on phase transitions in particle physics rather than a rigorous derivation.
As such this is not a theory or mechanistic model of how the interactions actually take place.  In Fig 4(a) the spectrum at the beginning of the Bragg peak (scaled by 0.1) is in still completely in the low LET regimen, while at the end of the Bragg peak a significant part of the dose depositing protons exhibits high LET characteristics. Fig.4(b) shows the RBE cd i (red line) together with the damage induced by a mono-energetic proton beam.
Carlo simulations fail to take these into account [18]. For the heavier multiply charged ions the discrepancies can be attributed to the increased probabilities of electron capture by the ions at lower energies. Something not accounted for in the simulations.

A. Geometrical interpretation
The expression of the Breit-Wigner function 1 is written as the inverse tangent of a ratio of energies. This ratio is dimensionless, allowing us to substitute this with a number of other quantities provided that we keep the dimensions of both numenator and denominator identical. One interesting quantity is the conversion of energy into a inelastic mean free path length (IMFP) of the charged particle in an appropriate medium . Values for the inelastic mean free path length (IMFP) for electrons are well known in the literature, not in the least as they are important in solid state physics and electron microscopy. They can be found in freely available databases for a variety of elements and compounds, even for organic molecules like DNA [21]. For protons values can be found in a publication by Zhen-Yu and colleagues [22]. For particles such as α-particles and carbon-ions the data is more difficult to find. We therefore opt not to use the data for these particles and restrict ourselves to electrons and protons.
In all current microdosimetric codes the Bethe formalism is used, which is valid for higher energies (i.e. above 500eV for electrons). This implies that changes in IMFP, denoted by λ, which impact the damage calculated using these codes also reflect the limitations of the Bethe formalism. From the theory the following expression is used: Note that the parameters in Table II  energies. This is due to the inability of Bethe's theory to take these into account directly.
Using this it is possible to provide a geometric interpretation of the model as outlined schematically in Figure 5. The probability of complex interaction is then related to the size of the volume in which the IMFP fits in a given volume. The relationship is then described as a Breit-Wigner resonance, which is mathematically equivalent to the Cauchy function wich describes the projection of a point source on a line [23].
FIG. 5. A schematic model of a source of charged particles with a given mean free path length (i.e. a given energy), which is comparable with the diameter of the sensitive cell volume. As the angle (θ) of the particles path with respect to the normal to the axis of the structure increases the chance that more than a single event will occur in the volume. This implies that if the angle for which the average length between interactions equals the diameter is smaller more high LET events will be registered. On the right hand side the abstracted version describes the distribution of horizontal distances at which a line segment tilted at a random angle θ cuts the x-axis.
From Figure 5 it follows that the contribution P for a given energy of charged particles to high LET events is proportional to: Where we denote H=λ(E 0 ) and r=λ(Γ/2) Performing the calculation, we obtain: The full function then has the following form: B. Characterizing the charged particle In this work the parameters H and r have this far not been linked to any physical property. An interesting proposition could be to link these to dimensions of the target structure. Indeed, the choice of a cylinder as a geometric representation is not an accident.
It is natural to use the diameter of a DNA-molecule as a measure of the cylinder's diameter.
The length of the cylinder is then related to the maximal distance we allow to classify two damage events as being part of the same cluster of complex damage. Both values can readily be found in the literature and text books [24]. For the most prevalent form of cellular DNA (B-DNA) the values are 3.4nm (i.e. the height of a spiral of 10 base pairs), and 2.37nm as the diameter, yielding a 1.185nm radius. We now define a dimensionless quantity k which is specific to the type of charged particle used. It is clear that this parameters acts as a scaling parameter but also depends on the ratio of both fixed parameters. Equation 12 now reads as follows: Reducing the impact of a charged particle on the induction of complex damage in a DNAmolecule to three parameters. Figure 6 illustrates the use of these parameters and shows that comparable results can be obtained using this formalism. We find values of k=3.68 for electrons and k=3.42 for protons.

C. Extending the model
In the work presented above as well as in the used Monte Carlo simulations, the Bethe formalism together with its flawed approach in the lower energy regions has always been used.
It is well established that the IMFP does not follow the expression outlined in Equation 7, where λ(E) keeps diminishing as the energy diminishes. Indeed, when the energy is lower than 200eV an increase in IMFP is observed. Ziaja et al [18] showed that it is possible to describe this behaviour analytically by extending Equation 7 with a second term as follows: In this equation the parameter t serves as a threshold separating the behavior as described by Bethe from the plasmonic interactions. Using the data provided in the work from Zhen-Yu and colleagues [22] it is straightforward to obtain parameters for the behavior of protons.
These are presented in table III.  To extend our model to incorporate the behavior of very low energy particles is suffices to replace the expression λ(E) by λ(E) Z in equation 13. In Figure 6 the modified curves show the difference with the calculations based on the Bethe formalism only. This also shows that there is an upper limit to the increase in double strand breaks which depends on the type of particle. It is conceiveable that this approach also works for the heavier particles as can be seen when using the IFMP's in water for these (not shown). [1] H. Bethe. Zur theorie des durchgangs schneller korpuskularstrahlen durch materie. Ann.