Atomic-scale insight into the interactions between hydroxyl radicals and DNA in solution using the ReaxFF reactive force field

To date, the treatment of cancer remains a big challenge. Although most cancer cells are killed after irradiation or by the use of chemicals, i.e. chemotherapy, a small concentration of tumorous cells may survive and prove more difficult to eliminate using conventional techniques [1, 2]. Additionally, it is well established that the above mentioned therapies lack selectivity towards cancerous cells, being harmful for both cancer and healthy cells. Moreover, many studies point to the risk of the growing resistance of cancer cells to these therapies and to the body's own defense system [2–6]. This resistance is a result of the very nature of tumorous cells: the capability to mutate due to the continuous stress, e.g. oxidation, on its own DNA, changing the biochemical pathways within the cell. For this purpose, new and more selective treatments need to be investigated, which could prove safer for the host body, treating the malignant cells without damaging the normal, healthy cells.

Plasma medicine is becoming a fast growing field of research, which attracts increasing interest for biomedical applications. Among these applications are sterilization, dental treatment, blood coagulation, wound and ulcer healing, inflammation as well as cancer treatment, which are widely investigated by applying cold atmospheric pressure plasmas (CAPPs) [7–16]. A wide range of studies have proven that plasmas work on a multilevel scale, where the synergy of the different plasma species leads to the envisaged results for each application. An example is found in the area of wound healing, where plasmas are able to stimulate living tissue by acting on the cellular biochemistry, while simultaneously invoke blood coagulation [13].

With respect to cancer treatment, many positive results are observed, both in vitro and in vivo, indicating that plasmas are able to attack a wide range of cancer lines without damaging the healthy cells [17, 18]. Among the investigated cancer lines are breast cancer [19], cervical cancer [20], lung cancer [21], gastric cancer [22], leukemia [23], pancreatic cancer [24], liver cancer [25], ovarian cancer [26], melanoma [27], neuroblastoma [28], glioblastoma and colorectal carcinoma [29, 30], etc. In these studies it was found that senescence and apoptosis were the leading, selective, anti-tumor mechanisms at low plasma doses and are the result of the synergetic activity of the produced reactive oxygen and nitrogen species in the plasma (ROS and RNS, respectively). Indeed, this synergy will affect multiple parts of the cell, both stimulating and inhibiting several biochemical pathways, which results in the inability of the cell to reproduce and which can eventually cause cell death [31]. Furthermore, plasma treatments have proven to be beneficial against drug resistance of these malignant cells, as the plasma species can interact with the specific transmembrane proteins, enhancing conventional anti-cancer treatments [32]. This way, plasma species act together with the chemotherapeutic drugs, in a synergistic manner, to overcome drug resistance and to achieve a higher apoptotic rate in tumors. In a recent study of Ishaq et al, this synergism was clearly demonstrated [33]. Here, an increase in apoptotic responses was observed in TRAIL-resistant HT29 colorectal cancer cells, without harming the healthy cells, using the combined treatment of CAPPs with TRAIL (tumor necrosis factor-related apoptosis-inducing ligand) [33].

One of the effects of CAPPs, seen in multiple independent studies, is the elevation of the ROS levels (e.g., OH, HO₂, H₂O₂ and O₂⁻) within cells [31, 34, 35]. Recent studies have shown that CAPPs are able to transport ROS up to 1.5 mm within the surface of targets and that ROS are observed to penetrate phospholipid bilayer vesicles, as a model system for the cell-membrane, up to 150 μm deep in the sample [36, 37]. Furthermore, in a study by Hong et al [37] no ruptures of the bilayer were observed, after the introduction of ROS in the vesicles, suggesting the formation of short-lived punctures allowing the ROS to transfer through this phospholipid bilayer, although the exact molecular mechanisms are still unknown. In a recent investigation performed within our group, the behavior of a phospholipid bilayer, at several rates of phospholipid peroxidation caused by ROS, was studied [38]. In this computational study, we observed the formation of aqueous pores due to the rearrangement of the oxidized phospholipids. This created a polar channel, enabling the ROS to enter the cell, at least in model systems with a reduced cholesterol content compared to the steady-state concentration of healthy cells (usually found in the 40–50% range). As it is stated that a wide range of cancer cell lines contain a significantly lower concentration of cholesterol in their cell-membrane, compared to their healthy counterparts, we expect further accumulation of ROS within the treated cancer cell lines [39]. Within the cell, reactive species interact or react with different parts of the cell which results in the further stimulation of the cellular production of ROS. The ROS then need to diffuse through the cell, in order to oxidize certain cellular parts, e.g., nuclear DNA, thus having a high chance of being scavenged by the multiple antioxidants present in the cell (glutathion, superoxide dismutase, catalase, cytoglobins, etc).

