key: cord-0783988-vnwzhjl2
authors: Friis, Ida; Verkhovtsev, Alexey V.; Solov'yov, Ilia A.; Solov'yov, Andrey V.
title: Lethal DNA damages caused by ion-induced shock waves in cells
date: 2021-03-18
journal: Physical review. E
DOI: 10.1103/physreve.104.054408
sha: a75d6988d1601b9af53adb72f1b08779a91895c3
doc_id: 783988
cord_uid: vnwzhjl2

The elucidation of fundamental mechanisms underlying ion-induced radiation damage of biological systems is crucial for the advancement of radiotherapy with ion beams and for radiation protection in space. The study of ion-induced biodamage using the phenomenon-based MultiScale Approach to the physics of radiation damage with ions (MSA) has led to the prediction of nanoscale shock waves (SW) that are created by ions in the biological medium at the high linear energy transfer (LET). The high-LET regime corresponding to energy losses higher than 1 keV/nm is typical for ions heavier than carbon in biological media at the Bragg peak region. This paper reveals that the thermomechanical stress of the DNA molecule by the ion-induced SW becomes the dominant mechanism of complex DNA damage at the high-LET ion irradiation. Damage of the DNA molecule in water caused by the ion-induced SW is studied by means of reactive molecular dynamics simulations. Five projectile ions (C, O, Si, Ar, and Fe) at the Bragg peak energies are considered. Simulations reveal that Ar and, especially, Fe ions induce multiple bond breakages in a DNA segment containing 20 base pairs. The DNA damage produced in segments of such size leads to complex irreparable lesions in a cell. This makes the SW-induced thermomechanical stress the dominant mechanism of complex DNA damage at the high-LET ion irradiation. A detailed theory for evaluating the DNA damage caused by ions at high-LET is formulated and integrated into the MSA formalism. The theoretical analysis reveals that a single ion hitting a cell nucleus at high-LET is sufficient to produce highly complex, lethal damages to a cell by the SW-induced thermomechanical stress. A good agreement of the calculated cell survival probabilities with experimental data obtained for the cell irradiation with iron ions provides strong experimental evidence of the ion-induced SW effect.

The elucidation of fundamental mechanisms underlying ion-induced radiation damage of biological systems is crucial for advancing radiotherapy with ion beams and for radiation protection in space. The study of ion-induced biodamage using the phenomenon-based MultiScale Approach to the physics of radiation damage with ions (MSA) has led to the prediction of nanoscale shock waves created by ions in a biological medium at the high linear energy transfer (LET). The high-LET regime corresponds to the keV and higher energy losses by ions per nanometer, which is typical for ions heavier than carbon at the Bragg peak region in biological media. This paper reveals that the thermomechanical stress of the DNA molecule caused by the ion-induced shock wave becomes the dominant mechanism of complex DNA damage at the high-LET ion irradiation. Damage of the DNA molecule in water caused by a projectile ion induced shock wave is studied by means of reactive molecular dynamics simulations. Five projectile ions (carbon, oxygen, silicon, argon and iron) at the Bragg peak energies are considered. For the chosen segment of the DNA molecule and the collision geometry, the number of DNA strand breaks is evaluated for each projectile ion as a function of the bond dissociation energy and the distance from the ion's path to the DNA strands. Simulations reveal that argon and especially iron ions induce the breakage of multiple bonds in a DNA double convolution containing 20 DNA base pairs. The DNA damage produced in segments of such size leads to complex irreparable lesions in a cell. This makes the shock wave induced thermomechanical stress the dominant mechanism of complex DNA damage at the high-LET ion irradiation. A detailed theory for evaluating the DNA damage caused by ions at high-LET is formulated and integrated into the MSA formalism. The theoretical analysis reveals that a single ion hitting a cell nucleus at high-LET is sufficient to produce highly complex, lethal damages to a cell by the shock wave induced thermomechanical stress. Accounting for the shock wave induced thermomechanical mechanism of DNA damage provides an explanation for the "overkill" effect observed experimentally in the dependence of cell survival probabilities on the radiation dose delivered with iron ions. This important observation provides strong experimental evidence of the ion-induced shock wave effect and the related mechanism of radiation damage in cells.

Experimental, theoretical and computational studies of radiation-and collision-induced processes with biomolecular systems are highly relevant nowadays in connection with the molecular-level assessment of biological damage induced by ionizing radiation [1] [2] [3] [4] . The scientific interest in obtaining a deeper understanding of radiation damage is motivated by the development of radiotherapy with ion beams [2, [5] [6] [7] and other applications of ions interacting with biological targets, e.g. radiation protection in space [8, 9] . Protons and carbon ions are currently used for cancer treatment, whereas the clinical implementation of other ions like helium and oxygen has been discussed as the next step [10, 11] . Heavier ions can be found in galactic cosmic rays, where such elements as iron are present, being potentially damaging for humans during space missions [8] .

The mechanisms involved in radiation damage at the nanoscale and molecular level are still not entirely understood and are thus a subject of fundamental multidisciplinary research [1] [2] [3] [4] 12] . The interaction of ion beams with biological materials has commonly been studied computationally by means of track-structure Monte Carlo simulations, which enable to follow the trajectory of each projectile, taking into account different physical interactions, such as elastic and inelastic scattering, electron transfer, nuclear fragmentation reactions, etc. [13] . Some Monte Carlo tools have recently included DNA models in the simulations of biodamage and the subsequent biological response (see [14, 15] and references therein). Despite the wide use of the Monte Carlo approach for modeling ion propagation through biological media, it is unable to simulate the dynamics of the molecular medium in the vicinity of ion tracks, thus missing important physical phenomena.

It has been shown in recent years that a detailed physical understanding of the fundamental processes underlying radiation damage is indeed possible due to recent advances in the theoretical methods and experimental tools developed in atomic and molecular physics [2] . The phenomenon-based MultiScale Approach to the physics of radiation damage with ions (MSA) has been formulated and elaborated during the past decade (see [2] [3] [4] 16] and references therein). This approach considers relevant physical, chemical and biological effects taking place on different scales in space, time and energy, and explores their manifestation in the biological damage. The key phenomena and processes treated by the MSA are ion stopping in the medium, production of secondary electrons and free radicals as a result of ionization and excitation of the medium, transport of secondary electrons and reactive molecular species, the interaction of secondary particles with biomolecules, radiation chemistry, thermomechanical effects caused by nanoscale shock waves induced by ions, and the analysis of induced biodamage. The important outcome of the MSA concerns the prediction of cell response to irradiation with ions on the basis of the assessment of complex DNA damage produced by a cascade of the aforementioned processes. The MSA also demonstrated great success in predicting cell survival probabilities as a function of the radiation dose in a wide range of the systems' parameters, including different cell types, ions with different values of linear energy transfer (LET), oxygenation level, as well as different cell repair conditions [3, 4, 17, 18] .

The important physical effect emphasized by the MSA concerns the manifestation of thermomechanical damage and related phenomena (e.g. transport of reactive secondary species) caused by nanoscale shock waves that are created by high-LET ions traversing biological medium [19] . The formation of ion-induced shock waves was predicted theoretically [19] and studied computationally in a series of subsequent papers [20] [21] [22] [23] [24] [25] [26] [27] [28] [29] . This phenomenon arises due to the fact that ions can deposit a large amount of energy on the nanometer scale resulting in the significant heating up the medium in the localized vicinity of ion tracks. The deposition of the energy lost by the ion into the medium occurs as a result of (i) production, transport and stopping of secondary electrons, and (ii) relaxation of the electronic excitation energy of the medium into its vibrational degrees of freedom through the electron-phonon coupling mechanism [30] .

The average kinetic energy of secondary electrons emitted in the vicinity of the Bragg peak is slightly below 40 eV [3] . Electrons of such energy, experiencing both elastic and inelastic collisions, propagate up to 1 − 2 nanometers away from the ion's path within ∼50 fs before they become solvated electrons [31] . The radial distribution of secondary electrons emitted in the vicinity of the Bragg peak, obtained from the solution of the diffusion equation, is in agreement with the outcomes of track-structure Monte Carlo simulations [32, 33] .

