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