Tumor-selective treatment is partly caused by the difference in the redox balance of cancerous cells compared to normal cells. In addition to the level of antioxidants, studies have pointed towards the plasma assisted activation of the antioxidant systems in both normal and cancer cells either directly by the presence of reactive species [40] or as a result of the increase in oxidative stress. As cancerous cells already contain a higher steady-state concentration of ROS [41–44], the increase in concentration, as a result of plasma treatment, can overcome the cellular antioxidant system. As a consequence, tumorous cells will suffer from damage caused by the oxidation of proteins, lipids and DNA, which may eventually lead to apoptosis or even necrosis if the ROS concentration and the consequent oxidative damage from the plasma treatment, become too high [30]. The reason for the elevated steady-state ROS concentration of cancer cells is not yet fully understood. Possible sources are changes in mitochondrial proteomics due to DNA mutations, hypoxia and anoxia, reduction in the antioxidative responses and errors in the electron transfer chain in the mitochondria as a result of the Warburg effect [7, 45].

Among the endogenous ROS, hydroxyl radicals (OH) are found to be the most reactive towards DNA, resulting in DNA strand breaks (both single strand breaks (SSBs) and double strand breaks (DSBs)), oxidation of the nucleotides and H-abstractions [46–48]. OH radicals are the resulting products after the Fenton reaction of the more stable H₂O₂. The latter forms one of the dominant plasma-generated ROS in aqueous solution and is furthermore partially produced in the mitochondria as a by-product in the electron transfer chain [49], resulting in a high concentration of H₂O₂ in the cytoplasma of cancer cells. The oxidation of DNA and their consequences, however, still remain the subject of many studies as outlined in a recent review paper [50].

Computer simulations are ideally suited for gaining knowledge and understanding of phenomena in biomolecular systems, complementary to experimental observations. In particular, classical molecular dynamics (MD) simulations provide an atomic scale insight in the chemistry and dynamics of relatively large systems, at least compared to quantum mechanical computational methods. In the past, our group already used reactive MD (rMD) simulations to elucidate the reactivity of plasma-generated ROS towards multiple biochemical structures [51], including parts of the cell wall of Gram-positive and Gram-negative bacteria (for sterilization) [52–54], a model system for the stratum corneum (for skin treatment) [55, 56] and a liquid layer, serving as a simple model system for biofilms [57, 58].

In the present work, rMD simulations are performed in a similar fashion, using the Reax potential [59], to investigate the interaction of OH radicals, as well as H₂O₂ molecules, with DNA, in order to provide atomic-scale insight into the dynamics of one of the leading mechanisms of plasma-driven apoptosis in mammalian cells, which will prove useful for plasma cancer treatment and even for oncology more generally. Using the ReaxFF potential, relatively large systems can be investigated, in a dynamically and computationally inexpensive manner [59]. In section 2, the computational setup will be explained, together with the treatment of the simulated system before impact simulations were performed. The results obtained for this system will be discussed in section 3, and finally, in section 4, a conclusion will be drawn.

In MD, the trajectories of all atoms in a given system are calculated by integrating the laws of motion. The forces, acting on every particle, are derived from an adequate interatomic interaction potential (force field), which governs all atomic interactions. In this work, we use the ReaxFF potential [59] based on the bond order/bond length relationship introduced by Abell [60]. In this potential, the total potential energy of the system is calculated as the sum over multiple partial energy contributions, including bond energy, valence and torsion angles, over- and undercoordination corrections, conjugation terms and the non-bonded van der Waals and Coulomb interactions. The ReaxFF potential is able to describe not only covalent bonds but also ionic bonds and a wide range of intermediate interactions, enabling us to investigate relatively large systems (up to 10⁵–10⁶ atoms) over a time scale of nanoseconds in a reactive manner. A more detailed description of the potential and the energy terms can be found in [60, 61]. A wide variety of parametrizations has been developed and optimized for ReaxFF, for the accurate representation of the chemistry of covalent materials [61–63], polymers [64, 65], combustion [61, 66], catalysis [67, 68], metals [69, 70], metal oxides [71–73], biomolecules [52–57], etc. For the simulation of the DNA system in this work, we have used a modified version of the force field developed by Monti [74], optimized for C/H/O/N/P systems. This force field is known as a reactive force field, as it enables to simulate bond formation and dissociation, providing fundamental insight in the reactivity of the simulated molecular system. The charge distribution is calculated based on geometry and connectivity, using the electronegativity equalization method [75, 76]. Simulations were carried out using ReaxFF as implemented in LAMMPS [77].

For the present study, a DNA string composed of 12 base pairs (i.e., a DNA dodecamer) was considered (see figure 1(a)). The dodecamer was placed in a periodic rectangular cuboid with dimensions 33 Å × 33 Å × 48 Å. The remaining volume of the simulation box was filled with water molecules (approaching a density of 1 g ml⁻¹), providing a simple model for DNA, consisting of only the specified ROS, water and DNA. Hydroxyl groups are added to the two ends of both DNA strands, replacing the missing phosphate group while preserving the O-atoms of the end nucleotides connected to either C5' or C3' (see figure 1(b)).

