The role of shock waves on the biodamage induced by ion beam radiation

Abstract

Background

Energetic ions irradiate the human body in different situations such as exposure to galactic cosmic rays and ion beam cancer therapy. Ions propagating through dense medium deposit energy on the nanometre scale, which cause large local doses in their tracks. This makes the mechanisms of the following biodamage different from that induced by photons. A distinctive feature of the scenario of radiation damage with ions is the potential hydrodynamic response of the liquid medium in the form of cylindrical shock waves around each ion's path.

Review

Recent studies of the predicted cylindrical shock waves produced by swift ions propagating in liquid water or biological environments along with their effects on biological media are discussed in this review. The reasons for the onset of cylindrical shock waves in the vicinity of the ion's path are analysed theoretically, following the hydrodynamic approach to their description. Then, the classical molecular dynamics simulations of ion-induced shock waves are presented and discussed.

Conclusion

The effects of shock waves on the scenario of radiation damage at different values of linear energy transfer are discussed in detail including the most recently simulated chemical effects. A number of ideas of experiments, in which the shock waves could be discovered, are presented, and ideas for the future developments in this field are conveyed.

Background

The human body is subjected to energetic ion radiation in a series of situations, which can lead either to health risks to be avoided or to desired therapeutic effects. Hence, radiation is usually referred to as a double-edged sword. Human irradiation with cosmic rays during manned space missions (Cucinotta and Durante 2006; Chancellor et al. 2018) or accidental exposure to radioactivity may lead to the so-called late effects that include a risk of cancer. However, ion beams are also used to treat cancer in the modern technique known as ion beam cancer therapy or hadron therapy (Schardt et al. 2010; Loeffler and Durante 2013; Surdutovich and Solov’yov 2014; Solov’yov 2017), which is gaining more and more popularity as a cancer treatment worldwide.

Ions with a wide range of linear energy transfer values (LET, the average amount of energy imparted to the tissue by an ion per unit path length of its trajectory) can be found in these situations. Depending on their velocity and charge state, their LET can vary from about 1 to 10,000 keV/μm. Ions with LET $\le 1000$ keV/μm, such as protons and carbon ions, are used in therapy. Ions with higher LET can be found among fission products, in cosmic rays, or produced in a laboratory. Depending on their LET, ions can induce different biological consequences, which need to be understood for the purposes of radiotherapy and radiation protection.

Even light ions such as protons feature higher LET than X-rays; as a result, the scenario of radiation damage and its dependence on dose differs from that of X-rays, mainly due to the non-uniform patterns of energy deposition by ions on the nanometre scale. A relative biological effectiveness (RBE) (Schardt et al. 2010) (a ratio of doses of X-rays and ions producing the same biological effect) reflects this difference, and since the RBE is usually larger than unity, it is one of the major assets of ion therapy giving an opportunity to kill tumour cells depositing lesser dose. This adds up to the advantageous depth–dose curve of ion beams (the so-called Bragg curve). Alas, on the other side, the RBE entails further complications to human safety in space exploration.

From 1960s, a series of biophysical models have been developed in order to relate biological outcomes (e.g. cell survival probabilities or RBE) to the physical dose. Popular examples are the local effect model (LEM) (Schardt et al. 2010; Scholz and Kraft 1996) and the microdosimetric kinetic model (MKM) (Hawkins 1996). More comprehensive approaches are pursued by the track-structure community, aspiring to include all the relevant physical, chemical, and biological mechanisms within Monte Carlo simulation packages (Stewart et al. 2015; Friedland et al. 2017). Recently, some biophysical and Monte Carlo approaches have started to explicitly include DNA models in the simulation of biodamage (Lampe et al. 2018; Sakata et al. 2019) and even some level of atomistic detail (Abolfath et al. 2019). Nonetheless, most of biophysical models, being generally based on the use of empirical relations between the physical dose and the cell survival probability obtained from X-rays experiments, are limited in the inclusion of fundamental atomic and molecular data. This limits the predictive power of these models and their extension to non-conventional situations [e.g. different radiation modalities—X-rays, electrons, ions—, use of radiosensitisers such as metallic nanoparticles (Haume et al. 2016), etc.].

In an attempt to avoid these shortcomings, a MultiScale Approach (MSA) to the assessment of radiation damage with ions has been developed during the last decade (Surdutovich and Solov’yov 2014, 2017a; Solov’yov 2017; Solov’yov et al. 2009; Scifoni et al. 2010; Surdutovich et al. 2009, 2011). The MSA has established a phenomenon-based approach, in which the relevant physical (energy deposition and secondary species generation), chemical (creation, propagation and effects of reactive species) and biological (DNA damage induction and repair by enzymes) mechanisms are incorporated and analysed. The MSA, instead of relying on empirical dose–response relations, predicts cell response to irradiation with ions on the basis of the assessment of complex DNA damage produced by the chain of all the above-mentioned events. The MSA has shown a great success in predicting cell survival curves as a function of radiation dose in a wide range of conditions, including different cell types, levels of LET and oxygenation, as well as different cell repair conditions (Verkhovtsev et al. 2016, 2018). One physical effect emphasised by the MSA and unnoticed by other approaches is the phenomenon of ion-induced shock waves predicted in Surdutovich and Solov’yov (2010). This effect will be at the focus of this review and we need to describe the scenario of events following the ion incidence on tissue in order to explain its origin.

The interaction of ions with tissue can be divided into the so-called physical, chemical and biological stages of irradiation. During the physical stage, ions ionise tissue very fast ($10^{-18}$–$10^{-17}$ s), producing large numbers of electrons (proportional to the LET), most of which have kinetic energies in the sub-50 eV-energy range. These (secondary) electrons propagate, lose their energy and thermalise within $\sim 10^{-13}$ s. Since these electrons lose their energy within 1–2 nm from the ion’s path, the energy deposition turns out to be highly non-uniform, with high density in the ions’ track cores. Large concentrations of free radicals and aqueous electrons are produced by $\sim 10^{-12}$ s after ion’s traverse within the tracks, as a result of excited water molecule dissociation and the stopping of the secondary electrons. These chemically reactive species then start to propagate and react, among them and with biomolecules in the medium, giving place to the so-called chemical stage, lasting up to $10^{-9}$–$10^{-6}$ s. The biological stage follows, in which DNA damage and repair and the subsequent biological consequences take place.

This high energy density in the vicinity of ions’ paths could be dissipated by the processes of molecular diffusion or heat conductivity, but these processes are much slower than the pressure buildup that is achieved by $10^{-13}$ s. Thus, a discontinuity in the initial conditions for pressure is formed, that starts propagating in the direction perpendicular to the ion’s path, also featuring discontinuities in density and collective flow velocity. This process occurs at the end of the physical and beginning of the chemical stages ($10^{-12}$–$10^{-11}$ s), thus having an impact on both of them.

The interest in the study of shock or blast wave phenomena, featuring a propagating surface of discontinuity, has an old history, dating back to the 19th century (Krehl 2001, 2009). However, it was during the 20th century, and particularly in the 1940s, when crucial contributions were made by renowned physicists such as Zel’dovich (1942), Bethe (1947), Sedov (1946), von Neumann (1947) or Taylor (1950). As an important milestone of these works, hydrodynamic solutions were derived to analytically describe spherical shock waves, as comprehensively summarised in Zel’dovich and Raiser (1996), Landau and Lifshitz (1987). Such research was mainly focused on macroscopic phenomena related to detonation. However, similar situations can be found on the nanometre scale as a result of energy deposition in ion tracks. An important example concerns the cylindrical shock waves induced by ions in liquid water, predicted in the context of the MSA (Surdutovich and Solov’yov 2010).¹ The discussion of this phenomenon is the object of this review. Similarly, thermo-mechanical phenomena induced by ion beams also happen in solid targets, where their effects are easier to detect and which present important technological applications (Komarov 2003, 2017).

Within the MSA, the production and physical properties of ion-induced shock waves in biological media were first theoretically studied by means of hydrodynamics, adapting the spherical solutions to the cylindrical symmetry relevant in this case (Surdutovich and Solov’yov 2010). However, it was not clear a priori whether the hydrodynamic model could be applied to describe phenomena on the nanometre scale. Moreover, it does not describe the widening of the wave front and smoothing discontinuities in pressure, velocity, and density, and neither it is capable of describing the shock wave interaction with biomolecules. Therefore, classical molecular dynamics simulations were performed in a large set of different systems and geometries (Yakubovich et al. 2011, 2012a, b; de Vera et al. 2016, 2017b, 2018; Surdutovich et al. 2013a, b), and they were instrumental in answering many questions. Both approaches were shown to predict the main features of ion-induced shock waves in mutual agreement (Yakubovich et al. 2011; de Vera et al. 2016, 2017b), which demonstrated the applicability of the hydrodynamic model to the nanometre scale and offered additional support to the theoretical description of the shock wave scenario. More recently, these approaches are being extended by means of more sophisticated reactive (Bottländer et al. 2015; de Vera et al. 2018) and ab initio (Fraile et al. 2019) molecular dynamics simulations, the development of which also allows considering the potential reactivity induced by radiation in the biological medium. All these advances will be analysed in this review.