The energy lost by electrons in the processes of ionization and excitation of the medium is transferred to its heating (i.e. vibrational excitation of molecules) due to the electron-phonon interaction, enabling the electronic deexcitation of the molecules from the energy levels forbidden for other channels of de-excitation (such as autoionization, fragmentation or Auger processes). As a result, the medium within the cylinder of the ∼1-2 nm radius surrounding the ion's path is heated up rapidly and the pressure inside this cylinder increases by several orders of magnitude (e.g. by a factor of 10 3 for a carbon ion at the Bragg peak [34] ) compared to the pressure in the medium outside the cylinder. High local temperature and pressure around the ion's path initiate a strong cylindrical explosion of the excited medium, resulting in the formation of a shock wave [19] . Note that this effect has been yet unnoticed in the track-structure models based on the Monte Carlo approach although the classical theory of shock waves was established long ago [35, 36] . Note also that the ion-induced shock wave effect has been completely disregarded in the adaptation of the MSA formalism by other groups [37] .

The two possible mechanisms of DNA damage originating from the ion-induced shock wave have been suggested [3, 22] . The shock wave may inflict damage by the thermomechanical stress and induce breakage of covalent bonds in the DNA molecule [20, 22-24, 26, 29] . Besides, the radial collective motion of the medium induced by the shock wave is instrumental in propagating the highly reactive molecular species, such as hydroxyl radicals and solvated electrons, to large radial distances (up to tens of nanometers) and preventing their recombination [25, 31] .

There are several strong evidences of the ion-induced shock wave effect. First of all, as the shock wave spreads out, it becomes weaker and eventually turns into an acoustic wave at large distances from the ion's path. Acoustic waves coming from the Bragg peak region of ions' trajectories were detected experimentally [38] [39] [40] . Second, a similar phenomenon arising on the micrometer scale was observed during irradiation of micron-sized water droplets with intense X-ray femtosecond pulses [41, 42] . Third, theoretical predictions for the radius and pressure on the shock wave front, based on the analytical solution of hydrodynamic equations [19] , were supported by a series of molecular dynamics (MD) simulations [20, 22, 24, 25, 29] . Finally, the inclusion of the shock wave effect in the multiscale scenario of biodamage with ions [3, 4] has enabled to reproduce experimentally measured cell survival probabilities and related radiobiological quantities such as oxygen enhancement ratio [17, 18] .

In the earlier investigations [20, 22, 24] , the DNA damage by ion-induced shock waves was studied by means of classical MD simulations using non-reactive molecular mechanics force fields. In those simulations the potential energy stored in a particular DNA bond was monitored in time as the bond length varied around its equilibrium distance [22, 24] . When the potential energy of the bond exceeded a given threshold value, the bond was considered broken. A more quantitative description of the phenomenon became possible by means of reactive MD simulations that permitted explicit simulation of covalent bond rupture and formation [43] . A recent study [29] presented a detailed computational protocol for modeling the shock wave induced DNA damage by means of the reactive CHARMM (rCHARMM) force field [43] .

In this paper the thermomechanical stress of the DNA molecule caused by the ion-induced shock wave is systematically explored by means of MD simulations with the rCHARMM force field following the aforementioned computational protocol [29] . Several projectile ions ranging from carbon to iron with different LET values corresponding to the Bragg peak region in liquid water are considered. The number of DNA strand breaks occurring in either one or both DNA strands is evaluated for each projectile ion as a function of the bond dissociation energy and the distance from the ion's path to the DNA strands. The simulations reveal that the shock wave induced thermomechanical stress by carbon and oxygen ions causes only a few isolated strand breaks within a DNA double twist containing 20 base pairs. At higher LET values the thermomechanical stress induced by the shock wave becomes the dominant mechanism of DNA damage. This investigation reveals that argon and especially iron ions produce highly complex DNA damage consisting of multiple localized DNA strand breaks.

The quantitative information obtained from the performed MD simulations has been utilized to evaluate (by means of the MSA formalism) the survival probabilities of cells irradiated with ions. It has been established that the shock wave affects the survival probabilities of cells irradiated with carbon ions mainly via the transport of reactive species away from the ion track. The shock wave induced by a single high-LET iron ion hitting a cell nucleus produces, in addition to the transport of reactive species, lethal damage to the cell due to the thermomechanical stress. The accounting for this DNA damage mechanism within the MSA permits explaining the "overkill" effect, which arises when high-LET ions produce more biodamage than needed for the cell inactivation. A good agreement of the calculated cell survival probabilities with experimental data obtained for the cell irradiation with iron ions provides strong experimental evidence for the ion-induced shock wave effect.

After setting up the all-atom model of a DNA molecule a series of reactive MD simulations have been performed while varying several parameters that characterize the interaction of an ion-induced shock wave with the target. The first part of this section describes the essentials of the computational protocol and introduces the different parameters used in the simulations. More details about this protocol are given in the recent study [29] . The second part of this section outlines the essentials of the MSA formalism regarding the evaluation of the number of DNA lesions produced by a projectile ion and the corresponding cell survival probability. The existing MSA formalism [3, 4] is then extended to account for the shock-wave induced thermomechanical stress in the DNA damage caused by ion irradiation. It should be noted that VMD [44] and MBN Studio [45] software have been used in the data analysis and visualization throughout the paper.

A. Setting MD simulations of the DNA system In order to conduct simulations of DNA damage induced by the shock wave the system must first be constructed and undergo an extensive, multi-step equilibration process to correctly introduce the reactive rCHARMM force field [43] and ensure the system's stability before the simulation of the shock wave propagation. The methodology of designing and equilibrating the system was described in detail in our earlier study [29] , and is therefore only briefly recapped below. The investigated molecular system is created by joining together three short DNA segments (PDB-ID 309D [46] ) resulting in a double-stranded DNA molecule containing 30 complementary base pairs. The molecule is placed in a water box padding of 17 nm from the DNA in the x-and y-directions. The coordinate system used in the simulations is illustrated in Fig. 1A . The x-axis of the coordinate system is oriented along the principal axis of inertia of the chosen DNA molecule with the largest moment of inertia at the initial time instance. The ion track is oriented along the z-axis. The y-axis is along the line defining the shortest distance between the ion track and the selected principal axis of inertia. One sodium ion is placed for every phosphate group present in the DNA to ensure a neutral charge of the entire system, resulting in a system with a total of 1,010,994 atoms. The whole system, including the DNA molecule and the water box, was equilibrated at 300 K temperature before the shock wave simulation. After an initial equilibration in NAMD [47] with the standard CHARMM force field [48, 49] , the system was transferred to the MBN Explorer software [50] , where the reactive rCHARMM force field [43] was used for further simulations.

rCHARMM is used to describe interatomic interactions in the C ′ 3 -O, C ′ 4 -C ′ 5 , C ′ 5 -O and P-O bonds in the DNA backbone, which connect the sugar ring of one nucleotide and the phosphate group of an adjacent nucleotide, see The ion track is oriented along the z-axis; the x-axis is oriented along one of the principal axis of inertia of the chosen DNA molecule, and the y-axis is along the line defining the shortest distance between the ion track and the selected principal axis of inertia. The collision parameter dgeo is defined as the displacement of the ion's path along the y-axis with respect to the principal axis of inertia. The specific collision parameters dA and dB are defined as the shortest distances from the ion's path to DNA strand A and strand B, respectively. Panel B illustrates C ′ 3 -O, C ′ 4 -C ′ 5 , C ′ 5 -O and P-O bonds in the DNA sugar-phosphate backbone and the corresponding potential energy curves obtained by means of DFT [29] . Bond dissociation energy, De, defined as the depth of the associated potential energy well of the covalent bond is considered in the simulations as a variable parameter. The values of De determined from the DFT calculations have been scaled by a factor of 2/3, 1/2, 1/3 and 1/6 (see main text for details). [49] which employs a harmonic approximation for the description of covalent interactions (thereby limiting its applicability to small deformations of the molecular system), rCHARMM treats the bonded, angular and dihedral interactions differently [43] , thus permitting an accurate description of the molecular dissociation process in complex molecular systems. The radial dependence of the bonded interactions is described in rCHARMM by means of the Morse potential. The bonded interactions are set to zero for interatomic distances greater than a user-defined cutoff distance, beyond which the bond is considered broken and the molecular topology of the system is changed. Once a bond starts to break, the associated angular and dihedral interactions involving the indicated atoms weaken and eventually disappear when the distance between the atoms reaches a critical value [43] . Once all the associated bonded, angular and dihedral interactions go to zero, they are automatically removed from the molecular topology of the system. The atoms that initially formed the broken bond are then considered unbound, leading to the formation of atoms with dangling bonds. Bond dissociation energies for the indicated bonds in the DNA strands and the cutoff distances for bond breakage/formation have been obtained by quantum chemistry calculations [29] . Note that the C ′ 3 -C ′ 4 bond was not parameterized by the rCHARMM force field because the dissociation energy of the C ′ 3 -C ′ 4 bond is much higher (9.6 eV according to our DFT calculations) than the dissociation energies of the aforementioned bonds (6.3 − 6.9 eV) [29] .