Zoom In Zoom Out Reset image size

Reactions at either one of these four end points are thus not considered in the results, i.e. reactions at the resulting OH groups (e.g. H-abstractions) as well as the adjacent atoms (e.g. bond breaks and formations) are neglected. For integrating the equation of motion, a time step of 0.25 fs was used. All impact simulations, 15 independent runs in total, were performed for 500 ps at 300 K.

Prior to the impact simulations, the molecular systems were equilibrated at room temperature for 300 ps in the canonical ensemble (temperature and volume were kept constant: NVT) using a Nosé–Hoover thermostat with coupling constant of 25.0 fs. We used the following procedure for the equilibration: the system was first slowly heated from 0 K to room temperature for 100 ps, after which the system was stabilized at 300 K for 200 ps. This simulation time was verified to be sufficient to stabilize the molecular structures at the given temperature. Additionally, test simulations have been conducted to verify the accuracy of the force field with respect to nucleotide stacking interactions, which are known to contribute significantly to the system stability and structure. For the guanine–adenine stacking, we find a stable intrastrand base distance of 3.4 Å, corresponding with a stacking interaction of −9.31 kcal mol⁻¹, in good agreement with the literature [78–80].

First of all, a clear difference between the reactivity of H₂O₂ molecules and OH radicals towards DNA was observed. While the OH radicals react upon impact, the H₂O₂ molecules were found to be non-reactive towards DNA. Instead, and only in rare cases, similar reactions have been observed as described in the work of Yusupov [57], where clusters of H₂O₂ and water result in the formation of OH and HO₂. More information can be found in [57]. This observation confirms the high reactivity of OH radicals, as well as the long-lived nature of H₂O₂. Additionally, we have investigated the interaction of HO₂ radicals with the dodecamer and found that they reacted with DNA in a similar fashion as the OH radicals, albeit at a much lower rate (14% of HO₂ resulted in an additional reaction, similar as depicted in figure 2, compared to 80% in the case of OH radicals). Therefore, in this work, we will only focus on the interactions of the OH radicals with DNA.

Zoom In Zoom Out Reset image size

3.1. OH-addition reactions

3.2. H-abstraction reactions

H-abstraction reactions at the nucleotides, caused by OH radicals, are also observed in our simulations, but to a lower extent, and they are presented in figures 3 and 4. Both reactions are a direct result from an H-abstraction at amines found on the nucleic acids, i.e., primary and secondary amines, resulting in 2-N-centered and 1-N-centered radicals, respectively. These reactions have only been encountered in 10% and 5% of the total number of observed reactions. The H-abstraction from a primary amine has been observed for both dAMP and dGMP, while the H-abstract from a secondary amine is only observed in the dGMP, found at the very end of the simulated DNA double helices, thus being more dissolved in water. This combined behavior has also been suggested by the works of Mundy and co-authors [90, 91] and Abolfath et al [92]. Using Car-Parrinello MD, Mundy et al investigated the H-abstraction reactions of guanine by OH radicals in gas phase [90]. Their observations pointed out that the H-abstraction reaction from 2-N is energetically significantly favored above 1-N, which is also partially reflected in our simulated reaction rates. In their later work by Wu and Mundy, they investigated the same H-abstraction reactions in aqueous solution [91], and they observed that the reaction rate for the H-abstraction at 1-N increases compared to their earlier work. On the other hand, Abolfath and co-workers concluded that the H-abstraction from the secondary amine is favored above the abstraction from the primary amine of guanine in aqueous solution, using QM/MM simulations [92]. These observations are in line with our results, as the abstraction from 1-N is only observed in the more water dissolved nucleotides (i.e., at the end points of the simulated DNA double helix), albeit a significantly lower rate (5% compared to 10%).

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

3.3. Partial unwinding of DNA at breaks

It is known from experimental studies that OH radicals can also lead to direct H-abstraction at the C5'-position of the phosphodiester bonds of the DNA backbones [93–95]. This will result in a ribose radical, which is able to react further, either with O₂ or another ROS, or by the direct desorption of a phosphate group [82], both giving rise to SSBs within the affected DNA-strand. A SSB is a reversible oxidation pathway in DNA, but combined with a second SSB at the opposite strand, in close vicinity, it results in irreversible DSBs, cleaving the DNA double helix, which in turn activates pro-apoptotic factors. Moreover, it has been observed that ribose radicals, after H-abstraction, are able to react with the opposite strand, again resulting in DSBs [93]. Single and DSBs are widely considered among the main antitumor responses, leading to apoptosis [9, 93–95], which forms the main goal of plasma cancer treatment.