Turning to the potential biological effects of ion-induced shock waves which can be studied by means of simulations, two main consequences have been pointed out (Surdutovich and Solov’yov 2010; Surdutovich et al. 2013a): the first one is the direct physical damage of DNA molecules by the large instantaneous temperatures and pressures generated in the medium, and the second one is the fast propagation of chemically reactive species by the induced collective flows, which changes the initial conditions for the induced radiation chemistry. These effects, each of them dominating at different values of LET, have been theoretically studied in a series of works (Surdutovich and Solov’yov 2010; Yakubovich et al. 2011, 2012a, b; de Vera et al. 2016, 2017b, 2018; Surdutovich et al. 2013a, b), and have been incorporated into the MSA, where they play an important role. However, as stated above, the complexity of the problem requires a further development of the theoretical studies. Also, despite the number of indirect signatures of the existence of ion-induced shock waves, the direct experimental discovery of this effect is still pending.

The purpose of the present review is to overview and summarise the work done in the last years on the study of shock waves induced by ions in biological media and their biological effects, and, by doing so, identifying current limitations in the models, knowledge gaps and future research directions. First, the origin of ion-induced shock waves and their biological significance will be discussed in the context of the MSA. The time and space evolution of ion’s track-structure will be analysed, which gives place to the initial condition for the development of shock waves, as demonstrated by hydrodynamics. The more detailed study of the problem by means of molecular dynamics simulations will be reviewed, and the studies on possible direct and/or indirect damage of DNA molecules will be revised. On the basis of the current knowledge on the topic, a discussion will be made on the potential experiments that could be conducted to verify these effects. Also, future research directions will be highlighted, identifying some important open questions in the field of radiation biodamage by ion beams and suggesting ways by which these can be tackled.

Theoretical prediction of shock waves within the MultiScale Approach to the radiation damage induced by ion beams

The MultiScale Approach (MSA) to the radiation damage induced by ions has been developed over the last 10 years, in an effort to provide a comprehensive predictive picture of the radiation damage mechanisms, incorporating relevant physical, chemical and biological effects, happening on the different space, time and energy scales involved (Surdutovich and Solov’yov 2014, 2017a; Solov’yov 2017; Solov’yov et al. 2009; Scifoni et al. 2010; Surdutovich et al. 2009, 2011).

The most complete descriptions of the approach can be found in Surdutovich and Solov’yov (2014, 2017a). The MSA is also reviewed in detail in another contribution within this special volume (Multiscale modelling for cancer radiotherapies). It considers the incident ion propagation in the biological medium, the production of secondary electrons in the process of ionisation of molecules of the medium by ions, and the subsequent electron cascade, as well as electron energy thermalisation and production of chemically reactive species. Biological damage is evaluated in terms of complex DNA damage by the produced secondary electrons and reactive chemical species.

As we shall discuss in the following, the analysis of the ion’s track-structure shows that large radial dose gradients arise around each ion’s path as the result of electronic excitations (Surdutovich and Solov’yov 2015), and that the transfer of this energy to the nuclear degrees of freedom of the target by electron–phonon coupling generates a “hot” medium (Toulemonde et al. 2009) which might prompt its hydrodynamic response in the form of cylindrical shock waves (Surdutovich et al. 2011). This will have a series of consequences on DNA damage.

In the next subsections, first the evolution of the ion’s track-structure will be analysed, which puts the initial conditions for shock waves. Their main physical characteristics will be reviewed from the point of view of hydrodynamics. More details on how DNA damage can be evaluated by means of molecular dynamics simulations are given in the next section.

Electron transport around the ion’s path and the onset of ion-induced shock waves

Although initially an ion traversing the biological medium loses its energy by ionisation and electronic excitation of the target’s molecules along its trajectory, this energy will not be deposited in the point of the interaction event. This is due to the transport of secondary electrons, which will give rise to the so-called radial dose $D(\rho )$. This represents the amount of energy (per unit mass of the target) deposited at a given radial distance $\rho$ from the ion’s path. From the radial dose, it is possible to calculate the energy stored inside a cylinder of radius $\rho$ around an ion’s path of length L, which is (Surdutovich and Solov’yov 2015) as follows:

