Dependence of Electronic Properties of DNA ... - ACS Publications

Jun 20, 2019 - Modeling deoxyribonucleic acid (DNA) is a challenge since this task typically ... As an example, a three-layer QM/QM/MM model was used ...
0 downloads 0 Views 1MB Size
This is an open access article published under a Creative Commons Attribution (CC-BY) License, which permits unrestricted use, distribution and reproduction in any medium, provided the author and source are cited.

Article Cite This: ACS Omega 2019, 4, 10791−10798

http://pubs.acs.org/journal/acsodf

Dependence of Electronic Properties of DNA Single Strands on Size and Environment Andreas Mauracher* and Stefan E. Huber Institute of Ion Physics and Applied Physics, Leopold-Franzens-Universität Innsbruck, Technikerstraße 25/3, 6020 Innsbruck, Austria

Downloaded via 5.189.202.212 on August 29, 2019 at 22:30:55 (UTC). See https://pubs.acs.org/sharingguidelines for options on how to legitimately share published articles.

S Supporting Information *

ABSTRACT: Modeling DNA at the electronic structure level remains a computational challenge. However, for a detailed understanding of possible reaction pathways, consideration of the electron distribution is essential. Here, we report a density functional theory study on the length dependence of electronic properties of DNA single strands consisting of up to eight cytosine units. We also study systems in which the cytosine either at the 5′-end or 3′-end is replaced by the anticancer agent gemcitabine. These systems are investigated without and with Na counterions, both in the gas phase and in a polarizable continuum to model an aqueous solvent. We calculate the vertical electron affinities, vertical ionization energies, and solvation energies. We find a pronounced odd−even oscillation for the vertical electron affinity and vertical ionization energy for all systems without counterions in the gas phase as well as in the polarizable continuum. These oscillations completely vanish by adding the Na counterions. We discuss the implications of our results with respect to comparisons between experimental and numerical investigations. For studies involving, for example, oligonucleotides, comparison of size-dependent properties derived from theory and experiment might lead to a better understanding of experimental conditions, which are sometimes difficult to determine.

1. INTRODUCTION Modeling deoxyribonucleic acid (DNA) is a challenge since this task typically requires a multiscale approach due to the intrinsic properties of DNA such as length, flexibility, and polymorphic structure. For instance, whereas the distance between base pairs is in the range of a few Ångströms, chromatin fibers have a length of several millimeters. Moreover, various temporal scales are involved in biological processes of DNA, which can range from femtoseconds up to several days and longer. A further complication comes from the interaction of DNA with its surroundings in biological tissue. Theoretical modeling of DNA depends thus on the desired spatial and/or temporal resolution as well as on the employed methods. These methods range from electronic structure approaches over atomistic descriptions employing force fields to more coarse-grained schemes at the meso- and/or macroscale and have been employed to gain a detailed understanding of DNA structure and function. A recent review summarizing the progress on multiscale simulation of DNA can be found in ref 1. Molecular mechanics involving biomolecular force fields such as CHARMM,2,3 AMBER,4,5 GROMOS,6 and OPLS7 among others has been proven useful to obtain various properties of DNA sequences such as structural information.8 However, to obtain details about the electronic structure of a molecular system, quantum chemical methods need to be employed.1 Since these methods are computationally rather costly, they are usually limited to small model systems. For example, density functional theory (DFT) © 2019 American Chemical Society

was used to study B-like conformations of deoxydinucleoside monophosphates,9 electron interactions with a DNA duplex,10 G-quartet systems (i.e., building blocks of G-quadruplex structures) of DNA,11 or d(A) 3d(T)3 and d(G) 3d(C)3 minihelixes.12 An approach to combine different methods constructively is the ONIOM model in which different topographic layers are defined for which then different methodological levels are employed using, for example, force fields at an outer, larger layer modeling a solvent and quantum chemical methods modeling the solute. A comprehensive overview on the ONIOM approach is found in ref 13. As an example, a three-layer QM/QM/MM model was used to study the interaction of proflavine with two B-DNA decamers, that is, [deca(dG-dC)]2 and [deca(dA-dT)]2.14 One important aspect of investigations concerning DNA is interactions with low-energy electrons (LEEs). They are produced in cells in copious amounts as secondary species originating from the interaction of biological matter with highenergy primary radiation stemming from, for example, cosmic radiation, the presence of radioactive Radon (222Rn), or in the course of medical treatment.15 The energy distribution of these secondary species is broad, reaching values up to 100 eV with a local maximum around 10 eV.16 In the energy range below the Received: March 25, 2019 Accepted: June 7, 2019 Published: June 20, 2019 10791

DOI: 10.1021/acsomega.9b00819 ACS Omega 2019, 4, 10791−10798

ACS Omega

Article

Figure 1. Schematics of GemCn (left) and CnGem (right). Gemcitabine is marked in red.

maximum, electrons formed in close proximity to DNA can efficiently cause damage to the DNA by dissociative electron attachment (DEA), which may lead to single-strand breaks (SSBs).17 Indeed, upon interaction between biomolecules and low-energy electrons, the formation of a stable radical anion is considered to be one of the most probable and crucial steps in order to result in DNA damage.18 For this reason, (vertical) electron affinities have been explored amply for DNA, its constituents, and closely related biomolecules both by experiments and theory.18 On the other hand, electron impact ionization processes constitute dominant processes for electron molecule scattering phenomena. In this respect, the (vertical) ionization energies are an important parameter for the description of these processes, especially concerning the eventual differences between the gas phase and the condensed state; for example, see ref 19. Several mechanisms are available to repair resulting DNA damage; otherwise, they might lead to mutation, cancer, or apoptosis.20 On the other hand, LEEs can lead to a synergistic effect between radiation therapy and chemotherapy. Typical substances which are used in such concomitant radiochemotherapy are halogenated molecules or platinum-containing anticancer drugs.21 For a better understanding of the resonant processes, which may occur upon electron attachment to DNA, many experimental studies investigated DEA to building blocks of DNA (such as nucleobases, phosphate groups, or sugar constituents) in the gas phase; for example, see ref 15 and references therein. In recent years, efforts have increased to advance from gas phase measurements of isolated molecules to microsolvated systems in order to mimic the effects of the environment in a cell to some extent.22−25 The effects of microsolvation on electronic properties of small molecules have been investigated from a theoretical point of view already for many years.26,27 In the case of uracil, it was found that the addition of several water molecules leads to an increase of the adiabatic electron affinity of the system.28−30 Indeed, it is expected that the surroundings in the cell have a notable impact on the interaction of DNA with LEEs.31