The DNA damage produced by the shock wave is systematically investigated for five different projectile ions propagating along the z-direction by varying the distance from the ion's path to the DNA strands and dissociation energies of the bonds in the DNA backbone. As described in detail in the following subsections, nine different values of the collision parameter d geo and five scaling factors for the bond dissociation energy D e have been considered for each projectile ion. For each simulation setup (i.e. for each ion type, values of d geo and D e ), two independent simulations of approximately 10 ps duration have been carried out. Hence, 450 independent simulations have been performed in total, and the total simulation time exceeded 2.5 million CPU hours.

The shock wave is induced by an energetic ion propagating through the aqueous environment, where the ion loses its energy mainly by electronic excitation and ionization of water molecules. For ions at the Bragg peak energies, ionization events result in the production of low-energy electrons (with the average kinetic energy of about 40 eV) which propagate radially on the nanometer scale away from the ion's path [3] . Theoretical analysis of secondary electron transport revealed [31] that sub-40 eV electrons lose most of their energy by ionizing and exciting molecules of the medium within approximately 1 nanometer from the ion's path in about 50 femtoseconds after the ion's passage through the medium. The electronic excitation energy of the medium is transferred into its vibrational degrees of freedom through the electron-phonon coupling mechanism [30] . The relaxation of the energy deposited in close proximity to the ion track leads to a rapid increase of the temperature and the pressure of the medium around the ion track, resulting in the dynamical response of the medium and the formation of a cylindrical shock wave that propagates radially away from the ion track [19] . In a continuous medium this phenomenon is characterized by the so-called self-similar flow and the discontinuities of pressure and density of the medium at the wave front as follows from the analytical solution of a set of corresponding hydrodynamic and thermodynamic equations [19] .

In the MD simulations, the energy lost by the propagating ion is deposited into the kinetic energy of water molecules located inside a "hot" cylinder of 1 nm radius around the ion's path. The radius of 1 nm is employed for all the ions considered in this study. The equilibrium velocities of all atoms inside the "hot" cylinder are increased by a factor α such that the kinetic energy of these atoms reads as [20] [21] [22] :