$$\begin{aligned} U(\rho ) = \int _0^{\rho } D(\rho ') 2\pi \rho '\text{d}\rho 'L \varrho \text{, } \end{aligned}$$

(1)

where $\varrho$ is the mass density of the medium. The slope of this radial profile of deposited energy may give rise to a pressure profile, since the force applied to the wall of the cylinder of radius $\rho$ is given by $F(\rho )=-\nabla U(\rho ) = 2\pi \rho L \varrho D(\rho )$. Thus, the pressure dependence on the radial distance is given by

$$\begin{aligned} P(\rho ) = \frac{F(\rho )}{2\pi \rho L} = \varrho D(\rho ) \text{, } \end{aligned}$$

(2)

i.e. it is directly proportional to the radial dose. This is an important relation, since it means that the radial dose will not only be related to the possible damage of biomolecules located in the vicinity of an ion’s path (by secondary electrons and chemically reactive species), but it will also determine the pressure buildup around an ion’s path. If this radial profile of pressure is steep enough, it will produce the hydrodynamic response of the medium, as we will review in the following.

The average energy of secondary electrons slowly increases with increasing energy of ions; it is close to 50 eV for a wide range of ion’s energies close to the Bragg peak region (de Vera et al. 2013a, b). The transport of electrons with energies $\le 45$ eV can be described by diffusion equations, since their cross sections are nearly isotropic (Nikjoo et al. 2006). This approach was used in a series of works (Surdutovich and Solov’yov 2012, 2014, 2015; Solov’yov et al. 2009, 2011). Particularly, in Surdutovich and Solov’yov (2015) the expression for the radial dose $D(\rho )$ was obtained. This analysis corresponded to ions in the Bragg peak region, where the production of more energetic $\delta$-electrons is suppressed by kinematic constraints (Rudd et al. 1992). However, the effect of $\delta$-electrons was later considered (de Vera et al. 2017b), which allowed to extend the analysis to more energetic ions out of the Bragg peak region.

Fig. 1

(Figure adapted from de Vera et al. (2017b), with kind permission of The European Physical Journal D)

Radial doses produced in liquid water by a 200-keV/nucleon and b 2-MeV/nucleon carbon ions. While thin lines represent the time evolution of radial dose produced by low energy electrons calculated using the diffusion equations, the thick lines depict radial doses at the end of the track-structure including the $\delta$-electrons contribution. Symbols correspond to Monte Carlo simulations results (Waligórski et al. 1986; Liamsuwan and Nikjoo 2013; Incerti 2014). The associated radial pressure profiles and OH radical densities produced by $\sim$ 0.1 ps are shown on the right axis (see the text for details)

Figure 1 shows the radial doses in liquid water, calculated by the previously described formalism, for (a) carbon ions in the Bragg peak region (energy 0.2 keV/nucleon) and (b) carbon ions with a higher energy of 2 MeV/nucleon, out of the Bragg peak region (de Vera et al. 2017b). Carbon is one of the most promising ions used in cancer therapy. Thin lines depict the temporal evolution of the radial dose deposited by low energy electrons, while thick lines are the total radial doses (including the $\delta$-electrons contribution) by the end of the track-structure, by $\sim$ 50 fs. As it can be seen, at the Bragg peak region, where the amount of $\delta$-electrons is small, most of the dose is accounted for by low energy electrons, and the tail due to the transport of more energetic electrons is almost negligible. On the contrary, for 2 MeV/nucleon carbon ions, there is a noticeable contribution of $\delta$-electrons to the dose at large radii, with a long tail of low dose. The results of the model calculations are compared to Monte Carlo simulation results, depicted by symbols (Waligórski et al. 1986; Liamsuwan and Nikjoo 2013; Incerti 2014), showing a good agreement.

As can be seen in Fig. 1, the pressure [proportional to the radial dose by virtue of Eq. (2)] rises by several orders of magnitude within a few nanometres distance (the GPa scale is shown on the right). These conditions can induce the hydrodynamic response of the liquid medium to give place to cylindrical shock waves, and it is reviewed in the next subsection. It should be noted that the radial dose/pressure profiles depicted in Fig. 1 correspond to the end of the electron transport process and before a shock wave has developed (see next section); this is why the comparison of radial doses with Monte Carlo simulations (that do not consider shock waves) is good at this time.

Apart from the radial dose and the radial profile of pressure built up around each ion’s path, the diffusion model allows us to obtain the distribution of chemically reactive species generated by the end of the electron track-structure (Surdutovich and Solov’yov 2015). Each ionised (H$_2$O$^+$) or electronically excited water molecule (H$_2$O$^*$) can fragment to yield OH radicals and other species:

$$\begin{aligned} \text{H}_2\text{O}^+\longrightarrow & {} \cdot \text{OH} + \text{H}^+ \text{, } \nonumber \\ \text{H}_2\text{O}^*\longrightarrow & {} \cdot \text{OH} + \cdot \text{H}. \end{aligned}$$

(3)

The calculated densities of OH radicals, one of the main species responsible for biodamage, produced by carbon ions in and out of the Bragg peak region (Surdutovich and Solov’yov 2015; de Vera et al. 2018), are also shown in Fig. 1 on one of the right axes, which follow the same functional form as the radial dose.

Hydrodynamic response of the liquid medium

The steep and intense radial doses deposited around each ion’s path, which are translated into large pressure gradients, produce a region which can be referred to as a “hot cylinder”. The energy initially deposited into the electronic excitations of water molecules is released to the nuclear system by electron–phonon coupling, building up the initial discontinuity in instantaneous pressure/temperature. The first attempt to study these effects was made within the framework of the thermal-spike model (Toulemonde et al. 2009). From this work, it was roughly estimated that the maximum in instantaneous pressure and temperature after ion traversal happens by $\sim 0.1$ ps, time at which the energy transfer from the electron to the nuclei system can be considered to be completed. This time is much shorter than the characteristic timescales needed for heat dissipation by molecular diffusion.

However, it should be kept in mind that the thermal-spike model considers the medium as static, and cannot account for its dynamical response. Indeed, the ratio of pressures between the inside and the outside of this hot cylinder is very large. According to hydrodynamics, the needed ratio for the production of a “strong explosion” is $P_\text{in}/P_\text{out} > (\gamma +1)/(\gamma -1)$ (Landau and Lifshitz 1987), where $\gamma =C_\text{P}/C_\text{V}$. For water molecules with all their degrees of freedom available $\gamma =1.222$ and then $(\gamma +1)/(\gamma -1)=10$. The ratio of pressures in a nanometric cylinder around the ion’s path is much larger than 10, so the formation of a cylindrical shock wave propagating away from the ion’s path is to be expected. This conclusion was obtained independently in Surdutovich and Solov’yov (2010) in the framework of the MSA, although curiously it was already pointed out much earlier in Gol’danskiǐ et al. (1975), and since then ignored in the context of ion radiation biodamage. This scenario has to be analysed by solving the hydrodynamic equations of the system.

The hydrodynamic study of strong explosions following sudden and concentrated energy deposition has been conducted since the 1940s, and solutions for the spherical case were obtained independently by Zel’dovich (1942), Sedov (1946), von Neumann (1947) and Taylor (1950), which are summarised by Zel’dovich and Raiser (1966) and Landau and Lifshitz (1987). Following these derivations, a solution for the cylindrical case, relevant to the present problem, was derived in Surdutovich and Solov’yov (2010).

Two important characteristics of the cylindrical shock waves shown by solving the hydrodynamic equations are the position of the wave front as a function of time:

$$\begin{aligned} R(t) = \beta \sqrt{t} \left[ \frac{\text {LET}}{\varrho }\right] ^{1/4} \text{, } \end{aligned}$$

(4)

and the pressure of the wave front for a given position:

$$\begin{aligned} P(R) = \frac{1}{2\left( \gamma +1\right) }\frac{\beta ^4 \text {LET}}{R^2} {\text{. }} \end{aligned}$$

(5)

In these equations, $\varrho = 1$ g/cm$^3$ is the density of unperturbed liquid water and $\beta$ is a parameter which value for liquid water is $\beta = 0.86$ (Surdutovich and Solov’yov 2010).

Fig. 2

(Figure adapted from de Vera et al. (2017b), with kind permission of The European Physical Journal D)

Wave front characteristics for the shock waves produced by 200-keV/u and 2-MeV/u carbon ions. a, b the wave front position as a function of time, and c, d pressure at the wave front position. Lines correspond to the hydrodynamic model, while symbols show molecular dynamics simulations. See the text for further details

Figure 2 shows by dashed lines the position of the wave front as a function of time (a, b) and the pressure of the wave front as a function of its radial position (c, d), as predicted by the hydrodynamic model when assuming that the LET is equal to the ion’s stopping power S (i.e. average energy lost per unit path length), for carbon ions in the Bragg peak region (0.2 MeV/nucleon) and out of it (2 MeV/nucleon). The radial doses of these ions are shown in Fig. 1. As can be seen, the front propagates fast at the beginning, with its speed progressively decreasing. Similarly, the pressure on the wave front is extremely high during the first instants of the shock wave development, but then rapidly decreases as the wave propagates.

It should be noted that the hydrodynamic model has several limitations. First, it considers that all the energy lost by the ion is deposited in a narrow cylinder around the ion’s path, so the specific characteristics of the radial dose, or the losses of energy as a result of water radiolysis, are not taken into account. Second, it is a model developed for macroscopic hydrodynamics, and it is not obvious whether it is applicable to this problem on the nanometre scale. Third, although the produced forces are strong enough to cause chemical bond dissociation, they only act for a short period of time. Third, even though high pressures are predicted, they cannot be directly used to assess the possible DNA damage. All these effects can be taken into account in molecular dynamics simulations. These will be analysed in the next section, where the meaning of the rest of lines and symbols in Fig. 2 will be explained. Still, the hydrodynamic model will be very useful for testing and benchmarking the molecular dynamics results, and to envisage some other potential effects of shock waves.

Apart from the possible direct thermo-mechanical damage of DNA molecules, the shock waves also prompt a collective mass transport of water molecules radially moving away from the ion’s path, the mass flux of which is proportional to $\sqrt{\text {LET}}$ (Surdutovich et al. 2013b). Since this mass transport follows the wave front, we can estimate the time needed by it to travel a radial distance $\rho =R$, which is, from Eq. (4), $t_\text{front} = \frac{\rho ^2}{\beta ^2}\sqrt{\frac{\varrho }{\text {LET}}}$. This can be compared with the time needed for a chemical species to reach the same distance by diffusion, $t_\text{diffusion} = \rho ^2/D_\text{r}$, with $D_\text{r}$ being the diffusion coefficient of a given species. The ratio of both quantities is given by Surdutovich et al. (2013a, b):

$$\begin{aligned} \frac{t_\text{front}}{t_\text{diffusion}} = \sqrt{\frac{\varrho }{\text {LET}}}\frac{D_\text{r}}{\beta ^2} \text{. } \end{aligned}$$

(6)

For the relevant chemical species present in the medium, the diffusion coefficients are $<10^{-4}$ cm²/s (LaVerne 1989), which makes the ratio of Eq. (6) inferior to $10^{-3}/\sqrt{\text {LET (keV/nm)}}$. This quantity is much lower than 1, even for low-LET ions such as protons, which indicates that the shock wave is a means of transport of chemically reactive species much more efficient than diffusion. This aspect will be reviewed in more detail by means of molecular dynamics also in the next section.

Classical molecular dynamics simulation of the biological effects produced by ion-induced shock waves

Through the last sections, we have seen how the evolution of the ion’s track-structure puts the initial conditions for the development of shock waves, the main characteristics of which have been analysed by means of hydrodynamics. However, the hydrodynamic model is not enough for evaluating in detail all the possible biological consequences of the ion-induced shock waves. For this purpose, classical molecular dynamics simulations can give much more insights, providing details on the atomic and molecular scale.

The classical molecular dynamics technique (Allen and Tildesley 1989; Solov’yov et al. 2017) consists on following the classical trajectories of all the atoms of a given system, which can give detailed information on the system’s temperature (related to atoms’ velocities), pressure, potential energy stored in chemical bonds and, even, their rupture. It is based on the Born–Oppenheimer approximation, which allows decoupling the nuclear and electronic motion, owing to the their large difference in mass. Atomic nuclei can then be treated classically, moving in the field created by the electronic system. Usually, the landscapes of potential energy between groups of atoms, as a function of their mutual distances and angles, can be parameterised in the form of convenient analytical expressions or numerical data tables, known as force fields.

In the case of water and biological molecules, their potential energies result from their geometrical configuration, owing to the symmetry of the different chemical bonds arising from molecular orbital hybridization. For these molecules, the popular CHARMM force field (Chemistry at Harvard Macromolecular Mechanics) (MacKerell et al. 1998) has been developed in the last decades, allowing an straightforward simulation of biological systems. The particles which interact through chemical bonds or Coulomb and van der Waals forces are determined by the topology of the system, which defines the connectivity of the atoms. It should be noted that, in the standard CHARMM force field, this topology is fixed and cannot vary during the simulation (i.e. chemical reactivity is not allowed). However, this limitation can be overcome by means of reactive force fields, as it will be discussed later on.

Simulation of shock waves and correspondence to hydrodynamics

For the purpose of studying ion-induced shock waves in biological media, simulations are mainly done in liquid water boxes. Biomolecules, such as DNA strands (de Vera et al. 2016) or nucleosomes (Surdutovich et al. 2013), can be also introduced, all of these molecules being described by the CHARMM force field, and with their structures being available at the Protein Data Bank (Berman et al. 2000). Previous to any production simulation of shock waves, these systems must be first equilibrated at body temperature ($T=310$ K). From this equilibration, each atom i of the system will get an equilibrium velocity $v_i^\text{eq}$, which will follow the Maxwell–Boltzmann distribution. Explanations on how to set up these simulations can be found elsewhere (de Vera et al. 2016; Surdutovich et al. 2013a). Simulation results shown here were mainly obtained by means of the MBN Explorer (Meso-Bio-Nano Explorer) software (Solov’yov et al. 2012, 2017).

Due to the development of the ion’s track-structure, a distribution of deposited energy arises, as a result of which atoms will acquire non-equilibrium velocities $v_{i}^\text{neq}=\alpha \, v_{i}^\text{eq}$, once the energy transfer from electrons to nuclei by electron–phonon interaction is completed. This will naturally originate a shock wave in the simulation box if the energy deposited is large and concentrated enough. The velocity scaling parameter $\alpha$ is obtained through the expression for the kinetic energy (de Vera et al. 2016, 2017b, 2018; Surdutovich et al. 2013a):

$$\begin{aligned} \sum _i^{N_j} \frac{1}{2}m_{i,j}(v_{i,j}^\text{neq})^2 = \sum _i^{N_j} \frac{1}{2}m_{i,j}(\alpha _j \, v_{i,j}^\text{eq})^2 = \frac{3N_j k_\text{B}T}{2}+ f(\rho _j) S L \quad \forall \, j \text{, } \end{aligned}$$

(7)

where sums go over all the $N_j$ atoms (of mass $m_{i,j}$) inside concentric cylindrical shells j of radius $\rho _j$ and width $\Delta \rho$ around the ion’s path. The first term at the right-hand side of this equation is the kinetic energy of the atoms inside the cylindrical shell j before energy deposition ($k_\text{B}$ is the Boltzmann’s constant), while the second term is the amount of energy deposited in this shell by the ion. S is the ion’s stopping power, L the distance travelled inside the simulation box and the fraction of deposited energy $f(\rho _j)$ accounts for the details of the radial dose distribution (de Vera et al. 2017b).

The application of the previous scheme was improved in successive works. Initially, the so-called “hot cylinder” approximation was used, which more closely resembles the assumptions of the hydrodynamic model described above (de Vera et al. 2016; Surdutovich et al. 2013a) (as well as the conditions at the Bragg peak region). In this case, all the energy was deposited within a single 1-nm-radius cylinder around the ion’s path, so $f(\rho _j)=1$ for $j=1$ and $f(\rho _j)=0$ otherwise. Later on, the details of the radial dose were considered, which enter in this equation through the radius-dependent fraction $f(\rho _j)$ (de Vera et al. 2017b, 2018). This allowed to simulate shock waves not only in the Bragg peak region, but also for more energetic ions out of the Bragg peak region, where the transport of energy far away from the ion’s path by $\delta$-electrons is more relevant, as illustrated in Fig.1.

It should be noted that Eq. (7) assumes that all the energy initially deposited by the ion in the electronic system of the target is transferred to the nuclear motion degrees of freedom by electron–phonon coupling. Actually, some fraction of energy may be lost, due to the energy needed to produce water radiolysis products. However, most of this energy will be again returned to the medium due to the kinetic energy gained by water molecule fragments as well as the solvation of any charged products. Thus, most of the energy deposited will end up in the motion of the target atoms, and the general physical situation will be the same. In any case, to take this situation into account, some simulations were done assuming that only a fraction of the energy is transferred to nuclear motion. In Surdutovich et al. (2013a), 75% of the energy was invested in nuclear motion as a conservative estimate. The effects of this reduction in the deposited energy on the simulation outcomes will be discussed in the next subsection.

Fig. 3

(Adapted by permission from Surdutovich and Solov'yov (2017b))

Time evolution of the pressure wave produced by a carbon ion in liquid water in the Bragg peak region, as obtained from a molecular dynamics simulation

Before we start discussing biological effects of shock waves, let us analyse their basic characteristics as obtained from the molecular dynamics simulations, namely the radial pressure distribution and its time evolution (de Vera et al. 2016). Figure 3 shows the radial distributions of pressure generated in liquid water by a carbon ion in the Bragg peak region (0.2 MeV/nucleon), for several times after ion traversal. For this case, the hot cylinder approximation was used, which resembles more closely the assumptions of the hydrodynamic model. As can be seen, the distribution starts in the form of a sharp peak, where the maximum is the wave front. However, as the wave front propagates radially, it also starts widening and weakening.

The positions and pressures of the wave front can be obtained from these distributions (de Vera et al. 2016), and they are represented by triangles in Fig. 2, both for 0.2 and 2 MeV/nucleon carbon ions. Since in this case the simulations were done using the hot cylinder approximation, the results compare very well with the predictions of the hydrodynamic model, shown by dashed lines. This good agreement serves as a benchmark for the molecular dynamics results, and indicates that the hydrodynamic approach is still valid at the nanometre scale.

If the radial doses are used to more accurately set up the initial conditions for the molecular dynamics simulations, the results slightly change. The results for the wave front produced by 0.2 and 2 MeV/nucleon carbon ions in water using the radial doses shown in Fig. 1 are represented in Fig. 2 by circles. As can be seen, for carbon ions in the Bragg peak, there is a small reduction in the shock wave strength, coming from the fact that the initial pressure distribution is smeared out with respect to the hot cylinder approximation. On the contrary, for the higher energy 2 MeV/nucleon carbon ion, where the radial dose is spread out over much longer distances due to the larger production of $\delta$-electrons, there is an appreciable reduction in the shock wave strength. Still, the hydrodynamic model reproduces well the time evolution of the front characteristics if, instead of assuming that LET$\simeq S$ in Eqs. (4) and (5), an effective LET is used. It has been shown that this LET corresponds to the amount of energy deposited within the first 1–2 nm around the ion’s path (de Vera et al. 2017b). When this effective LET is used, the hydrodynamic model yields the results shown by solid lines in Fig. 2.

Direct DNA thermo-mechanical damage by shock waves

The first potential biological consequence of ion-induced shock waves is the direct induction of single-strand breaks in DNA molecules as a consequence of the large instantaneous pressures and temperatures generated in the medium, which we will refer to as thermo-mechanical damage. These effects have been studied both in free DNA strands (de Vera et al. 2016; Bottländer et al. 2015), and more systematically in histone-wrapped DNA forming a nucleosome (Surdutovich et al. 2013a).

Fig. 4

(Reprinted by permission from de Vera et al. (2017a))

Snapshots of the structure of a DNA strand in liquid water, impacted by shock waves induced by a carbon and an iron ion in the Bragg peak region. a–c times 0, 5 and 10 ps, respectively, for carbon in water, and d–f times 0, 5 and 10 ps, respectively, for iron in water (de Vera et al. 2017a). Only a small fraction of the real simulation box is shown here for clarity, see the text

As a first illustration of these effects, Fig. 4 shows the evolution of the shock waves initiated by a carbon ion and an iron ion (LET eight times larger) in the Bragg peak region in liquid water, where a DNA strand has been placed so its axis is at 2 nm distance from the ion’s path (de Vera et al. 2016). The atoms explicitly shown in Fig. 4a, d depict the water molecules inside the hot cylinders around the carbon and iron ions’ paths (which trajectories are perpendicular to the figure plane) at time zero. Panels (b) and (c) for carbon, and (e) and (f) for iron, show the system at 5 and 10 ps after ion passage, respectively. In this figure, only a part of the whole simulation box is shown for clarity; the whole system needs to be much larger than shown here, in order to avoid interaction of the wave front with the boundaries of the simulation box, which would result in artifacts ( de Vera et al. 2016). As clearly seen, the shock wave produced by the iron ion is stronger and propagates faster, as the hydrodynamic model predicts (Eqs. (4) and (5)), since the LET is eight times larger in this case. The DNA strand located next to the ion’s track is distorted, especially in the case of iron. This happens as well in the wrapped DNA in nucleosome (Surdutovich et al. 2013a). However, these changes in the primary structure may be reversible, so it is more interesting to get insight into the potential irreversible damage, i.e. the bond breakage which might lead to single-strand breaks.

So far, the studies of possible bond breakage in DNA by ion-induced shock waves have been done using the CHARMM force field, which does not explicitly allow for bond breakage or formation. However, a study performed by means of first-principles molecular dynamics predicted similar effects as classical molecular dynamics (Fraile et al. 2019), and future works will benefit from the new reactive force fields implemented in various codes, such as in MBN Explorer (Sushko et al. 2016b) or ReaxFF ( van Duin et al. 2001) (see “Open questions, new tools and further research” section). Within the non-reactive simulations, the potential energy stored in a certain DNA bond as its length varies around its equilibrium distance can be monitored in time, and each local maximum stored as an energy deposition event (Surdutovich et al. 2013a; de Vera et al. 2016). The typical energy needed to dissociate a DNA backbone bond is around 3 to 6 eV (Range et al. 2004), although it can be much less due to the presence of reactive species in the medium (such as aqueous electrons, which can dissociatively attach to DNA molecules), even down to 0.3 eV (Smyth and Kohanoff 2012). When the potential energy of the bond overcomes this threshold, it can be considered to be broken.

Fig. 5

(Adapted by permission from de Vera et al. (2017a))

Distribution of energy deposition events in the DNA backbone bonds, for shock waves produced in the Bragg peak region by a carbon and an iron ion passing by at 1 nm from DNA. Filled symbols are the results for a short DNA duplex (de Vera et al. 2016, 2017a), while open symbols are results for nucleosome (Surdutovich et al. 2013a)

Figure 5 shows the distribution of energy deposition events lying between $U_\text{bond}$ and $U_\text{bond}+\text{d}U_\text{bond}$ in the backbone bonds of DNA which are located within 1 nm distance from the incident ion’s path. Results are shown for carbon and iron ions in the Bragg peak region, the cases previously shown in Fig. 4, both for free DNA strands (solid symbols) ( de Vera et al. 2016) and for nucleosome (open symbols) (Surdutovich et al. 2013a). The main trends are similar for both targets, with an increase in the number of high energy events with the increase in the LET of the ion.² This shows that the histones have a little role in protecting DNA from thermo-mechanical damage.

In Surdutovich et al. (2013a), a more systematic analysis of these trends as a function of the LET was done for the case of the nucleosome. A series of ions in the Bragg peak region were studied: carbon, neon, argon and iron, with LET ranging from 900 to 7195 eV/nm. The bond energy distributions for neon and argon lie in between of those shown for carbon and iron in Fig. 5. They are all characterised by a linear behaviour in the linear-log scale, indicating that they correspond to Boltzmann distributions with some characteristic temperature (Surdutovich et al. 2013a). From these Boltzmann distributions, it is possible to estimate the average number of bond breaks produced as a function of the LET (probability for above-threshold events) (Surdutovich et al. 2013a).

Fig. 6

(Adapted from: Surdutovich et al. (2013a))

Probability of SSB induction by shock waves produced by ions with different LET, assuming three different threshold energies for bond breakage $U_{\rm break}$. The results are compared with the probability of SSB induction by electrons and chemically reactive species (Surdutovich et al. 2013a)

The probability for a single-strand break (SSB) induction is represented in Fig. 6 as a function of the ion’s LET for three values of energy threshold for DNA backbone bond breakage $U_\text{break}$, namely 2, 2.5 and 3 eV, and compared to the number of SSB induced by secondary electrons and reactive species as obtained from the MSA (Surdutovich and Solov’yov 2014; Surdutovich et al. 2013a). For the conservative estimate of $U_\text{break} = 3$ eV, the shock wave becomes the leading mechanism for SSB induction for LET larger than 5 keV/nm, which corresponds to ions heavier than Ar in the Bragg peak region. Such heavy ions are not used for ion beam cancer therapy, but they can be present in other environments (e.g. as a result of nuclear fission in power plants or cosmic radiation), so the thermo-mechanical damage should be considered in these cases.

Impact of shock waves on the radiation chemistry and indirect DNA damage

As discussed earlier in the manuscript, DNA damage by ion irradiation is mainly accomplished by direct impact of the secondary electrons, as well as by the formation, propagation and reaction of chemically reactive species, i.e. free radicals and aqueous electrons. Indeed, the experimental evidence using radical scavenger’s points out to the fact that indirect DNA damage frequently accounts for a large fraction of DNA damage, between 40 and 90% (Hirayama et al. 2009). This is due to the rich radiation-induced chemistry, which strongly depends on LET. A comprehensive review on the experimental evidence available on the yields and reactions of chemically reactive species in ion tracks can be found in the current series (Baldacchino et al. 2019). The prevalence of indirect DNA damage by chemically reactive species is especially true for low-LET radiation (e.g. photons, electrons and protons), where it is largest, and it progressively reduces as the LET grows, although still being significant for high-LET ions (Hirayama et al. 2009).

This is explained by the fact that, the larger the LET is, the more concentrated the energy deposition around the ion’s path is, and thus the larger is the concentration of produced reactive species. For large concentrations, they have large recombination probabilities, so they will disappear before they can attack DNA molecules. To illustrate this point, the hydroxyl radical distributions produced by a carbon ion’s track in the Bragg peak region (a) is compared to that produced by a 500 keV proton (b), which has nearly thirty times less LET, in Fig. 7. Hydroxyl radicals’ atoms are explicitly shown as spheres and the ion trajectory as a straight line, while the liquid water medium is represented as a transparent box. The boundaries of the simulation box can be seen as walls. The periodicity that can be seen in Fig. 7a is simply due to the fact that this simulation box has been replicated for this figure, for clarity purposes; the simulation box for carbon ions was much smaller than that for protons in the original work (de Vera et al. 2018). See de Vera et al. (2018) for further details. Clearly, radicals will recombine with much more probability for the carbon’s track than for the proton’s one.

Fig. 7

(Figure adapted from de Vera et al. (2018), with kind permission of The European Physical Journal D)

OH radical distributions produced in liquid water around a a 200-keV/u carbon ion’s path and b a 500-keV proton path. The solid line represents the ion’s path, while spheres represent the OH radical atoms (red = oxygen, white = hydrogen), and the rest of water molecules in the environment are not displayed for clarity. See the text for details

It is interesting to note, however, that for LET values of $\sim$ 1000 eV/nm or larger, the indirect DNA damage estimated from radical scavenger studies can still contribute up to $\sim$ 40% (Hirayama et al. 2009), while the yield of OH radicals (one of the main species responsible for biodamage) estimated from radiochemistry experiments is almost zero at this LET (Taguchi and Kojima 2005). This point will be discussed later on in the manuscript, but it might be related to the formation of ion-induced shock waves.

Shock waves have been proposed as an effective means of reactive species transport by the initiated collective flow (Surdutovich and Solov’yov 2010; Surdutovich et al. 2013a), even for the lower LET of light ions. This was demonstrated by means of Eq. (6), where the ratio of characteristic times to reach a given distance from the ion’s path by collective flow and by diffusion was obtained. Thus, shock waves might play an instrumental role in the propagation of free radicals and indirect damage, especially for large LET ions for which radical recombination is otherwise extremely probable.

As for the thermo-mechanical effects, molecular dynamics simulations are capable of providing detailed insights into the effects of shock waves on the radiation chemistry at the atomic and molecular level. However, classical molecular dynamics is not often suitable for modelling chemistry, since electrons are not explicitly included, and because molecular topology (i.e. chemical bonding) is normally fixed when using empirical force fields such as CHARMM.

Even though, in the last years a number of so-called reactive force fields have been under development. Reactive force fields still are used within classical molecular dynamics, but they are parameterised in a way in which chemical bonds can be broken and formed under certain conditions. Well-known examples of such reactive force fields are ReaxFF (van Duin et al. 2001) or REBO (Brenner et al. 2002). Reactive force fields allow to simulate the chemistry in large systems, where performing ab initio simulations of the chemical reactions would be prohibitive.

Within the multiscale modelling package MBN Explorer (Solov’yov et al. 2012, 2017), which has been used to obtain many of the simulation results shown in this review, a reactive force field based on CHARMM has been recently developed (Sushko et al. 2016a). Such a force field is very useful in the sense that it allows the simulation of not only the dynamical processes involving biomolecules (such as shock waves in biological media), but also of their possible chemical reactivity, in the context of a force field that has been exclusively developed for biological systems and extensively tested.

In short, the reactive CHARMM force field implemented in the MBN Explorer software (Sushko et al. 2016a) uses Morse potentials to describe bonding interactions, so their strength gradually decreases as the separation between bound atoms increases. Additionally, a cutoff distance parameter can be defined to characterise the distance beyond which a given bond can be considered to be broken (or within which a new bond is considered to be formed). A number of chemical rules are defined in the simulation input, describing the admitted chemical transformations, and thus the topology of the system is allowed to evolve in the course of the simulation according to these rules.

One of the main chemically reactive species involved in the indirect damage of DNA molecules is the hydroxyl radical, OH$\cdot$ (Douki et al. 1998), formed as a consequence of water molecule splitting, Eq. (3). The hydroxyl radicals can damage DNA, but they can also annihilate, among other reactions (von Sonntag 1987), via a recombination forming hydrogen peroxide:

$$\begin{aligned} \cdot \text{OH} + \cdot \text{OH} \longrightarrow \text{H}_2\text{O}_2 \text{. } \end{aligned}$$

(8)

This reaction was studied in de Vera et al. (2018) by means of the reactive CHARMM force field implemented in the MBN Explorer code. For the sake of simplicity, this reaction was considered as a representative example of the induced radiation chemistry, so all other chemical reactions were disregarded. Pre-solvated electrons and hydrogen radicals were not included in these simulations.

As explained earlier in the manuscript, the radial distribution of OH radical density follows the radial dose profile, so the diffusion equations approach for the transport of low energy electrons allows us to obtain the initial conditions for the induced radiation chemistry. In Fig. 7, the OH densities produced by (a) a 200 keV/nucleon carbon ion and (b) a 500 keV proton are shown. These distributions were obtained by using the calculated densities, as shown in Fig. 1, as explained in de Vera et al. (2018).

The cases of these two ions serve as good examples to explore the dynamics of OH radicals after irradiation. For the low-LET case of 500 keV protons, shock waves are expected to be very weak. Thus, in this case, it was assumed that they are not present, and molecular dynamics simulations of the pure diffusion of hydroxyl radicals were done in order to compare and benchmark the results against the commonly used Monte Carlo track-structure codes. On the contrary, for carbon ions in the Bragg peak region, as it has been seen before, shock waves are particularly strong, so this case is a good case study of the effect of shock waves on the radiation chemistry scenario.

Fig. 8

(Figure adapted from de Vera et al. (2018), with kind permission of The European Physical Journal D)

Diffusion and reactivity of the OH radicals produced by 500-keV protons in liquid water, where a shock wave is not expected to occur due to the very low LET. a Mean square displacement of OH radicals, b number of OH and H$_2$O$_2$ molecules, as a function of time. MBN Explorer results (histograms with error bars) are compared to GEANT4-DNA simulations (dashed and dotted lines) (Karamitros et al. 2014). See the text for details

Figure 8 shows the simulation results for 500 keV protons, where it was assumed that shock waves are not present (de Vera et al. 2018). Panel (a) shows the mean square displacement (MSD) of OH radicals, in the pure diffusive situation. From this plot, the diffusion coefficient $D_\text{OH}$ was obtained by means of the Einstein relation, MSD$=6D_\text{OH}t$. A diffusion coefficient of 0.22 Å$^2$/ps was obtained from these simulations (de Vera et al. 2018), which corresponds relatively well with the results of other recent simulations giving 0.3 Å$^2$/ps (Pabis et al. 2011) and with the value used in the popular simulation packages PARTRAC (Kreipl et al. 2009) or GEANT4-DNA (Karamitros et al. 2011), 0.28 Å$^2$/ps. The numbers of OH and H$_2$O$_2$ molecules during simulations are shown by histograms with error bars in Fig. 8b, representing the average and standard deviations of three independent runs. These are compared with the results reported using the GEANT4-DNA package ( Karamitros et al. 2014): dashed lines depict the original results, in which many reactions among reacting species are possible, while dotted lines are scaled in a way in which only reaction (8) is included [see de Vera et al. (2018) for more details]. Even though it is computationally expensive to reach long times with molecular dynamics simulations, it is clear that its results (when ion-induced shock waves are not considered) agree very well with the evolution predicted from Monte Carlo simulations, where shock waves are not yet included. This validates the molecular dynamics model used in de Vera et al. (2018).

Similar characteristics were studied for the case of a carbon ion in the Bragg peak region, where the shock waves are substantial. The shock waves develop naturally in the molecular dynamics simulations, as a result of the energy deposited along the ion’s track. However, the computer simulations allow one to artificially “switch off” this process in order to distinguish its influence in the evolution of the radiation chemistry scenario. Of course, this switching off is something unnatural, but it was used in de Vera et al. (2018) to assess how the shock wave modifies the radiation chemistry as compared to the traditional picture of biophysical models, where shock waves have not yet been considered. This artificial switching off was achieved by setting initial atomic velocities to their equilibrium ones, instead of scaling them by means of Eq. (7).

Fig. 9

(Figure adapted from de Vera et al. (2018), with kind permission of The European Physical Journal D)

a Mean square displacement of the OH radicals produced around a 200-keV/u carbon ion’s path. b Time evolution of the number of OH radicals and produced H$_2$O$_2$ molecules. Results of simulations where the shock wave is allowed to develop (transport by the collective flow, solid lines) and where it is artificially “switched off” (transport by diffusion, dashed lines) are shown, in order to demonstrate the shock wave effects on the radiation chemistry

Figure 9a compares the transport of radicals in simulations where pure diffusion of hydroxyl radicals in the static medium is considered (by artificially switching off the shock wave, as explained above) and their transport by the collective flow induced by the shock wave. For the latter case, the radicals are transported at a rhythm almost 80 times faster than diffusion, clearly demonstrating the capacity of the collective flows initiated by shock waves to effectively transport reactive species in the medium.

The reactivity of the OH radicals in the presence of the shock wave is illustrated in Fig. 9b. It depicts the average number of OH and H$_2$O$_2$ molecules, together with the standard deviations obtained after three independent runs, in the two cases where the shock wave is artificially switched off and naturally allowed to develop, respectively. As clearly seen, the evolution of the number of molecules is different when the collective flow is present or absent. The transport of radicals by the collective flow does not only propagate the radicals much faster than diffusion, but also prevents their recombination, both by their spreading and by creating harsh conditions in which the formation of the O–O bond is suppressed (de Vera et al. 2018).

This effective transport of reactive species by collective flows is taken into account in the MSA, and it accounts for the fluence of species reaching DNA molecules at a given radial distance from the ion’s path for ions of sufficiently large LET. Unfortunately, a systematic study of the space and time distributions of reactive species by means of molecular dynamics has not yet been completed. Therefore, at present the transport of radicals by shock waves can only be estimated from considerations using the hydrodynamic model.

The distance to which reactive species will be transported by the collective flow depends on the time during which the hydrodynamic expansion of the shock wave is active. As the front expands radially, its pressure drops, and the inner medium of radius $\rho _\text{in}<\rho$ suffers a rarefaction. When the surface tension pressure $\sigma /\rho _\text{in}$ arising in the inner medium ($\sigma$ is the surface tension) overcomes the pressure of the front (given by Eq.(5)), the shock wave front will stop progressing. Equating the forces due to the pressure at the wave front with that at its inner surface, we get (Verkhovtsev et al. 2018)

$$\begin{aligned} \frac{1}{\gamma +1}\frac{\beta ^4}{2}\frac{\text {LET}}{\rho ^2}2\pi \rho L = \frac{\sigma }{\rho } 2\pi \rho L \text{, } \end{aligned}$$

(9)

where $\rho$ is considered to be the same on the left- and right-hand sides, the thickness of the wave front being much smaller than $\rho$, and where L is the length of an ion’s track segment. A linear dependence of the range of radical propagation R on LET can be derived from the above equation (Surdutovich and Solov’yov 2014, 2018; Verkhovtsev et al. 2018). At the end of the hydrodynamic process, when the pressure is again uniform, the reactive species are expected to be uniformly distributed inside the region of radius R. However, the precise value of R depends on $\sigma$, which is unknown in the conditions arising in the shock wave.

In practice, the number of DNA single-strand breaks produced by chemically reactive species at a given radial distance $\rho$ from the ion’s path is currently estimated in the MSA by Surdutovich and Solov’yov (2014, 2018), Verkhovtsev et al. (2018):

$$\begin{aligned} \mathcal{N}_\text{r}(\rho ) = \mathcal{N}_\text{r}\, \Theta \left( R(\text {LET})-\rho \right) \text{, } \end{aligned}$$

(10)

where $\Theta$ is the Heaviside function and $R(\text {LET})=10\,\text{nm} \times \text {LET}/\text {LET}_\text{CBP}$, with $\text {LET}_\text{CBP}$ being the LET of carbon ions in the Bragg peak, and where $R_\text{CBP} \simeq 10$ nm corresponds to a conservative estimate for the Bragg peak of carbon (Surdutovich and Solov’yov 2014). $\mathcal{N}_\text{r}$ has been estimated to be 0.08 in Verkhovtsev et al. (2016), based on experimental information (Dang et al. 2011). Still, more work is needed to derive a more precise relation for $\mathcal{N}_\text{r}(\rho )$ from molecular dynamics simulations.

Experimental evidence

As explained above, the shock waves induced by ions with sufficiently high LET propagating in liquid water were predicted in Surdutovich et al. (2011), and more detailed initial pressure distributions were used in de Vera et al. (2017b). In a number of works (Yakubovich et al. 2011, 2012a; b; de Vera et al. 2016, 2017b, 2018; Surdutovich et al. 2013a, b), the features of these waves were explored by means of molecular dynamics simulations. Works (Surdutovich and Solov’yov 2014; Verkhovtsev et al. 2016) demonstrated their importance in the scenario of radiation biodamage, showing that cell survival curves, calculated within the MSA including the shock wave scenario, agree with experiments. Yet, the shock waves have not been experimentally observed. In this section, we summarise the indirect evidence of their formation and effects and discuss feasible ways in which this phenomenon could be directly observed.

Indirect experimental observations

Some indirect observations appear to be related to the formation of the ion-induced shock waves in biological media.

First of all, acoustic waves coming from the Bragg peak region can be detected [see e.g. Baily (1992)]. These signals have even been suggested as a possible means to monitor the dose delivery (Hickling et al. 2018). These acoustic waves can be the remnants of shock waves induced by ion beams in the area of the Bragg peak, compatible with the region of high LET. Indeed, microshock waves have been identified among the possible sources of acoustic signals induced by radiation (Lyamshev 1992).

Apart from this, cavitation has been observed in liquid water samples by using pulsed X-ray radiation (Stan et al. 2016). This effect closely resembles the scenario of ion-induced shock waves in liquid water, in the sense that intense and concentrated patterns of energy deposition are produced.

Also, an experiment of water radiolysis for very high-LET Ni and U ions was reported by LaVerne and Schuler (1996). In this experiment, it was observed that a significant fraction of free radicals could escape the ion’s track. This observation is in qualitative agreement with the predicted effects of the ion-induced shock waves discussed above. It was hypothesised by the same author, indeed, that such an effect may originate from the broadening of the ion track (LaVerne 2000). However, the quantitative picture of this effect is absent and this may be a subject of further research.

However, perhaps one of the most promising arguments in favour of the shock wave scenario is the success of the MSA in predicting cell survival curves as a function of radiation dose in very good agreement with experimental data, for a wide range of conditions of LET, oxygenation and level of repair (Verkhovtsev et al. 2016, 2018). The effects of shock waves on the radiation chemistry, as explained in the previous section, are incorporated into the MSA, so a reliable theoretical prediction of probability of cell survival strongly corroborates the shock wave scenario.

Finally, a number of works such as McNamara et al. (2015) claim that the LET affects cell survival curves separately of the dose. Physically, this effect has to be explained and Surdutovich and Solov’yov (2017b) suggested that the dependence of the range of propagation of reactive species on LET can be responsible for that. This dependence can be due to the radial collective flow due to shock waves thus providing another indirect argument in favour of their existence.

Possible direct experimental verification of ion-induced shock waves effects

Apart from all the above-cited indirect evidence, the development of an experiment which could directly measure the appearance of shock waves and/or their effects in biological media would be extremely desirable. The following ideas have been considered.

The first possibility is related to the transport of chemically reactive species by the collective flows initiated by shock waves. As shown in the previous section, the collective flow leaves a characteristic signature in the early radiation chemistry, in the sense that the number of surviving reactive species in liquid water at a given time after ion passage will be larger in the case they are propagated by the collective flow as compared to diffusion. It is still unclear for how long this difference in the amount of surviving species in water would last, since molecular dynamics simulations are computationally expensive, and so far they have only been conducted for several tens of picoseconds (de Vera et al. 2018). The pulsed water radiolysis experiments, in which the amount of reacting species can be measured by chemical means as a function of time, are still limited to relatively long times (μs–ns), still challenging for molecular dynamics, and also to moderate LET values, for which these effects would be more difficult to distinguish. Still, some studies have been reported down to some hundreds of picoseconds (Jonah and Miller 1977). Fortunately, new techniques are becoming available, which might be able to monitor the evolution of the radiation chemistry on the picosecond scale. Specifically, pump-probe setups using laser-accelerated ion pulses have already been reported, capable of generating short ion pulses and monitoring their effects on the picosecond scale (Dromey et al. 2016; Senje et al. 2017). If the appropriate spectroscopic signal to monitor the evolution of reactive species (aqueous electrons, hydroxyl radicals...) is available, this technique could in principle be used to study water radiolysis at short times to compare with results obtained via molecular dynamics simulations.

Another experimental possibility arises from the irradiation of nanometric water droplets. In Surdutovich et al. (2017), different possible mechanisms of water droplet disintegration, when irradiated by ions, were analysed. A spherical water droplet of radius R will keep its stability as far as it is not subjected to forces overcoming its surface tension. An irradiated droplet can disintegrate for a number of reasons: it can explode due to the repulsion of electric charges in the droplet as a result of ionisation, it can evaporate as a result of energy deposition (or spontaneously after some time exposed to vacuum), or it can fragment into smaller droplets due to the pressure exerted by the shock wave front. In order to observe the effect of ion-induced shock waves, the droplet sizes should be in the right range so that their effect is dominant over other possible causes of multifragmentation. All these factors were carefully analysed in Surdutovich et al. (2017). It was concluded that, for carbon ions in the Bragg peak region, for water droplets having a radius $25 \le R \le$ 600–1400 nm (depending on the charge state of the impinging ion), the only mechanism capable of fragmenting it is the pressure at the front of the shock wave, while the rest of contributions would not be enough to break the droplet. Also, for these sizes, the characteristic time of evaporation in vacuum is orders of magnitude larger than the times in which the shock wave develops.

Still, other experimental possibilities can be imagined. Just to mention some of them, one can think that nanometric bioprobes, such as DNA origami structures (Bald and Keller 2014), could be used to detect the damage induced in them by shock waves. Also, it is expected that the channel temporally formed in the liquid medium by the wave expansion, as well as the cavitation effect following it, would temporally modify its optical properties. Such changes may be monitored by means of laser beams. Some similar experiments have been already performed, where the changes in optical properties of an inorganic liquid containing $^{235}$U isotopes fissioned by neutrons were measured, as a result of the formation of bubbles by the propagation of the fission fragments (Seregina et al. 2003). This case needs to be carefully analysed in the context of the shock wave scenario, and similar dedicated experiments may be planned.

Open questions, new tools and further research

Apart from the important need for experimental techniques to directly observe and characterise ion-induced shock waves in biological media, as discussed above, there are many research directions which need to be addressed, especially in order to answer some fundamental questions regarding the radiation biodamage mechanisms by ion beams.

In particular, let us lay out the following questions:

• Has the increased amount of radicals escaping high-LET ion’s tracks observed by LaVerne and Schuler (1996) been observed by other authors? Has some explanation to this effect been yet proposed? Such observation goes in line with the effects of shock waves on radiation chemistry explained in the previous sections.

• Why there is a substantial contribution of indirect DNA damage for very high LET ions? The data available on indirect effects point out to the fact that they still contribute up to $\sim 30$% of DNA damage for ions with LET larger than 2 keV/nm, despite the relatively small G-value for OH radicals at this LET (Hirayama et al. 2009). In connection to the previous question, may the shock wave be preventing radicals from recombination for times long enough to explain the reported indirect DNA damage?

• At which level of LET would the shock waves start to have the dominant role in ion-induced (direct and indirect) biodamage? Are they only relevant for high-LET ions such as fission fragments or high-atomic number and high energy ions present in cosmic radiation, or also for lower LET ions such as the protons commonly used in radiation therapy? Under which circumstances would the models which do not include shock waves in their mechanistic description of biodamage start to produce unacceptable errors in their predictions?

To answer these complex questions would need a vast amount of work, and goes beyond the scope of the present review. However, it is interesting to lay them out, since these questions seem to point out to a current lack of fundamental understanding, especially related to what is usually called the physico-chemical stage of irradiation, i.e. the processes occurring in the time frame from $\sim 100$ fs to $\sim 100$ ps after ion traversal. These processes include the decay of electronic excitations into the nuclear system, the production of hot molecular fragments leading to the nascent radiation chemistry and the hydrodynamic response of the liquid medium, and the further dynamics and reactivity of these species, among themselves and with the relevant biological molecules, e.g. DNA.

The physico-chemical stage is challenging both from the theoretical and experimental points of view. The processes mentioned above happen in a timescale which is difficult to monitor experimentally. Furthermore, they involve several space, energy and time scales, from the electron relaxation and dissociation of bonds in individual water and other biological molecules, to the collective movement of the medium and the transport and reaction of species with biological macromolecules. The modelling of this requires a multiscale approach, capable of simulating big and complex systems, but exploiting the information got from ab initio approaches and experimental information for the smaller molecules.

So far, sophisticated track-structure codes, such as GEANT4-DNA (Bernal et al. 2015) or PARTRAC (Kreipl et al.2009), have been developed, which successfully describe the physical stage (short times, < 0.1 ps) and chemical stage (long times, ns-μs scale) of irradiation in agreement with available experimental data. However, the still relatively scarce amount of experimental information on the intermediate stages, i.e. the physico-chemical stage (especially for very high-LET radiation), has made difficult the validation of these codes at these time frames. The connection of the initial physical events to the late chemistry is usually accomplished in track-structure codes by defining a series of parameters (e.g. branching ratios for molecular fragmentation, thermalisation distances of hot molecular fragments and low energy electrons, diffusion coefficients and reaction rates of reactive species...), determination of which ensures a reasonable transition from one to the other. However, this is exactly the time frame in which the dynamical response of the medium might play an important role, as explained in the previous sections. Thus, it is clear the importance of understanding the role that ion-induced shock waves might have on the determination of the parameters needed to model the physico-chemical stage.

Fortunately, nowadays new experimental and theoretical tools are being developed, which are allowing to start exploring and better understand the physico-chemical stage of radiation:

• Experimentally, the development of laser-accelerated ion sources is allowing to perform water irradiation experiments with femto and picosecond resolution (Dromey et al. 2016; Senje et al. 2017; Gauduel et al. 2010). Modern pulsed-beam time-of-flight techniques are also allowing to monitor in detail the fragmentation pathways of water molecules, as well as the kinetic energy distributions of the ejected fragments (de Barros et al. 2009; Ferreira et al. 2017a; b). Nanomolecular probes, such as DNA origami (Bald and Keller 2014), nowadays allow to explore damage mechanisms on the nanometre scale.

• From the theoretical point of view, classical molecular dynamics is incorporating the new tools of reactive force fields (Sushko et al. 2016a; van Duin et al. 2001; Abolfath et al. 2011, 2013) and irradiation-driven molecular dynamics (Sushko et al. 2016b). These methods allow the simulation of the dynamics and reactivity of large and complex systems [and even the consideration of electrons within classical simulations (Su and Goddard 2007; Islam et al. 2016)], taking advantage of the detailed information which can be obtained from ab initio approaches, either pure quantum calculations or first-principles molecular dynamics for the smaller systems (Fraile et al. 2019; Smyth and Kohanoff 2012; Gaigeot et al. 2010; Kohanoff et al. 2017).

All these new tools will allow a deeper understanding of the processes connecting the initial physical processes of energy deposition by ions to the induced dynamical and chemical effects, as well as their link to the final biological outcomes. All these will be useful for a more accurate determination of the parameters utilised in track-structure codes. In particular, they will serve to advance on the study of the biological effects produced by ion-induced shock waves at high-LET and their inclusion in track-structure simulations.

In this context, much work is to be done with regard to the direct physical and the indirect chemical effects of shock waves. With respect to the direct effects

• The new tool of reactive force field has already been started to be applied to the evaluation of direct damaging effects of shock waves in DNA. Bottländer et al. (2015) used ReaxFF to simulate direct proton impact in DNA strands. Curiously, even though these authors considered LET values much lower than those considered by Surdutovich et al. (2013a), as reviewed in previous sections, they found larger probabilities to produce single-strand breaks. These differences may be related to the use of the reactive force field, to the direct vs indirect ion-impact situations, or to the disregard of a detailed energy deposition profile around the ion’s path. This case deserves a more detailed analysis, with an in-depth study of DNA damage by shock waves using reactive force fields, such as ReaxFF (van Duin et al. 2001) or reactive CHARMM (Sushko et al. 2016a).

• Related to the previous point, so far the analysis of DNA damage by shock waves has been limited to specific impact parameters and orientations between the ion’s path and the DNA molecule. A more precise evaluation of these effects needs a comprehensive analysis, where the effects of direct vs indirect ion impact, impact parameter, orientation and LET-level are taken into account. Such studies, moreover, may be very useful to interpret potential experiments of irradiation of DNA probes, such as DNA origami (Bald and Keller 2014).

With respect to indirect effects

• There is an urgent need to build up an atomic and molecular level model for the early radiation chemistry. Initial attempts are being done by several groups in connection to low- (Abolfath et al. 2011) and high-LET radiation (de Vera et al. 2018; Abolfath et al. 2013), as well as for plasma medicine applications (Bogaerts et al. 2014, 2016).

• This requires that reactive force fields (Sushko et al. 2016a; van Duin et al. 2001; Abolfath et al. 2011, 2013) include all the relevant chemical species in water radiolysis, not only free radicals, but also aqueous electrons. Some force fields are starting to be developed, where electrons can be explicitly considered (Su and Goddard 2007; Islam et al. 2016).

• New tools where the radiation (quantum) effects can be coupled to the system’s (classical) dynamics, such as irradiation driven molecular dynamics (Sushko et al. 2016b), need to be further developed. For these purposes, it is useful to count on with first-principles molecular dynamics simulation techniques (Fraile et al. 2019), which can yield the molecular input data needed for larger scale simulations.

• By using these tools, a detailed analysis of the effect of shock waves on the radiation chemistry and indirect DNA damage, for a wide LET range, is needed. This information is essential for improving the MSA and other biophysical models and, particularly, to understand the transition from the low- to high-LET regimes.

Summary and conclusions

In this review, the research done so far on the production of cylindrical shock waves on the nanometre scale around energetic ions’ paths, as well as their biological effects, have been discussed and summarised. The formation of ion-induced shock waves in biological media was predicted in the context of the MultiScale Approach (MSA) to the radiation damage induced by ions (Surdutovich and Solov’yov 2014, 2017a). They are a consequence of the fast energy transfer from the ion to the electronic excitations of water molecules in the close vicinity to the path, and to the transfer of these large energies to the motion of the nuclei by means of electron–phonon coupling. The energy is thermalised in the form of hydrodynamical response of the liquid medium, which has a series of consequences on the direct physical and indirect chemical damage of biomolecules.

The formations of these cylindrical shock waves and their effects have been analysed by means of hydrodynamics (Surdutovich and Solov’yov 2010) and molecular dynamics simulations (Yakubovich et al. 2011, 2012a, b; de Vera et al. 2016, 2017b, 2018; Surdutovich et al. 2013a, b). It has been found that shock waves are capable of producing DNA single-strand breaks (SSB), overcoming the action of secondary electrons and chemically reactive species, when the ion’s LET is large enough, for ions heavier than Ar in the Bragg peak region. Even though this situation is not typically found in therapy, such heavy ions can be found under other circumstances, such as accidental irradiation by nuclear fission products or by cosmic rays in space.

Still, even lower LET ions can have a large impact on the induced radiation chemistry. It has been predicted that the collective flows induced by shock waves can propagate chemically reactive species much faster than diffusion (Surdutovich and Solov’yov 2010; de Vera et al. 2018; Surdutovich et al. 2013a). Indeed, for the large concentrations of radicals produced along a carbon ion’s track in the Bragg peak region, it is hardly possible that reactive chemical species can escape the track by simple diffusion without recombination (Surdutovich and Solov’yov 2015). Thus, the shock waves might explain why indirect DNA damage by reactive species is experimentally observed even for very large LET ions, where the probability of recombination is extremely large (Hirayama et al. 2009).

Still, ion-induced shock waves have yet to be directly observed experimentally, and more work needs to be done to fully understand and model the physico-chemical stage of radiation and, in particular, the biological consequences of ion-induced shock waves. On the basis of the knowledge gathered on the topic, new research lines have been suggested for the experimental detection of shock waves and a more detailed modelling of the relation between physics and chemistry in ion beam irradiation of biomaterials.