One frequently encountered difficulty when experimental results are compared with theoretical investigations, especially in nucleic acid research, is the different size range of molecules under consideration. At the electronic structure level, usually, small model systems are considered, which then form the basis for predictions concerning larger systems studied in experiments. Such an example is a DFT study of 3′,5′-dXDP with X = T, C, G, and A, which were seen as the smallest subunits of a DNA single strand.32 In that study, vertical attachment energies as well as vertical electron affinities were calculated to approach reasonable predictions of electron attachment to DNA single strands. Another approach to infer relations between properties at the microscopic, that is, molecular, level and the macroscopic one are (quantitative) structure−activity relationship (QSAR) models, which often make use of regression models. In the framework of pharmacology, QSAR models relate certain descriptors, which are typically physicochemical parameters derived from theoretical investigations, to response variables, such as anticancer activity, which was investigated for various Pt(II) compounds.33 The success of any QSAR model depends on the selection of the proper descriptors. In general, it seems reasonable that the comparison of certain trends in computational studies and experimental investigations might lead to a better understanding of the underlying physical and chemical processes rather than the reproduction of single quantities obtained from experiment in theory and vice versa.34 Inspired by the experimental work of Rackwitz and Bald on the interaction of LEEs with oligonucleotides of different lengths placed on DNA origami platforms,35 we decided to investigate electronic properties of specific DNA single strands as a function of strand length by means of DFT. An unexpected and hence interesting outcome of the presentations and discussions at the Fifth International Conference “Dynamics of Systems on the Nanoscale” (DySoN 2018) held in Potsdam (Germany; for further details, see http:// dyson-conference.org) was that experimental conditions concerning the extent to which results are affected by solvation and counterions of deposited DNA single strands on surfaces 10792

DOI: 10.1021/acsomega.9b00819 ACS Omega 2019, 4, 10791−10798

ACS Omega

Article

Although atomic fluorine has a high electron affinity, replacement of the two hydrogen atoms at the 2′ carbon in C1 by fluorine atoms leads only to a slightly higher VEA in the case of Gem. Note that the adiabatic electron affinities derived at the PBE0/DZVP level of theory are both slightly positive and yield 0.05 eV for C1 and 0.43 eV for Gem. Inspection of the highest occupied molecular orbitals (HOMOs) for C1− and Gem− yields similar results; that is, the HOMOs are mostly localized in the region of the cytosine. The presence of the two fluorine atoms in Gem leads only to a slightly higher probability of the additional electron to be located at the sugar ring in the vicinity of the 2′ carbon. As a representative case, the HOMO of Gem− is shown in Figure 3 among other selected HOMOs.

are apparently less clear than might be assumed. By variation of DNA single-strand lengths, which can, in principle, be realized both experimentally and theoretically, a comparison of trends of electronic properties is envisaged to determine which conditions may fit the experimental situation better. In the present study, we aim to explore the interplay of these contributing factors to some extent. As a first step, we investigate the vertical electron affinities (VEAs) and vertical ionization energies (VIEs) of DNA single strands in the Bconformation consisting of up to eight deoxycytidine monophosphate molecules (dCMP), in short, Cn with n = 1,...,8. These molecules were altered by replacing either the deoxycytidine molecule (dC) with gemcitabine at the 5′- or 3′-end, denoted as GemCn−1, n = 1,...,8 or Cn−1Gem, n = 1,...,8, respectively. In Figure 1, a schematic of GemCn and CnGem is shown on the left-hand and right-hand sides, respectively. Gemcitabine is an anticancer substance in which the hydrogen atoms of the 2′ carbon of dC are replaced by fluorine atoms. Furthermore, we added as many Na counterions to these three systems as there are phosphate groups present. All systems were investigated in the gas phase, and employing a polarizable continuum model (PCM) in order to mimic an aqueous solvent.

Figure 3. Representative HOMOs of anionic molecules in the gas phase. For better visibility of the HOMOs, the 3′- to 5′-ends of the molecules are aligned from top to bottom. The colors light gray, dark gray, red, blue, orange, and purple indicate hydrogen, carbon, oxygen, nitrogen, phosphorus, and sodium atoms, respectively. Positive and negative isovalues of the shown HOMOs are depicted by red and green meshes.

2. RESULTS AND DISCUSSION The VEAs (based on the neutral geometries optimized with the AMBER99 force field) for Cn, Cn−1Gem, GemCn−1, CnNan−1, Cn−1GemNan−1, and GemCn−1Nan−1 for n = 1, ..., 8 calculated at the PBE0/DZVP level of theory in the gas phase and employing a PCM are shown in Figure 2 in the left- and

In the case of C2, C1Gem, and GemC1, we find positive values for the VEAs yielding 4.71, 4.77, and 4.62 eV, respectively. In all three molecules, the additional electron is located at the phosphate group. Connecting two C1 molecules to a C2 molecule via a PO2 group leads to an open-shell, doublet configuration, that is, a system exhibiting an unpaired electron. This is also the case for Cn with n = 4, 6, 8. By adding an additional electron, the afterward negatively charged molecule represents a closed-shell, singlet configuration, that is, a system in which all electrons are paired. The HOMO of C2− indicates that the additional electron is mainly located around the phosphate group and will most likely affect the 5′ and 3′ bonds in the case of geometrical relaxation of the molecule; see Figure 3. In general, we find a slight increase in the VEAs for C4, C3Gem, and GemC3 in comparison to C2, C1Gem, and GemC1 followed by an almost linear decrease of VEAs for C6, C5Gem, and GemC5 and our longest singlestrand oligonucleotides C8, C7Gem, and GemC7 reaching values of approximately 3.9 eV. The higher VEAs of neighboring single strands with an odd number of bases can be explained by the nature of the closed-shell, neutral system. An open-shell anion has a higher gain in energy due to a more flexible electronic relaxation in comparison to a closed-shell anion, which increases with increasing size of the molecule. In contrast to the dimers C2−, C1Gem−, and GemC1−, the HOMOs of the tetramers C4−, C3Gem−, and GemC3− are delocalized around the two exterior phosphate groups as well as around the nucleobase at the 3′-end. Also in these cases, the presence of the fluorine atoms hardly influences the shape of the HOMOs. For the HOMOs of the anions C5−, C4Gem−, and GemC4−, we find differences. For C5− (see Figure 3) and C4Gem−, two adjunct central bases yield the highest probability to find the additional electron, with some extension

Figure 2. Vertical electron affinities (VEAs) in eV for all molecules under consideration as a function of the strand length, n = 1,...,8. Left panel: VEAs in the gas phase; right panel (light blue background): VEAs employing a polarizable continuum model. The lines between the data points were added to guide the eye.

right-hand panels, respectively. For the discussion, we divide our results into four parts, single strands without and with counterions in the gas phase and single strands without and with counterions in the polarizable continuum. We begin with the results of the gas phase calculations. The VEAs for C1 and Gem are not depicted in Figure 2 since these values are both negative, that is, −0.85 eV for C1 and −0.80 eV for Gem. Note that these numbers obviously depend strongly on the size of the basis set employed in the calculations. Further details of the basis set dependence of the results can be found in the Methods and in the Supporting Information. 10793

DOI: 10.1021/acsomega.9b00819 ACS Omega 2019, 4, 10791−10798

ACS Omega

Article

The shift of the HOMOs from the phosphate groups toward the bases upon solvation is most pronounced in the case of the dimers, that is, C2−, C1Gem−, and GemC1− (see Figure 4). For C2−, the additional electron is localized at both bases as well as at the parts of the sugar ring at the 3′-end. Similar findings can be noted for C1Gem−, whereas here, the HOMO also partially involves the sugar ring at the 5′-end. In the case of GemC1−, the HOMO extends solely over the base at the 3′-end and the adjunct sugar ring; see Figure 4. The VEAs for C2, C1Gem, and GemC1 yield again very similar values of 6.50, 6.60, and 6.52 eV, respectively, and are approximately 1.8 to 1.9 eV larger than those in the gas phase. In the fourth case, that is, single strands with counterions in the polarizable continuum, we find similar VEAs to those in the gas phase ranging from 1.29 to 1.46 eV for C4Na3 and C8Na7, respectively. The similarity of the VEAs in the gas phase and aqueous solution seems surprising especially since the location of the additional electrons appears to be rather different. Whereas in the case of gas phase molecules, the additional electron was located at the Na counterions, here, the HOMOs are distributed over up to three nucleobases. In particular, it is found that the HOMOs form bridges between different layers of nucleobases. Considering experiments, this would point toward a selective chemistry depending on the conditions in which DNA is prepared. In the gas phase, an additional electron seems more likely to interact with the sugar−phosphate backbone, whereas in the aqueous solution, reaction pathways rather involving the DNA bases may be accessible. In Figure 5, we show the VIEs of all molecules under consideration calculated in the gas phase and employing a

of the HOMO to the connecting phosphate group. In contrast, in the case of GemC4−, the highest probability to find the additional electron is at the first and third phosphate groups (counting them beginning at the 5′-end). The HOMOs of C6− and C5Gem− show a similar pattern with a significant portion of the orbital appearing at two phosphate groups separated by one phosphate group in between. For GemC5−, the HOMO is localized simply at the phosphate group next to the 5′ end. In the case of the largest considered single strands, that is, for n > 6, the HOMOs become more delocalized occupying up to three phosphate groups as well as bases to a minor fraction. The next step in our analysis was to add Na counterions according to the number of phosphate groups present in the respective molecule, that is, CnNan−1, Cn−1GemNan−1, and GemCn−1Nan−1 with n = 2,...,8. These structures were optimized as well with the AMBER99 force field. In the left panel of Figure 2, we see that the VEAs are hardly affected by the length of the single strands and that the odd−even oscillation of the VEAs observed in the case of molecules without counterions has vanished. The VEAs range from 1.06 to 1.58 eV for GemC2Na2 and C8Na7, respectively. Inspection of the HOMOs reveals that the additional electron is mainly located at the position of the counterions. In the case of larger molecules, that is, n > 2, the highest probability to find the additional electron is at the Na atom next to the 5′-end, and the probability steadily decreases toward the 3′-end. As a representative case, we show the HOMO of C5Na4− in Figure 3. We now turn to the employment of a PCM in the case of single strands without counterions. The respective VEAs are shown in Figure 2 (right-hand panel). We observe a similar odd−even oscillation as found in the gas phase; however, there is a significant increase in the VEAs, which has also been found for smaller systems.29,30 In addition, the oscillation magnitude is also smaller than that in the gas phase, that is, the difference in the VEAs of, for example, C8 to C7 is 4.60 eV in the gas phase and 3.21 eV in the polarizable continuum. In particular, the VEAs of C1 and Gem are positive and read 1.33 and 1.40 eV, respectively. Also here, the VEAs are hardly affected by substituting the respective hydrogen atoms in C1 with fluorine atoms in Gem. In general, the HOMOs appear to be similar to those in the gas phase; however, in most cases, there is a higher likelihood to find the additional electrons at the bases. A representative selection of HOMOs is presented in Figure 4 (i.e., for Gem−, GemC1−, and C5−).

Figure 5. Vertical ionization energies (VIEs) in eV for all molecules under consideration as a function of the strand length n = 1,...,8. Left panel: VIEs in the gas phase; right panel (light blue background): VIEs employing a polarizable continuum model. The lines between the data points were added to guide the eye.

PCM in the left-hand and right-hand panels, respectively. As in the case of the VEAs, a pronounced odd−even oscillation is found for single strands without counterions both in the gas phase as well as in the aqueous solution. However, in contrast to the anionic systems, the single strands with an even number of bases have higher VIEs in comparison to their next neighbors with an odd number of bases. This is due to the fact that single strands with even-numbered bases are open-shell systems. Single strands with an even number of bases show a

Figure 4. Representative HOMOs of anionic molecules embedded in a polarizable continuum in order to model an aqueous solvent. For better visibility of the HOMOs, the 3′- to 5′-ends of the molecules are aligned from top to bottom. The colors light gray, dark gray, red, blue, orange, and purple indicate hydrogen, carbon, oxygen, nitrogen, phosphorus, and sodium atoms, respectively. Positive and negative isovalues of the shown HOMOs are depicted by red and green meshes. 10794

DOI: 10.1021/acsomega.9b00819 ACS Omega 2019, 4, 10791−10798

ACS Omega

Article

rather linear increase in the VIEs, whereas those with an odd number of bases show an almost linear decrease of VIEs. The VIEs of C1 and Gem in the gas phase read 8.12 and 8.27 eV, respectively. The lowest VIEs are found for C7, C6Gem, and GemC6, which yield 5.22, 5.10, and 5.00 eV, respectively. As mentioned in the Introduction, at electron energies below 10 eV, electron attachment involves resonant formation of anions. However, the size dependence of the VIE implies that electron impact ionization remains to be a possible channel accessible at those low electron energies. On the other hand, for the dimers C2, C1Gem, and GemC1, the VIEs are approximately 0.7 eV higher than for the monomers. For the largest open-shell, neutral systems C8, C7Gem, and GemC7, the VIEs rise up to 11.45, 11.63, and 11.63 eV, respectively. At such energies, electron attachment is rather unlikely to occur. Adding Na counterions to single strands in the gas phase also affects the VIEs, yet not as pronounced as in the case of VEAs. The VIEs of the closed-shell, neutral systems, that is, single strands with an odd number of bases, are affected only slightly. The largest differences here are observed for C7Na6, C6GemNa6, and GemC6Na6 where the VIEs increase by 0.87, 0.65, and 0.96 eV, respectively. For single strands with an even number of bases, the VIEs decrease substantially with the largest differences observed for C8Na7, C7GemNa7, and GemC7Na7 of 2.90, 3.12, and 3.17 eV, respectively. A difference to the gas phase concerning the VIEs is also observed when single strands without Na counterions are considered in a polarizable continuum. In comparison to the gas phase, the VIEs of closed-shell, neutral molecules decrease, whereas the VIEs of open-shell, neutral molecules increase. Only in the case of the pentamers, that is, C5, C4Gem, and GemC4, the differences are below 0.13 eV. The VIEs of openshell, neutral single strands are below 9.55 eV and thus may also be candidates for electron attachment in contrast to their counterparts in the gas phase. Also in the polarizable continuum, the strand lengthdependent odd−even oscillations of the VIEs vanish when counterions are added to the single strands. The VIEs range from 6.37 to 6.62 eV, which is a narrower range than that in the gas phase where the VIEs take values in the interval from 5.86 to 7.27 eV. Interestingly, the VIEs of the trimers C3, C2Gem, and GemC2 with or without counterions and in the gas phase or in a polarizable continuum differ at most by only 0.64 eV. We note that the VEAs of single strands without counterions in the gas phase are generally lower than the corresponding VIEs. In contrast, in aqueous solution, the values of VEAs are close to those of VIEs. However, it is questionable if the probability to channel the energy released in a molecule upon electron attachment into a certain reaction mode, which then leads to the release of an electron elsewhere in the molecule, is substantial enough to observe this process in experiment. Adding counterions to single strands in the gas phase as well as in a polarizable continuum leads to rather similar VEAs and VIEs, whereas the VEAs are significantly lower than the VIEs, and thus a process as described above will certainly not occur. In Figure 6, we show the solvation energies Esolv as derived by eq 3 for all systems under consideration. In most cases, Esolv increases with the length of the molecule regardless of its charge state. We start with the discussion of single strands without counterions. Esolv is the largest for most of the anionic systems (see Figure 6, central panel with light blue background) followed by cationic systems (see Figure 6, right panel

Figure 6. Solvation energies Esolv in eV for the neutral, anionic, and cationic systems under consideration are shown in the left, middle (light blue background), and right (light red background) panels, respectively. For a better comparison, the scales of Esolv are the same for all three panels. The lines between the data points were added to guide the eye.

with light red background) and neutral systems (see Figure 6, left panel). The sole exceptions are C6+, C8+, C7Gem+, and GemC7+, and their solvation energies are 0.24, 0.22, 0.30, and 0.57 eV higher, respectively, than those in the case of anions. In the case of molecules without counterions, Esolv shows small deviations from a linear increase with the strand length for all charge states. It appears most notable to us that for cationic systems, Esolv is smaller for the heptamers C7+, C6Gem+, and GemC6+ in comparison to the hexamers C6+, C5Gem+, and GemC5+ by 0.20, 0.28, and 0.16 eV, respectively. For single strands with counterions, Esolv rises almost linearly with the strand length for all charge states but with an apparently steeper slope than in the case without counterions.

3. CONCLUSIONS We studied electronic properties of DNA single strands in Blike conformations consisting of one up to eight cytosine bases, that is, Cn with n = 1,...,8, at the PBE0/DZVP level of theory, both in the gas phase and employing a polarizable continuum model to mimic an aqueous solvent. The input geometries were optimized with the AMBER99 force field. These structures were altered by replacing a cytosine either at the 5′- or 3′-end with gemcitabine, that is, GemCn−1 or Cn−1Gem with n = 1,...,8. For these systems, we calculated the vertical electron affinities (VEAs), vertical ionization energies (VIEs), and solvation energies of all three charge states. For both, VEAs and VIEs, we find pronounced odd−even oscillations, which can be explained by the nature of the underlying either closed-shell or open-shell neutral system. Open-shell, neutral single strands form closed-shell anions provided that they are in their lowest electronic configuration. This leads to lower electron affinities than those in the case of closed-shell neutral systems. Although electronic relaxation is allowed after the attachment of the electron, the singly occupied molecular orbital of the neutral molecule restricts the space of the additional electron. On the other hand, electron attachment to closed-shell, neutral systems do not have such a constraint, and the electronic relaxation leads to higher VEAs. In the case of VIEs, the situation is reversed. Here, closed-shell, neutral molecules are associated with lower VIEs than those of open-shell systems since upon ionization, the single electron remaining in the HOMO is more strongly bound than each of 10795

DOI: 10.1021/acsomega.9b00819 ACS Omega 2019, 4, 10791−10798

ACS Omega

Article

developed for proteins, nucleic acids, and organic molecules.4 We chose AMBER99 among other options (as outlined in the Introduction) for its known performance to reproduce structural properties of DNA in good agreement with experimental data.37,38 Indeed, it is our intention to investigate electronic properties of DNA single strands, which are close to experimental conditions. Using these input structures, vertical electron affinities (VEAs) and vertical ionization energies (VIEs) were then determined by calculating the differences of the total electronic energies, Eel, between the neutral (M) and singly charged molecules (M+/−), that is

the electrons in the pair occupying the HOMO of the neutral molecule. The solvation energies do not show this odd−even oscillation and increase typically with the strand length. The influence of some environmental conditions was investigated by introducing Na counterions to the systems according to the number of phosphate groups present in the molecules. In this case, both VEAs and VIEs depend hardly on strand length, and particularly, the odd−even oscillations vanish. Moreover, the influence of an aqueous solvent was investigated by means of a polarizable continuum model. In the case of bare single strands, we found similar results to those in the gas phase concerning the shape of odd−even oscillations. However, we found that the respective VIEs are considerably lower, such that dissociative ionization could become a competing process to dissociative electron attachment, which describes the resonant anion formation at electron energies below 10 eV. For single strands with counterions in a polarizable continuum, we found only slight changes in comparison to the gas phase. Although the identified trends for the considered strand lengths partially exhibit an intriguing dependence on strand length (like the steep increase and decrease of the VIEs for even- and odd-numbered Cn, Cn−1Gem, and GemCn−1, respectively, and other similar features, see Figures 5 and 6), we would advise to be cautious with straightforward extrapolation of these results to much larger systems. Such length dependence would be interesting to investigate experimentally, which might be more easily accessible than in theory. We think that testing the length dependence of certain observables in experiment and comparison to numerically obtained trends as presented in this work could lead to a better understanding of the experimental conditions in studies concerning electron interaction with DNA single strands.

VEA(M) = E el(M) − E el(M−)

(1)

and VIE(M) = E el(M+) − E el(M)

(2)

For calculating the total energies, we employed the PBE0 density functional39 in conjunction with the split-valence DZVP basis set augmented with polarization functions,40 which was specifically optimized for DFT calculations; see also ref 9. Note that the calculation of electronic properties using a higher level of theory than that used for the determination of structural properties is a typical approach in quantum chemistry followed by various composite schemes aiming specifically at very high accuracies.41−43 To estimate the influence of the input geometries, we reoptimized C1 and Gem at the PBE0/DZVP level of theory and found only slightly different VEAs and VIEs (within 0.10 eV). Furthermore, to estimate the dependence of VIEs and VEAs on the basis set size, we calculated VIEs and VEAs for Cn with n = 1, ..., 4 and for CnNan−1 for n = 2, ..., 4 employing the Pople-type basis sets 6-31+G(d) and 6-311G(d).44 We chose these two basis sets to take into account diffuse functions and valence triple-zeta functions. In comparison to DZVP, the VIEs and VEAs change the most by 0.06 and 0.09 eV (given that the anion is stable at all, which is not the case for C1), respectively, when the 6-31+G(d) basis set is employed. In the case of using the 6-311G(d) basis set instead of DZVP, the VIEs and VEAs change the most by 0.15 and 0.22 eV, respectively. The tabulated values are listed in the Supporting Information. Based on these results, we assume that the choice of the DZVP basis set is appropriate to investigate trends of electronic properties for the systems under consideration. We also compared the highest occupied molecular orbitals (HOMOs) for these seven molecular systems (by visual inspection). Here, we also find small differences for the stable anions for the three basis sets DZVP, 6-31+G(d), and 6-311G(d). We show the HOMOs of Cn− with n = 1,...,4 and CnNan−1− with n = 2,...,4 in Figures S1 and S2 in the Supporting Information. More importantly, the inevitable small quantitative differences do not affect the conclusions of our study. We also calculated VEAs and VIEs employing the PCM of Tomasi and co-workers45 with a dielectric constant of water (ε = 78.3553) to explore possible influences of an aqueous solvent. Whereas the use of an implicit solvent model bears the obvious advantage of being computationally inexpensive and easily available and applicable, the explicit molecular resolution of solvent constituents is obviously sacrificed. Explicit solvent models would become important when interactions beyond polarization effects on the solute are considered, such as (partial) charge transfer from the solute to its environment. On

4. METHODS We investigated DNA single strands in the B-conformation consisting of up to eight deoxycytidine monophosphate molecules (dCMP), in short, Cn with n = 1,...,8. These molecules were altered by replacing either the dC with gemcitabine at the 5′- or 3′-end, denoted as GemCn−1, n = 1,...,8 or Cn−1Gem, n = 1,...,8, respectively. Gemcitabine is a derivative of deoxycytidine in which the hydrogen atoms of the 2′ carbon are replaced by fluorine atoms. This leads to 23 different molecules in total. Cn with n = 2, 4, 6, 8 and GemCn−1 and Cn−1Gem with n = 2, 4, 6, 8 all have an odd number of electrons and exhibit a doublet configuration (one unpaired electron, open-shell system) in their respective electronic ground state. Cn with n = 1, 3, 5, 7 and Cn−1Gem and GemCn−1 with n = 1, 3, 5, 7 have an even number of electrons and exhibit a singlet configuration (only paired electrons, closedshell system) in their respective electronic ground state. For Cn, GemCn−1, and Cn−1Gem with n > 1, Na counterions were added according to the number of phosphate groups present in the respective molecule leading to 21 more molecules. All molecules with Na counterions have an even number of electrons and exhibit a singlet configuration in their respective electronic ground state. The input structures of the neutral molecules used in the quantum chemical calculations specified above were optimized by minimizing the potential energy employing the AMBER99 force field5 via the steepest-descent algorithm implemented in the HyperChem 8.0 package.36 This is a revised version of the fixed charge additive force field by Cornell et al., which was 10796

DOI: 10.1021/acsomega.9b00819 ACS Omega 2019, 4, 10791−10798

ACS Omega

Article

and Its Parameterization. In Encyclopedia of Computational Chemistry; Schleyer, P. V. R.; Allinger, N. L.; Clark, T.; Gasteiger, J.; Kollman, P. A.; Schaefer, III, H. F.; Schreiner, P. R., Eds. John Wiley & Sons: Chichester, 1998; pp 271−277. (4) Cornell, W. D.; Cieplak, P.; Bayly, C. I.; Gould, I. R.; Merz, K. M.; Ferguson, D. M.; Spellmeyer, D. C.; Fox, T.; Caldwell, J. W.; Kollman, P. A. A Second Generation Force Field for the Simulation of Proteins, Nucleic Acids, and Organic Molecules. J. Am. Chem. Soc. 1995, 117, 5179−5197. (5) Wang, J.; Cieplak, P.; Kollman, P. A. How well does a restrained electrostatic potential (RESP) model perform in calculating conformational energies of organic and biological molecules? J. Comput. Chem. 2000, 21, 1049−1074. (6) van Gusteren, W. F.; Berendsen, H. J. C. Groningen Molecular Simulation (GROMOS) Library Manual; BIOMOS: Groningen, The Netherlands, 1987; pp 1−221. (7) Jorgensen, W. L.; Tirado-Rives, J. The OPLS [optimized potentials for liquid simulations] potential functions for proteins, energy minimizations for crystals of cyclic peptides and crambin. J. Am. Chem. Soc. 1988, 110, 1657−1666. (8) MacKerell, A. D., Jr. Empirical force fields for biological macromolecules: Overview and issues. J. Comput. Chem. 2004, 25, 1584−1604. (9) Poltev, V. I.; Anisimov, V. M.; Danilov, V. I.; Deriabina, A.; Gonzalez, E.; Jurkiewicz, A.; Leś, A.; Polteva, N. DFT Study of B-like Conformations of Deoxydinucleoside Monophosphates Containing Gua and/or Cyt and their Complexes with Na+ Cation. J. Biomol. Struct. Dyn. 2008, 25, 563−571. (10) Gu, J.; Wang, J.; Leszczynski, J. Electron interaction with a DNA duplex: dCpdC:dGpdG. Phys. Chem. Chem. Phys. 2016, 18, 13657−13665. (11) Villani, G. Quantum mechanical investigation of G-quartet systems of DNA. New J. Chem. 2017, 41, 2574−2585. (12) Gorb, L.; Zubatiuk, T. A.; Zubatyuk, R.; Hovorun, D.; Leszczynski, J. d(A)3d(T)3 and d(G)3d(C)3 B-DNA mini-helixes: the DFT/M06-2x and DFT/B97-D3 comparison of geometrical and energetic characteristics. J. Mol. Model. 2017, 23, 289. (13) Chung, L. W.; Sameera, W. M. C.; Ramozzi, R.; Page, A. J.; Hatanaka, M.; Petrova, G. P.; Harris, T. V.; Li, X.; Ke, Z.; Liu, F.; Li, H.-B.; Ding, L.; Morokuma, K. The ONIOM Method and Its Applications. Chem. Rev. 2015, 115, 5678−5796. (14) Ruiz, R.; García, B.; Ruisi, G.; Silvestri, A.; Barone, G. Computational study of the interaction of proflavine with d(ATATATATAT)2 and d(GCGCGCGCGC)2. J. Mol. Struct.: THEOCHEM 2009, 915, 86−92. (15) Baccarelli, I.; Bald, I.; Gianturco, F. A.; Illenberger, E.; Kopyra, J. Electron-induced damage of DNA and its components: Experiments and theoretical models. Phys. Rep. 2011, 508, 1−44. (16) Pimblott, S. M.; LaVerne, J. A. Production of low-energy electrons by ionizing radiation. Radiat. Phys. Chem. 2007, 76, 1244− 1247. (17) Boudaïffa, B.; Cloutier, P.; Hunting, D.; Huels, M. A.; Sanche, L. Resonant Formation of DNA Strand Breaks by Low-Energy (3 to 20 eV) Electrons. Science 2000, 287, 1658−1660. (18) Gu, J.; Leszczynski, J.; Schaefer, H. F., III Interactions of Electrons with Bare and Hydrated Biomolecules: From Nucleic Acid Bases to DNA Segments. Chem. Rev. 2012, 112, 5603−5640. (19) Huber, S. E.; Mauracher, A. Electron impact ionisation cross section of fluoro-substituted nucleosides. European Physical Journal D 2019, in production. (20) Sancar, A.; Lindsey-Boltz, L. A.; Ü nsal-Kaçmaz, K.; Linn, S. Molecular Mechanisms of Mammalian DNA Repair and the DNA Damage Checkpoints. Annu. Rev. Biochem. 2004, 73, 39−85. (21) Rezaee, M.; Hill, R. P.; Jaffray, D. A. The Exploitation of LowEnergy Electrons in Cancer Treatment. Radiat. Res. 2017, 188, 123− 143. (22) Neustetter, M.; Aysina, J.; da Silva, F. F.; Denifl, S. The Effect of Solvation on Electron Attachment to Pure and Hydrated Pyrimidine Clusters. Angew. Chem., Int. Ed. 2015, 54, 9124−9126.

the other hand, methods that take into account solvent molecules explicitly become easily computationally unfeasible at the electronic structure level. Even for more coarse-grained approaches, difficulties for an appropriate treatment of (aqueous) solvents remain, such as the comprehensive description of water per se as well as the combination of different methodological levels in composite schemes such as ONIOM,13 which cannot be considered as a straightforward task. Due to the exploratory aims of the present study, we hence focused on general trends already detectable at the continuum level, therefore sacrificing the electronic resolution of the solvent. Solvation energies were determined as the differences of the total electron energies of the systems in the gas phase (Eelg ) and embedded in the polarizable continuum (Eelaq) of the same charge state, el Esolv (M(0/ +1/ −1)) = Egel(M(0/ +1/ −1)) − Eaq (M(0/ +1/ −1))

(3)

The quantum chemical calculations were carried out using Gaussian 16.46



ASSOCIATED CONTENT

S Supporting Information *

The Supporting Information is available free of charge on the ACS Publications website at DOI: 10.1021/acsomega.9b00819. Values of vertical electron affinities, vertical electron affinities employing a polarizable continuum model, vertical ionization energies, vertical ionization energies employing a polarizable continuum model, solvation energies, vertical ionization energies for different basis set sizes, and vertical electron affinities for different basis set sizes, highest molecular orbitals of selected anionic species, and Cartesian coordinates of all molecules under consideration in this article (PDF)



AUTHOR INFORMATION

Corresponding Author

*E-mail: [email protected]. Phone: +43 512 507 52670. ORCID

Andreas Mauracher: 0000-0003-3418-789X Author Contributions

The authors contributed equally to this article. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by the Austrian Ministry of Science BMWF as part of the University infrastructure program of the scientific computing platform LFU Innsbruck. This work was supported by the Austrian Science Fund (FWF, P 30355).



REFERENCES

(1) Dans, P. D.; Walther, J.; Gómez, H.; Orozco, M. Multiscale simulation of DNA. Curr. Opin. Struct. Biol. 2016, 37, 29−45. (2) Brooks, B. R.; Bruccoleri, R. E.; Olafson, B. D.; States, D. J.; Swaminathan, S.; Karplus, M. CHARMM: A program for macromolecular energy, minimization, and dynamics calculations. J. Comput. Chem. 1983, 4, 187−217. (3) MacKerell, A. D., Jr.; Brooks, B. R.; Brooks, C. L., III; Nilsson, L.; Roux, B.; Won, Y.; Karplus, M. CHARMM: The Energy Function 10797

DOI: 10.1021/acsomega.9b00819 ACS Omega 2019, 4, 10791−10798

ACS Omega

Article

Doubles Variations of W1 Theory. J. Chem. Theory Comput. 2009, 5, 2687−2693. (44) Frisch, M. J.; Pople, J. A.; Binkley, J. S. Self-consistent molecular orbital methods 25. Supplementary functions for Gaussian basis sets. J. Chem. Phys. 1984, 80, 3265−3269. (45) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3094. (46) Frisch, M. J.; Trucks, G. W.; Schlegel, H. B.; Scuseria, G. E.; Robb, M. A.; Cheeseman, J. R.; Scalmani, G.; Barone, V.; Petersson, G. A.; Nakatsuji, H.; Li, X.; Caricato, M.; Marenich, A. V.; Bloino, J.; Janesko, B. G.; Gomperts, R.; Mennucci, B.; Hratchian, H. P.; Ortiz, J. V.; Izmaylov, A. F.; Sonnenberg, J. L.; Williams; Ding, F.; Lipparini, F.; Egidi, F.; Goings, J.; Peng, B.; Petrone, A.; Henderson, T.; Ranasinghe, D.; Zakrzewski, V. G.; Gao, J.; Rega, N.; Zheng, G.; Liang, W.; Hada, M.; Ehara, M.; Toyota, K.; Fukuda, R.; Hasegawa, J.; Ishida, M.; Nakajima, T.; Honda, Y.; Kitao, O.; Nakai, H.; Vreven, T.; Throssell, K.; Montgomery, J. A., Jr.; Peralta, J. E.; Ogliaro, F.; Bearpark, M. J.; Heyd, J. J.; Brothers, E. N.; Kudin, K. N.; Staroverov, V. N.; Keith, T. A.; Kobayashi, R.; Normand, J.; Raghavachari, K.; Rendell, A. P.; Burant, J. C.; Iyengar, S. S.; Tomasi, J.; Cossi, M.; Millam, J. M.; Klene, M.; Adamo, C.; Cammi, R.; Ochterski, J. W.; Martin, R. L.; Morokuma, K.; Farkas, O.; Foresman, J. B.; Fox, D. J. Gaussian 16; Rev. A.03, Gaussian, Inc.: Wallingford, CT, 2016.

(23) Kočišek, J.; Pysanenko, A.; Fárník, M.; Fedor, J. Microhydration Prevents Fragmentation of Uracil and Thymine by Low-Energy Electrons. J. Phys. Chem. Lett. 2016, 7, 3401−3405. (24) Poštulka, J.; Slavíček, P.; Fedor, J.; Fárník, M.; Kočišek, J. Energy Transfer in Microhydrated Uracil, 5-Fluorouracil, and 5Bromouracil. J. Phys. Chem. B 2017, 121, 8965−8974. (25) Kočišek, J.; Sedmidubská, B.; Indrajith, S.; Fárník, M.; Fedor, J. Electron Attachment to Microhydrated Deoxycytidine Monophosphate. J. Phys. Chem. B 2018, 122, 5212−5217. (26) Markovich, G.; Pollack, S.; Giniger, R.; Cheshnovsky, O. Photoelectron spectroscopy of Cl−, Br−, and I− solvated in water clusters. J.Chem. Phys. 1994, 101, 9344−9353. (27) Arnold, D. W.; Bradforth, S. E.; Kim, E. H.; Neumark, D. M. Study of I−(CO2)n, Br−(CO2)n, and I−(N2O)n clusters by anion photoelectron spectroscopy. J. Chem. Phys. 1995, 102, 3510−3518. (28) Oyler, N. A.; Adamowicz, L. Electron attachment to uracil: theoretical ab initio study. J. Phys. Chem. 1993, 97, 11122−11123. (29) Bao, X.; Sun, H.; Wong, N.-B.; Gu, J. Microsolvation Effect, Hydrogen-Bonding Pattern, and Electron Affinity of the Uracil− Water Complexes U−(H2O)n (n = 1, 2, 3). J. Phys. Chem. B. 2006, 110, 5865−5874. (30) Kim, S.; Schaefer, H. F., III Effects of microsolvation on uracil and its radical anion: Uracil·(H2O)n (n=1−5). J. Chem. Phys. 2006, 125, 144305. (31) Kumar, A.; Sevilla, M. D. Low-Energy Electron (LEE)-Induced DNA Damage: Theoretical Approaches to Modeling Experiment. In Handbook of Computational Chemistry; Leszczynski, J.; KaczmarekKedziera, A.; Puzyn, T.; Papadopoulos, M. G.; Reis, H.; Shukla, M. K., Eds. Springer International Publishing: Switzerland, 2017; pp 1741− 1802. (32) Gu, J.; Xie, Y.; Schaefer, H. F., III Electron attachment to DNA single strands: gas phase and aqueous solution. Nucleic Acids Res. 2007, 35, 5165−5172. (33) Cleare, M. J.; Hoeschele, J. D. Studies on the antitumor activity of group VIII transition metal complexes. Part I. Platinum (II) complexes. Bioinorg. Chem. 1973, 2, 187−210. (34) Varbanov, H. P.; Jakupec, M. A.; Roller, A.; Jensen, F.; Galanski, M.; Keppler, B. K. Theoretical Investigations and Density Functional Theory Based Quantitative Structure−Activity Relationships Model for Novel Cytotoxic Platinum(IV) Complexes. J. Med. Chem. 2012, 56, 330−344. (35) Rackwitz, J.; Bald, I. Low-Energy Electron-Induced Strand Breaks in Telomere-Derived DNA SequencesInfluence of DNA Sequence and Topology. Chem. − Eur. J. 2018, 24, 4680−4688. (36) HyperChem(TM) Professional 8.0; Hyperchem Inc.: 1115 Gainsville, Florida, USA. (37) Priyakumar, U. D.; MacKerell, A. D. Base Flipping in a GCGC Containing DNA Dodecamer: A Comparative Study of the Performance of the Nucleic Acid Force Fields, CHARMM, AMBER, and BMS. J. Chem. Theory Comput. 2005, 2, 187−200. (38) Kruse, H.; Banás,̌ P.; Š poner, J. Investigations of Stacked DNA Base-Pair Steps: Highly Accurate Stacking Interaction Energies, Energy Decomposition, and Many-Body Stacking Effects. J. Chem. Theory Comput. 2018, 15, 95−115. (39) Adamo, C.; Barone, V. Toward reliable density functional methods without adjustable parameters: The PBE0 model. J. Chem. Phys. 1999, 110, 6158−6170. (40) Godbout, N.; Salahub, D. R.; Andzelm, J.; Wimmer, E. Optimization of Gaussian-type basis sets for local spin density functional calculations. Part I. Boron through neon, optimization technique and validation. Can. J. Chem. 1992, 70, 560−571. (41) Curtiss, L. A.; Redfern, P. C.; Raghavachari, K. Gaussian-4 theory. J. Chem. Phys. 2007, 126, 084108. (42) Ochterski, J. W.; Petersson, G. A.; Ontgomery, J. A., Jr. A complete basis set model chemistry. V. Extensions to six or more heavy atoms. J. Chem. Phys. 1996, 104, 2598−2619. (43) Barnes, E. C.; Petersson, G. A.; Montgomery, J. A., Jr.; Frisch, M. J.; Martin, J. M. L. Unrestricted Coupled Cluster and Brueckner 10798

DOI: 10.1021/acsomega.9b00819 ACS Omega 2019, 4, 10791−10798