Although this particular H-abstraction reaction was not yet observed in our simulations, a strong affinity was found between OH radicals and C-5', suggesting that this reaction would occur on a longer time-scale, and it is important to mention that the subsequent SSBs result in an increased contact surface for water-interactions. This is important because our simulations point out that oxidation reactions almost exclusively occur at positions, which are in direct contact with the solution (DNA backbone and the outer side of the nucleic bases). As SSBs result in an increased contact surface with the solution, we investigated the effects of DNA strand cleavages at the dodecamer ends, mimicking the situation as it is observed at DSBs and SSBs. Indeed, at such points a partial opening of the double helix has been observed in our simulations, as displayed in figure 5. It is seen that the H-bonds between opposite nucleic bases at such positions are broken and replaced with H-bonds with the surrounding water molecules. We believe that due to the introduction of these migrated nucleotides in the solution, more oxidation reactions at these locations can be expected (e.g., formation of 1-N centered radicals, figure 4). This may prove very difficult to repair by the cell's repair proteins, which, again, may activate pro-apoptotic pathways in the cell, ultimately leading to controlled cell death. Thus, we expect that if enough time is given to the simulated structure, more and more H-bonds between opposite nucleic bases may be replaced with H-bonds with water molecules. In a similar fashion it is expected that two independent SSBs in the vicinity of each other (i.e., within several tens of base pairs [93]) may eventually lead to irreversible DSBs.

Zoom In Zoom Out Reset image size

As a final remark, it needs to be mentioned that a number of oxidation reactions described in literature have not been encountered in our simulations, like the addition of OH radicals on sp² carbon atoms of pyrimidines (dTMP or dCMP) [96, 97]. Ji et al investigated this with thymine using ab initio DFT calculations and concluded that the addition of OH on 6-C and 5-C are the energetically most favored reactions between OH and pyrimidine nucleic acids, with yields of 60 and 30%, corresponding to experimental values [97]. These observations suggest a somewhat lower overall reactivity in our simulations compared to experimental and ab initio studies. However, it was not the aim of this work to describe all possible reactions, but rather to give dynamic insight in DNA oxidation based on (reactive) MD. The reactions observed in this work point towards the most relevant oxidation products, which can be expected when considering an aqueous solution containing only OH or H₂O₂ and a dodecamer.

We have investigated the initial interactions of OH radicals and H₂O₂ molecules with a DNA dodecamer in solution, by means of rMD simulations using the ReaxFF potential. This enables us to obtain atomic-scale insight in the dynamic nature of the interactions. From these calculations, several observations were made.

(i) H₂O₂ does not react directly with the biomolecule.

(ii) OH radicals react either through H-abstraction resulting in the formation of N-centered radicals, or by OH addition leading to the formation of 8-OH-adduct radicals. The latter oxidation products are known as the first steps towards the widely known and mutagenic 8-oxo-Gua, 8-oxo-Ade, FapyG and FapyA, which are known as markers for DNA oxidation.

(iii) A strong affinity between OH radicals and the C5' position of the DNA backbone is found, although this did not yet result in H-abstraction within the simulation time, leading to SSBs.

(iv) Partial opening of the DNA molecule as a function of time, at positions of SSBs, can result in increased oxidation at the more dissolved nucleotides, with might be very difficult to repair by the cell's repair proteins, thus activating pro-apoptotic pathways in the cell, ultimately leading to controlled cell death.

These observations point out that OH-radicals, and indirectly H₂O₂ molecules, are pro-apoptotic species. The presented results provide insight in the interaction dynamics of OH radicals, as one of the most aggressive ROS in the treated cells, as well as H₂O₂, with DNA. This is of particular importance for understanding the interaction of CAPPs with living tissue, in particular in the context of plasma cancer treatment. With respect to the use of CAPPs for cancer treatment, the discussed oxidations are expected to happen only in malignant cells. In this way, a selective anti-tumor treatment is provided, where damage is induced on cancerous cell lines, while leaving their healthy counterparts unharmed or even stimulating the defense mechanisms of the latter cell lines and their environment. However, more investigations are needed in the field of plasma-based oncology to elucidate the role of plasma species in triggering tumor specific apoptosis.

CCWV, ECN and AB acknowledge the contribution of J Van Beeck who is investigating the interaction between H₂O₂ and DNA using rMD simulations. Furthermore, they acknowledge financial support from the Fund for Scientific Research—Flanders (project number G012413N). The calculations were performed using the Turing HPC infrastructure at the CalcUA core facility of the Universiteit Antwerpen, a division of the Flemish Supercomputer Center VSC, funded by the Hercules Foundation, the Flemish Government (department EWI) and the Universiteit Antwerpen. TJ and DF gratefully acknowledge support from the European Research Council through the ERC-Starting Grant THEOFUN (Grant Agreement No. 259608).