Here S e is the LET of the simulated ion, l is the length of the simulation box in the z-direction (parallel to the ion's path), and N is the total number of atoms within the "hot" cylinder. The first term on the right-hand side of Eq. (1) is the kinetic energy of the 1-nm radius cylinder at the equilibrium temperature, T = 300 K, whereas the second term describes the energy loss by the ion as it propagates through the medium. Note that the 1-nm radius for the energy deposition by low-energy secondary electrons was evaluated [22] as the average distance at which secondary electrons lose most of their energy, according to the random walk approximation. The dispersion of the deposited energy due to more energetic secondary electrons (with the kinetic energy above 40 eV) and its impact on the dynamics of the ion-induced shock wave were addressed in the earlier study [51] . It was demonstrated that, for ions at the Bragg peak, accounting for more energetic secondary electrons makes only a small correction to the results obtained for the uniform energy deposition within the cylinder of 1 nm radius around the ion's path. Since the present study is focused on the effects produced by ions at the Bragg peak region, the utilized "hot" cylinder model captures all the relevant phenomena correctly.

DNA damage caused by the ion-induced shock wave is simulated for five different ions varying the distance from the ion's path to the DNA molecule and dissociation energies of bonds in the DNA backbone, see Fig. 1 . The choice of the specific parameters is explained and justified below.

In the simulations each ion propagated along the z-axis orthogonal to the principal axis of inertia of the DNA molecule with the largest moment of inertia at the initial time instance, below called simply the principle axis of inertia. The collision parameter d geo , defined as the displacement of the ion's path with respect to the principal axis of inertia, varied from 0 to 12Å with an increment step of 3Å. The ion's path was considered at the positive and the negative directions along the y-axis resulting in the positive and negative values of d geo . To account for the orientation of DNA strands with respect to the ion's path, the collision parameter was related to the shortest distance to strand A, d A , and the shortest distance to the strand B, d B . As such, an increase of the displacement d geo could result simultaneously in an increased distance to one strand and a decreased distance to the other strand. Geometry of the system is illustrated in Fig. 1A , whereas the values of the considered parameters d geo , d A and d B are listed in Table I .

The number of DNA strand breaks induced by the shock wave impact may depend on the energy required to break covalent bonds. The typical dissociation energy of covalent bonds in the DNA backbone varies from about 3 to 6 eV [52] . The deposition of such an amount of energy into a given bond would most likely lead to its instantaneous rupture. However, it has also been established that the threshold energy for bond dissociation can be several times smaller due to the presence of solvated electrons in the molecular medium surrounding the DNA. For instance, it was shown [53, 54] that attachment of a solvated electron to a DNA molecule decreases the dissociation energy of covalent bonds in the backbone down to ∼1 eV and leads predominantly to cleavage of a phosphodiester bond. In the present study the bond dissociation energy D e is considered as a variable parameter to account for different possible scenarios that happen on the femto-to sub-picosecond time scales preceding the shock wave formation. A detailed analysis of the DNA damage events created by secondary electrons on the indicated time scales is beyond the scope of this study. Dissociation energies for several bonds along the DNA backbone, shown in Fig. 1B , were determined from density functional theory (DFT) calculations [29] . The obtained values are scaled by a factor of 2/3, 1/2, 1/3 and 1/6 to account for the weakening of the bonds, which may happen e.g. upon the attachment of solvated electrons. The resulting bond dissociation energies thus vary from about 1 eV to 6 eV; this range corresponds to the range of values reported in [52] [53] [54] .

The number of shock wave induced DNA strand breaks also depends on the type of ions irradiating the biological target. Carbon ions are presently used as radiation modality in ion-beam cancer treatments [2, [5] [6] [7] , whereas the interaction with heavier ions (up to iron) is particularly relevant for the radiation protection of astronauts during manned space missions [8] . In the present study shock waves induced by five different projectile ions (C 6+ , O 8+ , Si 14+ , Ar 18+ and Fe 26+ ) with energies corresponding to the Bragg peak region in liquid water are analyzed.

The LET S e as a function of projectile's kinetic energy E has been calculated using the analytical MSA model described in detail in earlier studies [3, 55, 56] . The model is based upon the Rudd's formalism [57] which is extended to account for relativistic corrections and an effective charge of the projectile that arises when a bare ion picks off electrons while propagating through a medium. The dependence of LET on E then reads as:

where n is the number density of water molecules in the medium and W is the kinetic energy of ejected electrons. The sum on the right-hand side is taken over all electron shells of the water molecule with I i being the binding energy of the ith electron shell and dσ i /dW the partial single differential ionization cross section of the corresponding shell. Parameters of the analytical MSA model for liquid water are taken from [58] . Although these parameters were originally derived for proton-water interactions, they are also applicable for evaluating the LET of heavier ions, as illustrated below. Solid lines in Figure 2A show the S e (E) dependence for C 6+ , O 8+ , Si 14+ , Ar 18+ and Fe 26+ ions calculated using Eq. (2). The results are compared with the values compiled in the ICRU73 report [59] (open symbols) and the results

The LET for C 6+ , O 8+ , Si 14+ , Ar 18+ and Fe 26+ ions as a function of ion's kinetic energy, calculated using the analytical MSA model [3] (solid lines). The results are compared with the values compiled in the ICRU73 report [59] (open symbols) and the results of Monte Carlo simulations [60] performed using the Geant4-DNA software package [61] (closed symbols). B: Comparison of the calculated LET for carbon ions (thick solid line) with experimental measurements [62] (open triangles) and other theoretical calculations performed using the widely-used SRIM [63] (thin solid line) and CasP [64] (dashed line) codes. of Monte Carlo simulations [60] performed using the Geant4-DNA software package [61] (closed symbols). Figure 2B shows a detailed comparison of the calculated LET for carbon ions with the results of recent experiments [62] (open triangles) as well as with other theoretical calculations performed using the popular SRIM [63] and CasP [64] codes. The S e (E) dependence calculated using Eq. (2) (thick solid line) gives the best agreement with the experimental data for carbon ions [62] in terms of both the position of the Bragg peak and its magnitude. Reportedly, there is no experimentally measured S e (E) dependence for ions heavier than carbon, and the comparison can only be made with the results of other calculations or Monte Carlo simulations. As shown in Fig. 2 , there is some deviation (about 10 − 15% in the Bragg peak region) between the results obtained with different theoretical methods. The results of the present analysis fit nicely into this range of values. Table II lists the values of LET for each ion at the Bragg peak in liquid water and provides the respective ion's kinetic energy. The values from Table II have been used in Eq. (1) to scale the velocities of atoms of the medium located within the "hot" cylinder for the MD simulations of shock wave propagation.

In order to quantify the impact of the shock wave on the transport of reactive molecular species, additional MD simulations of a shock wave propagating in pure water have been performed following the computational protocol described in [29] . The water box dimensions were set to 49.5 nm × 49.5 nm × 8 nm. No DNA molecule or neutralizing ions were included, so that the shock wave propagated in liquid water. Simulations for the shock wave induced by silicon, argon and iron ions were carried out for ∼10 ps, while the simulations for the lighter (carbon and oxygen) ions were performed for ∼30 ps.

The MSA formalism has been developed to describe survival probabilities of cells irradiated with ion beams on the basis of detailed physical understanding of the fundamental processes underlying radiation damage by ions [3, 4, 17] . As described above, all the relevant physical, chemical and biological processes and phenomena are interlinked within the MSA into a unified multiscale scenario of ion-induced biodamage. A comprehensive description of the MSA formalism is presented in [2] [3] [4] . This section outlines the formalism for evaluating the number of lesions of the DNA molecule produced upon its irradiation with ions and the corresponding cell survival probabilities. The case study is focused on the projectile ions in the vicinity of the Bragg peak. The previously developed formalism [2] [3] [4] is extended towards accounting for the DNA lesions produced by the thermomechanical stress imposed on the DNA molecule by the propagating shock wave.

The starting point for this theory is the calculation of N (r, S e ) -the total average number of simple lesions, i.e. single-strand breaks (SSBs), produced in a DNA double convolution (a DNA double twist) located at distance r from the ion's path. This number depends on the ion's type and its LET S e . According to the MSA analysis [2, 3] , N (r, S e ) is equal to

Here

is the number of simple lesions produced by secondary electrons. The function F e (r, S e ) is the number of electrons incident on the DNA segment located at distance r from the ion's path. The quantity Γ e is the average probability of producing a SSB per electron hit. For ions in the Bragg peak region the probability Γ e does not depend on S e since the average kinetic energy of produced secondary electrons is about 40 eV for different ions with different LET values [3] . For ions outside the Bragg peak region the dependence of Γ e on S e should arise as the kinetic energy of produced δ-electrons is LET-dependent. The analysis of this regime goes beyond the scope of the present study. The second term on the right-hand side of Eq. (3),

is the number of lesions produced by free radicals that are uniformly spread over the distances r < R r (S e ) defined by the radius of shock wave propagation. θ(x) on the right-hand side of Eq. (5) is the Heaviside step function. A linear dependence R r ∝ S e was explored in the earlier study [65] , and a conservative estimate R r ≈ 10 nm was derived for carbon ions in the Bragg peak region [3] . In the present paper the R r value for carbon ions is evaluated more precisely on the basis of MD simulations, and the R r values for heavier ions are estimated according to the R r ∝ S e dependence from the analysis of the pressure at the shock wave front, see Sect. IV B. The value N r,0 (S e ) depends on the number of formed free radicals, which in turn is proportional to the number of generated secondary electrons and hence proportional to LET. N r,0 (S e ) depends also on the degree of oxygenation of the medium since the concentration of oxygen dissolved in the medium affects the number of formed radicals and, consequently, the creation of DNA lesions. For carbon ions at the Bragg peak, the value N r,0 = 0.08 for the environment with the normal concentration of oxygen was derived earlier [17] from the comparison of the experimental results [66] for plasmid DNA, dissolved in pure water and in a scavenger-rich solution, and irradiated with carbon ions at the Bragg peak region. A number of cell survival experiments performed at hypoxic conditions were reproduced with the twice smaller value of N r,0 = 0.04 [2, 17] .

The third term on the right-hand side of Eq. (3), N SW (r, S e ), is the number of DNA lesions produced by the thermomechanical stress imposed on the DNA molecule by the propagating shock wave.

The creation of DNA lesions by secondary electrons, free radicals and the shock wave are statistically independent events taking place at different time scales after the ion passage [2, 3] . Therefore, the total average number of simple lesions in a DNA double twist, N (r, S e ), is a cumulative quantity derived by integrating all the events over time. N e (r, S e ) and N r (r, S e ) were worked out earlier within the MSA [2-4, 17, 31] , whereas N SW (r, S e ) is quantified in the present study by means of MD simulations.

Knowing N (r, S e ) at a given distance r and for a given ion's LET S e , one can use the Poisson statistics to calculate probabilities of different independent events. The probability to produce k lesions in a DNA double twist placed at a distance r from the ion track is equal to

A lethal DNA lesion is defined within the MSA framework as one double-strand break (DSB) plus at least two additional single lesions occurring within a DNA double twist [3] . This definition relies on earlier findings [67] [68] [69] that complex DNA damage is irreparable for a cell if the damage occurs in a localized DNA segment, which typically consists of two helical turns containing 20 base pairs. Lesions within the DNA double twist may occur on one DNA strand or be present on both strands. As shown in Fig. 1 each nucleotide in the DNA molecule has four vulnerable covalent bonds in the sugar-phosphate backbone. Therefore the number of such covalent bonds in one strand in a DNA double twist is equal to 80, and the total number of such bonds in both strands in the DNA double twist is 2n = 160.

The total number of events N ν for ν = 0, 1, ..., 2n lesions occurring within the DNA double twist is equal to the number of combinations for ν choices taken out of 2n places:

On the other hand, the number of independent events of k lesions occurring in a single DNA strand of the length n is equal to C k n . Therefore, N ν can be calculated as follows:

This relationship is well-known, see e.g. Eq. (0.156) in [70] . In the case ν = n + 1, n + 2, ..., 2n the number of events N ν is equal to

Here k = ν − n is the minimum number of lesions on a DNA strand if the other strand within the DNA double twist possesses n lesions. Substituting k = k ′ + ν − n in Eq. (9), one derives

Noting that C n−k n = C k n and using Eq. (0.156 (2)) from [70] , one derives the same relationship as in Eq. (8), but now valid for ν = n + 1, n + 2, ..., 2n. This proves that counting of the lesion events occurring on both DNA strands leads to the same result for N ν as given by Eq. (7).

Similarly, the number of events N (1) ν of ν lesions being all located on one strand within the DNA double twist can be calculated as

n! (n − ν)! ν! , ν = 1, 2, ..., n 0 , ν = n + 1, n + 2, ..., 2n .

The number of events N (1) 0

corresponding to the absence of lesions (ν = 0) is naturally equal to one. For ν = 1, 2, ..., n lesions the factor 2 accounts for the two strands within the DNA double twist. The larger number of lesions (ν = n + 1, n + 2, ...2n) will necessarily occur on both DNA strands, thus the corresponding numbers N ν are equal to zero.

One can also calculate the number of events N (2) ν when ν lesions result in at least one DSB within the DNA double twist:

C k n C ν−k n , ν = 2, 3, ..., n (2n)! (2n − ν)! ν! , ν = n + 1, n + 2, ..., 2n .

The numbers N ν , N (1) ν and N

(2) ν from Eqs. (7), (11) and (12) obey the obvious relationship

Knowing N (1) ν and the total number of events for ν lesions, N ν , one derives the probability P (1) ν to create ν SSBs located on one DNA strand within the double twist:

Substituting here N ν and N (1) ν from Eqs. (7) and (11) respectively, one derives

Analogously, the probability P

(2) ν that ν lesions result in at least one DSB within the DNA double twist reads as

Substituting N

(2) ν from Eq. (12) and using Eqs. (13)-(15) one derives

Now following the above introduced criterion for a lethal DNA lesion, one can derive the probability of such event as follows:

Here λ is the probability that a SSB can be converted to a DSB and ν max = 2n. Accounting for λ relies on the experimental findings [71, 72] that the DSBs caused by low-energy electrons with energies higher than ∼5 eV happen in one hit. In that case the subsequent break in the second DNA strand occurs due to the action of debris generated by the first SSB. Following [71, 72] λ is set equal to 0.15 within the MSA framework [3] . The first term on the right-hand side of Eq. (18) describes the sum of probabilities to have all ν (ν = 3, 4, ..., 2n) lesions on one DNA strand with the subsequent conversion of one SSB into a DSB. The second term is the probability of three lesions with at least one DSB among them and the subsequent conversion of one SSB into a DSB, i.e. creating two DSBs. The third term is the sum of probabilities of ν lesions (ν = 4, 5, ..., 2n) with creation of at least one DSB.

After simple algebraic transformations Eq. (18) can be rewritten in the form:

with P 

and rewrite Eq. (19) in the form P l (r, S e ) = λ Γ(ν max + 1, N (r, S e )) ν max ! − e −N (r,Se) 1 + N (r, S e ) + 1 2 N 2 (r, S e ) The dependence of P l on N (r, S e ) calculated according to Eq. (21) is shown in Fig. 3A .

At small LET values when the number of lesions in a DNA double twist N (r, S e ) < ∼ λ ≪ 1, the probability of lethal events P l (r, S e ) is simplified to

If the characteristic number of lesions is much smaller than the total number of bonds in the DNA double twist, ν ≪ ν max = 2n, the probability P (2) ν , Eq. (17), is reduced to

Then one derives from Eq. (21) the following expression P l (r, S e ) ≃ λ Γ(ν max + 1, N (r, S e )) ν max ! − e −N (r,Se) 1 + N (r, S e ) + 1 2 N 2 (r, S e )

Now let us consider the region N (r, S e ) ≫ 1. Using the definition of the function Γ(n + 1, x), Eq. (20), the fact that

at 1 ≪ N (r, S e ) ≪ ν max , and keeping only the leading terms in Eq. (24), one derives

This means that the probability of lethal lesions P l (r, S e ) → 1 within the entire region where 1 ≪ N (r, S e ) ≪ ν max . Knowing P l (r, S e ) one can now calculate the number of lethal events in a cell nucleus traversed by a projectile ion. Equation (19) represents the probability to create a lethal lesion in a DNA double twist located at the distance r from the ion track. Integrating P l (r, S e ) over the area perpendicular to the ion's trajectory and convoluting the result with the number density of DNA double twists in a cell nucleus one derives the average number of lethal lesions per unit length of the ion's trajectory:

Here n s is the number density of DNA double twists in a cell nucleus which is equal to the number of DNA base pairs accommodated in a cell nucleus, N bp , divided by the number of DNA base pairs in one double twist and by the nuclear volume V n [17] ,

The function σ l (S e ) is the cross section of producing lethal DNA damage in a cell nucleus, which depends on LET and the concentration of oxygen in the target. The σ l (S e ) dependence originates from the dependence of N (r, S e ) on LET; this dependence is discussed further below in this section. The number of lethal events in a cell nucleus at a given dose d produced by N ion ions is equal to [3] :

wherez is the average distance traversed by N ion ions through the cell nucleus. The average number of ions hitting the nucleus, N ion , depends on the nucleus area A n , the dose and LET:

where ρ is the mass density of the irradiated medium taken equal to the density of liquid water, ρ = 1 g/cm 3 . The probability of cell survival is given by the probability of zero lethal lesions occurrence [3] . According to the Poisson statistics it is equal to

Substituting N ion into Y l and taking the logarithm of Π surv one obtains

where dN l dxz is the average number of lethal events created by a single ion in a cell nucleus. Now let us analyze the dependence of the cross section of a DNA lethal lesion σ l and the number of lethal lesions in a cell nucleus Y l on LET. First, consider the irradiation with low-LET ions at the Bragg peak region. A representative case for this scenario is irradiation with protons. In this case the number of lesions in a DNA double twist N (r, S e ) ≪ 1 and hence P l (r, S e ) ∼ N 3 (r, S e ) according to Eq. (22) . For protons at the Bragg peak, the number of lesions in a DNA double twist is proportional to LET, N (r, S e ) ∝ S e . This dependence can be explained as follows. The number of lesions created by secondary electrons incident on a DNA double twist, N e (r, S e ), is proportional to the number of such electrons (see Eq. (4)), which in turn is proportional to LET. The number of lesions created by free radicals, N r (r, S e ), is determined by the number of such species, which is proportional to the number of secondary electrons [31] . The shock wave induced by protons at the Bragg peak does not transport free radicals and other reactive species to the distances much larger than the secondary electron propagation range R e . Note, however, that the diffusion of free radicals on the picosecond time scale might be affected by the temperature increase in the vicinity of the ion tracks. As follows from the analysis described below in Sect. IV B, the free radicals propagation range R r is smaller than R e ∼ 1−2 nm in the S e region up to 70-140 keV/µm. Combining Eqs. (22) and (27) and using the N (r, S e ) ∝ S e dependence one obtains that in this case σ l depends on LET as

The number of lethal lesions in a cell nucleus Y l , Eq. (29), thus increases with LET as

The quantity N r (r, S e ) might grow with the growth of LET due to the increase of the SW radius and correspondingly R r . The growth of R r results in lowering the density of free radicals and thus their recombination rate constant. The additional growth of N r (r, S e ) with S e will result in the faster growth of σ l and Y l with increasing S e . Even steeper dependencies of σ l and Y l on LET may arise at higher S e values when the number of lesions in a DNA double twist N (r, S e ) < ∼ 1 due to a steeper dependence of P l (r, S e ) on N (r, S e ), see Fig. 3 . Finally, let us consider the case N (r, S e ) ≫ 1 when the probability P l (r, S e ) → 1. This case describes iron and heavier ions at the Bragg peak. In this case multiple lesions are created by the shock wave induced thermomechanical stress of the DNA double twist within the distance range r < R SW (S e ) from the ion track. The number of lesions produced by the shock wave in the region r < R SW (S e ) is much bigger than the number of lesions produced by secondary electrons and free radicals, i.e. N SW (r, S e ) ≫ N e (r, S e ) and N SW (r, S e ) ≫ N r (r, S e ). As described in detail in Sect. IV A, the number of lesions N SW (r, S e ) has been evaluated from the MD simulations for the five ions with different LET values at the Bragg peak region. The critical distance R SW is analyzed below in Sect. IV C. These results suggest the following stepwise dependence of N (r, S e ) on distance r from the ion track:

Since at large LET values N ≈ N SW ≫ 1 within the range r < R SW (S e ), the probability P l (r, S e ) → 1 at r < R SW (S e ). Then for high-LET irradiation one obtains

In this case the number of lethal events in a cell nucleus at a given dose d produced by N ion ions, Eq. (29), transforms into:

with the probability of cell survival being given by Eq. (31) . As demonstrated below in Section IV C, R SW depends on LET as R SW = b S 

where α = π b 2 n sz A n ρ d .

This means that the number of lethal events in a cell nucleus at a given dose d decreases slowly with high LET, which corresponds to the experimental observations for iron and heavier ions at the Bragg peak region, see Sect. IV D.

The first part of this section presents the results of the reactive MD simulations of the shock wave induced damage occurring in a 30 base pairs long DNA segment introduced in Fig. 1 . The simulations revealed that most bond breaks in the DNA backbone are produced within the central segment consisting of two helical turns and containing 20 base pairs. Therefore, the number of bond breaks in the central DNA double twist, being the target DNA segment considered within the MSA formalism [3] , has been quantified. The shock wave induced dynamics of the liquid water medium is analyzed next to evaluate the range of shock wave propagation and hence the range of shock wave driven propagation of reactive species. The analysis concludes with the evaluation of survival probabilities of cells irradiated with high-LET ions within the MSA formalism. This analysis reveals the significant role of the shock wave induced thermomechanical mechanism of DNA damage in the cell inactivation.

A. Quantification of the number of bond breaks in the DNA double twist MD simulations of the shock wave induced damage of a 30 base pairs long DNA molecule reveal that the projectile ion propagating in close proximity to the principal axis of inertia of the molecule produces significant damage within the central segment containing 20 DNA base pairs. The DNA damage produced in segments of such size may lead to complex irreparable lesions in a cell [3, 17, 67] . The number of bond breaks in the DNA double twist was counted after each completed MD simulation and analyzed as a function of the distance d geo from the ion's path to the principal axis of inertia of the DNA molecule, see Fig. 1A . The results of the performed analysis are shown in Fig. 4 .

The figure shows that the number of bond breaks produced in the DNA double twist by the ion-induced shock wave increases with the LET of a projectile ion (see Table II ). Simulation results obtained for the scaled bond dissociation energies D e /6 ( Fig. 4A) reveal that up to two DNA backbone bonds break due to the shock wave induced by the carbon ion whereas up to 50 bonds may be broken due to the iron ion impact. For every combination of the bond Table III . As the ion passes at larger distances from the main axis of the DNA molecule the average number of bond breaks within the DNA double twist gradually decreases. Figure 4 shows that the shock wave induced thermomechanical stress of the DNA mostly occurs at < ∼ 1 nm from the ion's path for ions lighter than iron. A more systematic and precise analysis of the threshold distance from the ion's path for inducing DNA strand breaks for each projectile ion is possible, but it would require a significantly larger number of additional simulations aiming at decreasing statistical uncertainties and considering larger values of the collision parameter d geo . Such an analysis might be considered in the future.

The spatial distribution of the total number of bond breaks occurring in the DNA double twist has been analyzed as a function of the shortest distances d A and d B (see Fig. 1 ) from the ion's path to DNA strand A and strand B, respectively. Figure 5 shows the results obtained with the scaled bond dissociation energies D e /6. Due to the geometry of the studied system, not all combinations of d A and d B are accessible at the same time; the inaccessible spatial regions in Fig. 5 are marked with dashed lines. The total number of strand breaks in the DNA double twist decreases with the simultaneous increase of the distance from the ion track to both DNA strands. The results shown in Fig. 5 indicate that the number of bond breaks drops sharply when the ion track passes at distances larger than 5Å to both DNA strands. In contrast, the largest number of bond breaks in the sugar-phosphate backbone is observed when the ion passes in close proximity to at least one of the DNA strands. Fig. 1 ), computed for Fe, Ar, Si, O and C ions in the Bragg peak region. The spatial region inaccessible for the given combination of collision parameters is marked with dashed grey lines. Note that the scale of the color bars for Fe, Ar and Si ions is twofold larger than for the O ion and tenfold larger than for the C ion.

The front of the shock wave propagates radially away from the ion's path. The dependence of the radius of the wave front R on time reads as [19] :

where t is the time from the start of the shock wave propagation, S e is the ion's LET, ρ is the density of the unperturbed medium (ρ = 1 g/cm 3 for liquid water) and β = 0.86 is a dimensionless parameter determined in [19] . The position of the wave front calculated for different projectile ions is illustrated in Fig. 6 .

To evaluate the range of shock wave driven propagation of reactive species, a shock wave propagation was simulated in a pure water box with dimensions of 49.5 nm × 49.5 nm × 8.0 nm. The evolution of radial density of water around the tracks of the five projectile ions is shown in Fig. 7 . Water molecules located in the vicinity of the ion's path are transported away from their initial positions, which results in the formation of a cylindrical cavity around the ion's path. The radius of the cavity grows with time up to the values of about 6 nm for carbon and oxygen ions, while the density of water increases at larger distances from the ion's path. Following the mass conservation law, the mass of water molecules transported from the region in the vicinity of the ion track should be equal to the mass of excess water molecules at larger distances from the track. 40) is indicated with dots. The position of the wave front separates the mass transported by the shock wave, that is considered as displaced (behind the wave front, i.e. at r < R) and excess (beyond the wave front, i.e. at r > R). B: Time evolution of the linear mass density, defined as the water mass transported by the carbon ion-induced shock wave and normalized by the unit length of ion's trajectory. The displaced mass calculated behind the wave front is shown by the dashed black line whereas the excess mass beyond the wave front is shown by the solid green line.

The position of the shock wave front at different time instances, determined using Eq. (40) , is depicted in Fig. 8A with symbols for the case of a carbon ion induced shock wave. The mass of all water molecules behind the wave front is calculated and normalized to the unit length of ion's trajectory. The obtained value is compared to the normalized excess mass beyond the wave front. The comparison shown in Fig. 8B illustrates that the linear mass density is indeed conserved within the simulation time range considered.

The simulation of the propagation of the carbon ion induced shock wave reveals that at a certain time instance the shock wave has stopped propagating away from the ion's path and started to move slowly in the inward direction. This happens when the pressure of the shock wave front drops below a certain value determined by the balance of the pressure at the wave front and the water surface tension pressure [27, 65] . Fig. 9A shows the radial position of the maximal density of water as a function of simulation time for the case of the projectile carbon ion. The radial displacement of the maximal density from the ion track axis increases rapidly during the first 15 ps of the simulation, then reaches the maximal value and starts to decrease at later time instances. The analysis shown in Figure 9A suggests that the maximal radial displacement of the density corresponds to the time instance t = 16.7 ps. Note that for t > 20 ps the radial displacement of the maximal density stops decreasing but fluctuates around the value of 21 nm. This behavior is attributed to interference with the outer part of the shock wave front, which reaches the simulation box boundary and gets reflected. The behavior of the system within the simulation time range t ≤ 20 ps is nevertheless physically meaningful as the shock wave has not yet reached the simulation box boundary within this time interval.

According to Eq. (40), the front of the carbon ion induced shock wave propagates by the time t = 16.7 ps to the distance R = 11.9 nm from the ion track. This characteristic distance defines the propagation range of free radicals, R r in Eq. (5), which are transported by the shock wave driven collective flow. To evaluate the range of shock wave propagation for heavier ions one would need to run longer simulations and consider much larger simulation boxes than the one used in the present study. Alternatively, the range of the shock wave propagation induced by high-LET ions can be estimated from the analysis of the pressure on the shock wave front [19] :

where β = 0.86 is a dimensionless parameter determined in [19] , γ = C P /C V is the heat capacity ratio (γ = 1.222 for liquid water), and R is the radius of the shock wave front. Figure 9B shows the pressure induced by the shock wave front generated by the different ions in the Bragg peak region. Colored lines correspond to the results derived using Eq. (41) . A red dot depicts the pressure P = 0.115 GPa at the distance R = 11.9 nm from the ion track, that is the maximal distance of the wave front propagation for a carbon ion. The indicated value of R corresponding to the instant t = 16.7 ps has been evaluated using Eq. (40) . The shock wave propagation in the radial direction away from the ion's path causes cavitation in its wake, leading to the formation of a rarefied cylindrical region [19, 65] . This effect has been observed in MD simulations described in Fig. 7 and Fig. 8A . In the course of the shock wave propagation, the pressure at the shock wave front becomes balanced by the surface tension pressure building up at the border of the wake region. As a result, the growth of the wake region stops and this region shrinks after the instant when the pressure of the wave front becomes equal to the water surface tension pressure on the surface of the wake region. The latter can be estimated as

where ξ is the coefficient of surface tension and R is the distance from the ion track. Using the aforementioned values P = 0.115 GPa and R = 11.9 nm for the carbon ion at the Bragg peak, one obtains the surface tension coefficient ξ = 1.37 N/m. Note that the medium in the vicinity of the shock wave front is far from equilibrium, and the density of the medium is significantly higher than the density of water at ambient conditions. The high water density in the vicinity of the shock wave front and the large amount of energy deposited into the medium explain the large value of the corresponding surface tension coefficient. The dependence of the surface tension pressure on the distance from the ion track, calculated using Eq. (42) , is shown in Fig. 9B by a dashed line. Assuming that ξ depends weakly and smoothly on LET (or does not depend at all) at the pressures balance point one can evaluate the radii R r for different ions. The radii R r define the propagation ranges of free radicals transported by the shock wave induced by different ions. Equating the pressure on the shock wave front, Eq. (41) , and the surface tension pressure, Eq. (42), one obtains a linear dependence of R r on LET:

The calculated values of R r for the five studied ions at the Bragg peak region are summarized in Table IV . The results indicate that for an iron ion the free radicals are transported by the shock wave to the distance of ∼60 nm. This value exceeds by an order of magnitude typical distances that radicals can diffuse in the medium being at the equilibrium during the time corresponding to the duration of formation of the shock wave wake region with the maximum radius. The characteristic range of the shock wave induced thermomechanical damage, R SW , can be evaluated by analyzing the pressure on the shock wave front, Eq. (41), and the corresponding force exerted by the shock wave on covalent bonds in the DNA backbone.

Equation (41) arises as a solution of the hydrodynamic Euler equation, the continuity and entropy conservation equations [19] . As such, it is applicable to the volumes with sizes being much larger than the characteristic size of a single molecule or a molecular bond. Nevertheless, the medium pressure induces forces applied to single molecular bonds. Such forces should be treated as a result of averaging over an ensemble of molecules exposed to the pressure P (r) at a distance r from the ion track.

Let us evaluate forces acting on molecular bonds in the presence of the ion-induced shock wave. For the sake of simplicity, let us consider a molecular bond oriented parallel to the direction of the shock wave propagation. The force stretching the bond is given by:

where πa 2 0 is the transverse area of the molecular bond exposed to the pressure created by the shock wave front, ∂P ∂r is the pressure gradient, and l is a characteristic interatomic distance in the medium on which the pressure gradient is evaluated. The calculations described below are performed using the value a 0 = 0.15 nm corresponding to the van der Waals radius for atoms forming the DNA backbone [73] and l ≈ 0.15 nm, that is a characteristic length of covalent bonds in the DNA backbone. Let us consider potential energy of a DNA backbone bond being under the pressure created by the shock wave front:

The first term on the right-hand side of Eq. (45) is the bond potential energy described by the Morse potential, D e is the bond dissociation energy, r 0 is the equilibrium bond length, and κ defines the steepness of the potential energy curve. These parameters for the C ′ 3 -O, C ′ 4 -C ′ 5 , C ′ 5 -O and P-O bonds in the DNA backbone (see Fig. 1B ) are determined from the potential energy curves obtained by means of DFT [29] and listed in Table V . The second term on the right-hand side of Eq. (45) describes the work against the force F , caused by the pressure gradient on the distance from r 0 to r. The force F should not depend on the interatomic distance r as it originates from the medium and is defined by its properties at given thermodynamic conditions and at any given location of the bond in space. Considering the geometry when the bond is oriented along F one derives

Stretching the DNA backbone bond by the force F results in lowering the energy barrier for bond rupture, see Fig. 10 . The energy barrier height reads as

where ∆r 1 = r 1 − r 0 is the shift of the potential energy minimum with respect to r 0 , and r 2 = r 0 + ∆r 2 is the position of the potential energy maximum, see Fig. 10 .

The threshold value of the external force at which the bond rupture becomes possible depends on the amount of energy accessible for atoms forming the bond at a given temperature T and on the amount of energy deposited into the TABLE V. The bond dissociation energy De, the equilibrium bond length r0 and steepness of the potential energy curve κ for the C ′ 3 -O, C ′ 4 -C ′ 5 , C ′ 5 -O and P-O bonds in the DNA backbone (see Fig. 1B ). These parameters have been determined from the potential energy curves obtained by means of DFT [29] . medium by the projectile ion. The condition for the bond rupture due to the shock wave induced thermomechanical stress of the DNA can be formulated as follows:

The first term on the left-hand side of Eq. (48) is the average energy available for one degree of freedom in a thermodynamic system being at the equilibrium at T = 300 K. The factor 2 arises since both kinetic and potential energies of the bond are equal to 1 2 k B T according to the equipartition theorem. The second term on the left-hand side is the kinetic energy of the relative interatomic motion caused by the shock wave; µ = m i m j /(m i + m j ) is the reduced mass of a pair of atoms i and j. The analysis described below is performed for the C ′ 3 -O, C ′ 4 -C ′ 5 , C ′ 5 -O and P-O bonds in the DNA backbone which are shown in Fig. 1B .

The gradient of the pressure created by the shock wave on the radial distance l results in the variation of the relative velocity between the atoms forming a molecular bond in the DNA backbone. This variation can be calculated as follows:

where v ′ i and v ′ j are atomic velocities and v is the velocity of the medium. Let us denote ∆

Then Eq. (49) can be written as:

Assuming that the molecular bond is oriented parallel to the direction of the shock wave propagation, one derives

Using the Euler equation

where ρ is the density of the medium, Eq. (50) is rewritten as

where n r = r r is the unit vector. Let us assume that the variation of atomic velocities ∆ v ′ ij ∼ kB T µ caused by the thermal motion of atoms in the molecule is much smaller than the variation of the velocity of the medium caused by the shock wave propagation, i.e.

In this case the condition for the bond rupture, Eq. (48), can be written in the form:

The energy barrier ∆E can be easily evaluated. Thus, equating the derivative of U (r), Eq. (46), over r to zero one derives

By solving this equation one obtains the values ∆r 1 = r 1 − r 0 and ∆r 2 = r 2 − r 0 :

is the dimensionless parameter. Substituting Eq. (57) into Eq. (47) and performing simple algebraic transformations, one derives the expression for the energy barrier ∆E in the following form:

Combining Eqs. (40) and (41), one obtains the relationship between the pressure at the shock wave front and the velocity of the medium caused by the shock wave propagation [19] :

Substituting Eqs. (44) , (58) , (59) and (60) into Eq. (55), one derives the following condition for bond rupture:

where is the dimensionless parameter. The left-hand side of Eq. (61) is a function of distance from the ion track, r, and the position of the shock wave front, R, at a given time moment t.

The parametric inequality (61) is fulfilled in the region α ≥ᾱ, where the threshold valueᾱ is determined by equating the left-and right-hand sides of Eq. (61) . The pressure gradient ∂P ∂r in the vicinity of the shock wave front can be related to ∂P ∂R , that is the derivative of the pressure on the shock wave front P (R), Eq. (41), with respect to the shock wave front radius R. As demonstrated earlier for the carbon ion at the Bragg peak [19] and derived in the present study for the heavier ions at the Bragg peak, the following relation applies:

where the proportionality factor ν = 5.95 is independent on ion's LET and time. Differentiating Eq. (41) over R and combining Eqs. (44) and (58), the condition defining the solution of Eq. (61) can be expressed as:

From this expression one derives the threshold distance from the ion track, R SW , below which the bonds in the DNA backbone can be broken by the shock wave imposed thermomechanical stress:

where the pre-factor b reads as

The distance R SW depends on the parameters of a specific covalent bond (D e and κ) and on the ion's LET. Figure 11 shows the dependence of R SW on S e for the five studied ions in the Bragg peak region. Symbols show the R SW values for cleavage of the C Fig. 11 and Table VI indicate that the threshold distance from the ion's path for the bond rupture by the pressure gradient on the shock wave front varies from 0.7 nm for a carbon ion at the Bragg peak to 1.3 nm for an iron ion at the Bragg peak. The estimated R SW values are consistent with the results of MD simulations shown in Fig. 4 . Note however that no strand breaks have been observed in the simulations for carbon and oxygen projectile ions for the default bond dissociation energies D e (Fig. 4C) , which can be attributed to a small number of simulated trajectories and low number of the events. Accounting for different possible orientations of the molecular bonds in the DNA backbone should also increase their average stability. Simulations performed with the scaled bond dissociation energies D e /2 and D e /6 ( Fig. 4B and Fig. 4A , respectively) indicate the formation of bond breaks in the DNA backbone by the carbon-and oxygen-ion induced shock wave within the range of distances from the ion track, which are consistent with the R SW values determined by Eq. (65) and listed in Table VI. D. Shock wave induced DNA lethal damage of cells irradiated with high-LET ions Figure 12 shows the average number of simple lesions per a DNA double twist as a function of radial distance from the ion's path for irradiation with a carbon ion (Fig. 12A) and with an iron ion (Fig. 12B ) in the vicinity of the corresponding Bragg peaks. The average number of simple lesions created by secondary electrons and free radicals (N e (r, S e ) and N r (r, S e )), is calculated according to Eqs. (4) and (5), respectively. The average number of lesions created by the shock wave induced thermomechanical stress of the DNA, N SW , is taken from the MD simulations described in Sect. IV A (see Table III ). The number of breaks corresponds to the bond dissociation energies D e obtained from the DFT calculations [29] .

As follows from the MD simulations (see Fig. 4 and Table III ) the thermomechanical stress by the carbon ion induced shock wave does not produce any lesions within the DNA double twist for the bond dissociation energies D e . In the case of irradiation with a carbon ion, the lesions are created by secondary electrons, free radicals and other reactive species which are spread over the large distance range by the shock wave, see Fig. 12A . This is in agreement with the results of earlier studies [3, 22] which demonstrated that at the values of LET typical for a single carbon ion at the Bragg peak (S e = 830 keV/µm), most of ion-induced DNA damage occurs via the chemical effects involving interactions of DNA molecules with secondary electrons, free radicals, solvated electrons, etc. In contrast, the number of lesions produced by the thermomechanical stress caused by the iron ion induced shock wave outweighs the number of lesions produced by the chemical effects at distances r ≤ R SW from the ion's path, as shown in Fig. 12B . The analysis described below has been performed using the effective radius R SW = 0.4 nm for the iron ion at the Bragg peak, which is lower than the values reported in Sect. IV C. One should stress that the model presented in Sect. IV C gives the maximal values of R SW for the five studied ions at the Bragg peak, corresponding to the ideal orientation of the molecular bond parallel to the direction of the shock wave propagation. Accounting for different orientations of the bonds in the DNA backbone with respect to the direction of a shock wave propagation should lead to lowering the R SW values. The value R SW = 0.4 nm has been obtained by averaging the number of lesions produced due to direct thermomechanical damage by the iron ion-induced shock wave, N SW , over the range of distances from the ion track to the principal axis of the DNA molecule, see Fig. 4C .

On the basis of the non-reactive MD simulations and subsequent estimates for the energy deposited into the DNA backbone bonds, it was concluded earlier [22] that the bond breaking due to the shock wave induced thermomechanical stress becomes dominant for ions heavier than argon, propagating in liquid water. This result is confirmed in the present study by means of reactive MD simulations.

The calculated probabilities P l (r, S e ) of producing the lethal DNA damage in a DNA double twist located at distance r from the ion's path, Eqs. (19) and (21) , are shown in Fig. 13 for carbon and iron ions. In the case of irradiation with iron ions (see Fig. 13B ), accounting for the shock wave induced thermomechanical stress results in a significant increase of the probability of lethal DNA damage within the characteristic distance r ≤ R SW from the ion track. A conservative estimate for the number of bond breaks produced by the iron ion induced shock wave thermomechanical stress, corresponding to the largest bond dissociation energy D e (see Fig. 4 and Table III ), reveals that five or more bond breaks within the DNA double twist are created when the iron ion propagates at distances smaller than R SW = 0.4 nm from the principal axis of inertia of the DNA molecule. The indicated number of breaks exceeds the minimal number of lesions needed to produce the lethal DNA damage, and hence the probability P l (r, S e ) = 1 at r ≤ 0.4 nm from the iron ion track. This means that even a single hit of a cell nucleus by a high-LET ion will be sufficient to inactivate the cell. Figure 14 shows survival probabilities for two human fibroblast cell lines irradiated with carbon ions at high values of LET; the probabilities were evaluated within the MSA using Eqs. (19)- (32) . Lines show the survival curves obtained with accounting for the DNA damage produced by the secondary electrons, free radicals and the shock wave mechanism. For carbon ion irradiation, the shock wave mechanism enhances transport of radicals and thus reduces their fast recombination thereby increasing the damaging effect of projectile ions. However, the direct thermomechanical DNA damage by the shock wave plays a minor role in the case of carbon ion irradiation. One should stress a good agreement of the calculated survival probabilities with experimental data [74, 75] . These calculations were performed using the range of shock wave driven propagation of reactive species, R r = 11.9 nm, which was determined from the reactive MD simulations described in Sect. IV B (see Table IV ). The shock wave mechanism plays even bigger role in producing lethal damage to cells by high-LET ions as demonstrates in Fig. 13B . Figure 15 shows survival probabilities for two normal rodent cells, V79 and CHO, irradiated with high-LET iron ions in the vicinity of the Bragg peak. Solid red lines show the probabilities calculated with accounting for the shock-wave induced thermomechanical damage. These probabilities were calculated within the MSA using the number of lethal lesions Y l , defined by Eq. (37) , and R SW = 0.4 nm, as discussed above. Dashed black lines show the probabilities calculated with accounting for DNA damage produced only by secondary electrons and free radicals. It is apparent that for the irradiation with high-LET ions the shock wave induced thermomechanical stress of the DNA has a significant impact on the cell survival probabilities. If this mechanism is not taken into consideration, the calculated survival probabilities deviate by orders of magnitude from the experimental values [76, 77] . Indeed, according to Eqs. (27)- (32) , the number of lethal lesions Y l produced by secondary electrons and free radicals in a cell nucleus grows with an increase of LET. As a consequence, the slope of cell survival curves would monotonically increase with an increase of LET. This behavior contradicts with experimentally observed phenomenon known as the "overkill" effect, which manifests itself when cells are irradiated with high-LET ions. At higher LET a given dose can be delivered with the smaller number of ions. This increases chances that some cells remain non-targeted, i.e. the cell survival probability should increase. This leads to a less steep dependence of cell survival probability on the deposited dose [7] .

Different approaches have been adopted in existing radiobiological models to account for the "overkill" effect. For instance, empirical saturation corrections due to non-Poisson distribution of lethal lesions in the cell nucleus were introduced in the commonly used LEM and MKM models to describe the radiobiological response to high-LET irradiation [78, 79] . In contrast to other models, the MSA describes quantitatively the "overkill" effect through accounting for the shock wave induced thermomechanical stress of the DNA.

As follows from Eq. (37) , the quantification of the number of lethal lesions produced by the ion-induced shock wave in a cell requires data on nucleus area for a particular cell line. Solid red curves in Figs. 15(A,B) are obtained with the values A n (V79) = 88 µm 2 and A n (CHO) = 127 µm 2 taken, respectively, from the experimental studies [80, 81] . As it was reported by Konishi et al. [81] , the distribution of nucleus areas for the CHO cells is characterized by a rather broad Gaussian-like profile, and the measured nucleus areas varies from about 80 µm 2 up to 160 µm 2 with the average value of 127 µm 2 . The variation of the calculated cell survival probabilities related to the variation of the nucleus size is illustrated in Fig. 15 by the shaded areas. Note also that no data on the experimental uncertainties of the measured cell survival probabilities were provided in the earlier experimental studies [80, 81] . Therefore, the characteristic uncertainties for the cells irradiated at doses up to about 10 Gy have been estimated based on the typical experimental uncertainties arising in such measurements with carbon ions (Fig. 14) . The estimated uncertainties for the iron ion irradiation are shown in Fig. 15 by gray color. One may thus conclude that, within the experimental uncertainties, the calculated survival probabilities for cells irradiated with iron ions are in a very good agreement with the experimental results [76, 77] . This agreement provides a strong experimental evidence for the biodamage effects caused by ion induced shock waves upon irradiation of biological targets with high-LET ions.

The thermomechanical stress of the DNA molecule caused by the ion-induced shock wave was explored using the reactive molecular dynamics (MD) simulations performed by means of high-performance computing. Five projectile ions with different values of LET, ranging from carbon to iron at the Bragg peak energies in liquid water, were considered. The number of bond breaks in the DNA backbone was systematically evaluated for each projectile ion as a function of bond dissociation energy and the distance from the ion's path to the principal axis of inertia of the DNA molecule.

Reactive MD simulations revealed that argon and, especially, iron ions induce rupture of multiple bonds in a DNA double twist containing 20 DNA base pairs. The DNA damage produced in segments of such size lead to complex irreparable lesions in a cell [3, 17, 67] . This makes the thermomechanical stress of the DNA molecule caused by the ion-induced shock wave the dominant mechanism of complex DNA damage at the high-LET ion irradiation. In contrast, the shock wave induced by lighter ions, such as carbon and oxygen, causes only a few isolated bond breaks within a DNA double twist, but plays an important role in the transport of reactive species to larger distances away from the ion track.

A detailed theory for evaluating the DNA damage caused by ions at high-LET was formulated and integrated into the multiscale approach to the physics of radiation damage with ions (MSA). The theoretical analysis revealed that a single high-LET ion hitting a cell nucleus is sufficient to produce highly complex, lethal damages to a cell by the shock wave induced thermomechanical stress. Using the parameters of the ion-induced shock wave propagation in liquid water, obtained numerically from MD simulations, survival probabilities of cells irradiated with high-LET iron ions were evaluated by means of the MSA. Accounting for the shock wave induced thermomechanical mechanism of DNA damage within the MSA provides an explanation for the "overkill" effect observed experimentally in the dependence of cell survival probabilities on the radiation dose delivered with iron ions. A good agreement of the calculated cell survival probabilities with experimental data obtained for the cell irradiation with iron ions provides strong experimental evidence of the ion-induced shock wave effect and the related mechanism of radiation damage in cells.

Radiation Damage in Biomolecular Systems

Nanoscale Insights into Ion-Beam Cancer Therapy

Ion Beam Therapy: Fundamentals, Technology, Clinical Applications

Course of Theoretical Physics

Physics of Shock Waves and High-Temperature Hydrodynamic Phenomena

International Comission of Radiation Units and Measurements, Stopping of Ions Heavier Than Helium

Table of Integrals, Series, and Products

The authors are grateful for financial support from the Lundbeck Foundation, Danish Councils for Independent Research, the Volkswagen Foundation (Lichtenberg Professorship to IAS), the German Research Foundation DFG (Projects no. GRK1885, SFB 1382 and 415716638), and the European Union's Horizon 2020 research and innovation programme -the Radio-NP project (GA 794733) within the H2020-MSCA-IF-2017 call (Individual Fellowship to AVV). Computational resources for the simulations were provided by the DeiC National HPC Center, SDU and the CARL Cluster at the Carl-von-Ossietzky University Oldenburg, which is supported by the DFG and the Ministry for Science and Culture of Lower Saxony.