Inverse Temperature Dependence of Nuclear Quantum Effects in DNA Base Pairs.

Despite the inherently quantum mechanical nature of hydrogen bonding, it is unclear how nuclear quantum effects (NQEs) alter the strengths of hydrogen bonds. With this in mind, we use ab initio path integral molecular dynamics to determine the absolute contribution of NQEs to the binding in DNA base pair complexes, arguably the most important hydrogen-bonded systems of all. We find that depending on the temperature, NQEs can either strengthen or weaken the binding within the hydrogen-bonded complexes. As a somewhat counterintuitive consequence, NQEs can have a smaller impact on hydrogen bond strengths at cryogenic temperatures than at room temperature. We rationalize this in terms of a competition of NQEs between low-frequency and high-frequency vibrational modes. Extending this idea, we also propose a simple model to predict the temperature dependence of NQEs on hydrogen bond strengths in general.

It has been said that without hydrogen bonds (HBs), all wooden structures would collapse, cement would crumble, oceans would vaporize, and all living things would disintegrate into inanimate matter.[1] While the concept of the HB dates back to at least the 1920s[2] and is now well-defined,[3] the small mass of the proton means that HBs are intrinsically quantum mechanical and that zero-point energy and tunneling can be of critical importance. Quantum fluctuations involving HBs are crucial, for example, in biological processes such as DNA tautomerization[4−6] and enzyme reactions.[7−12] It is also known that n class="Chemical">hydrogendated and deuterated chemicals can have different biochemical potencies, a fact that is now enthusiastically being exploited within the pharmaceutical industry through the development of novel classes of deuterated drugs. Nonetheless, fundamental understanding of the quantum nature of HBs is far from complete, with a general understanding of how and to what extent nuclear quantum effects (NQEs) influence the strengths of HBs yet to be established. Given that strength is arguably the most important property of any bond, this seems to represent a fairly significant gap in understanding.

Indirect information about the role of NQEs on HB strengths can be made through isotopic substitution experiments. These have shown that upon replacing hydrogen with n class="Chemical">deuterium, the lengths of HBs can change, a so-called secondary geometric isotope effect.[13,14] HBs can get shorter or longer depending on the material, and the extent of the change can vary from one material to the next and with temperature (see, e.g., refs (15−18)). Interestingly, very small changes in structure can lead to large variations in physical properties. For example, the lattice constants of the hydrogenated and deuterated versions of the ferroelectric material squaric acid differ by 1%, yet their ferroelectric to paraelectric transition temperatures differ by ∼200 K.[16,17] Secondary geometric isotope effects such as these have been explained with the help of theory and simulation,[18−23] and notably, it has been argued that the direction and extent of the change upon isotopic substitution depends on the strength of the HB.[18] This, in turn, has been rationalized with a theory of competing quantum effects where it is said that quantum delocalization along the HB helps to shorten the bond but delocalization out of the plane acts to lengthen it.[19,22,24,25] Indeed, this concept has proved to be useful in explaining a host of phenomena in, for example, liquid water, ice, and biomolecules[25−31] and has seen recent experimental verification.[32] Nonetheless, most work to date has focused on geometrical properties, and direct information on how and to what extent NQEs influence the strengths of HBs[33] is desirable. Although the total contribution of NQEs to HB strength is likely to be small, small energies are often important. This is particularly true in biology, where structures and processes are governed by a delicate balance of interactions.[31] The cost to unzip double-stranded DNA in solution is, for example, only ∼20–100 meV (∼0.5–2 kcal/mol or 1–4kBT) per base pair duplex at room temperature.[34−36] Similarly, the melting temperature of DNA strands is so exquisitely sensitive that substituting just one out of a thousand base pairs leads to a measurable change in melting temperature.[37]

In this work, we use computer simulations to directly evaluate the influence of NQEs on the binding free energy of DNA base pair complexes. These hydrogen-bonded complexes are not only crucial to life but have also gained significant interest in nanotechnology.[38] We specifically examined Watson–Crick n class="Chemical">hydrogen-bonded base pair complexes of adenine-thymine (AT) and cytosine-guanine (CG) in the gas phase. The focus is on understanding how quantum effects alter the duplex hydrogen-binding interaction in the dimers; this is, of course, an integral interaction to DNA binding and is, for example, the key parameter in nearest-neighbor[36] and coarse-grained models[39] of DNA binding. Although stacking interactions and solvent effects are relevant to the unzipping and melting of real double-stranded DNA, by focusing exclusively on the duplex hydrogen-bonding interaction, we are able to unambiguously understand the role played by NQEs. The particular computational techniques employed involve density functional theory (DFT) for a description of the potential energy surface in conjunction with path integral molecular dynamics (PIMD), which together enable equilibrium thermal properties including NQEs to be rigorously accounted for. This methodology has been widely applied to tackle a host of chemical problems in the gas and condensed phases (see, e.g., refs (17), (18), (29−31), (33), and (40−46)). Furthermore, by combining PIMD with thermodynamic integration, we can explicitly calculate the binding free energy change upon transforming the system from one containing classical nuclei to one containing quantum nuclei. With this approach, we find, at room temperature, that NQEs increase the interaction strength in both complexes by ∼0.5 kcal/mol or ∼1kBT. At 100 K, a temperature appropriate for preserving DNA information and a temperature at which NQEs are generally expected to be more significant, we find that the impact of NQEs on the binding energy is smaller. Analysis reveals that this surprising temperature dependence arises from a competition of quantum effects associated with low- and high-frequency vibrational modes. Extending our findings from DNA base pairs, we use our physical understanding of competing quantum nuclear effects to propose a simple model to estimate the temperature dependence of NQEs on binding free energies of hydrogen-bonded complexes in general.

Our PIMD simulations were performed with the CP2K[47,48] code connected to the i-PI wrapper.[49] A full account of the computational setup is given in Supporting Information (SI) section II, and here, we outline the key features. The PIGLET thermostat[50] was used for an efficient sampling of the imaginary time path integrals. At 300 (100) K, 6 (16) replicas of the systems were taken to sample the imaginary time path integral, which has been shown to yield converged quantum kinetic energies.[50,51] Molecular dynamics trajectories were generally 10 ps long, which we show in the SI (section II) are sufficiently long to obtain converged binding free energies. For DFT, we used the optB88-vdW functional,[52] which is a revised version of the van der Waals density functional (vdW-DF) of Dion et al.[53] This functional is particularly appropriate for the base pair duplexes under consideration as it yields very accurate interaction energies for the complexes in comparison to quantum chemistry reference calculations (see SI section II). In the SI (section III), we also report results from the hybrid PBE0 functional[54] for the binding free energy change of the AT base pair. Within the statistical error bars, the results obtained with PBE0 and optB88-vdW are the same at room temperature.

The well-known structures of the Watson–Crick base pair complexes are shown in Figure . From this, it can be seen that the AT complex is held together by two HBs (an NH–O and an NH–N bond), whereas the CG complex is held together by three (two NH–O bonds and an NH–N bond). The HBs have a range of lengths, with both the NH–O and NH–N bonds varying from 1.7 to 1.9 Å in the ground-state (geometry-optimized) structure.

Figure 1

Structures of the Watson–Crick AT and CG base pairs. Black: carbon; red: oxygen; blue: nitrogen; white: hydrogen.

Structures of the Watson–Crick AT and CG base pairs. Black: n class="Chemical">carbon; red: oxygen; blue: nitrogen; white: hydrogen.

As a first step to understand the role of NQEs, we ran a set of ab initio molecular dynamics (AIMD) simulations as well as a set of ab initio PIMD simulations. We concentrated on room temperature (300 K) and a cryogenic temperature (100 K); room temperature is of obvious interest to biology, and cryogenic temperatures are relevant to, for example, DNA preservation and DNA-based devices.[56,57] At both temperatures, the dimers remain n class="Chemical">hydrogen-bonded, and no intermolecular proton transfer is observed. Simply by comparing the structures obtained from the simulations with the classical and quantum nuclei, we can gain an initial indication of the role that NQEs play. Interestingly, we find that different HBs respond in a different manner to the inclusion of NQEs; some HBs get longer, some get shorter, and some remain unchanged. Previously, it was shown for a range of systems that how a HB responded to the inclusion of NQEs depended on its strength, with relatively strong HBs becoming shorter and relatively weak HBs becoming longer.[18] In Figure a, we explore this issue for DNA base pairs by plotting how the intermolecular separation (specifically the N–O and N–N heavy atom distances) changes upon going from classical to quantum nuclei. HB strength is estimated with a simple and standard criterion involving the red shift in the harmonic stretching frequency of the covalent NH bond involved in the HB.[18,58,59] The larger the red shift of this stretching frequency, the stronger the HB. We find for the individual HBs in the DNA complexes considered that the correlation seen before also holds here; the strong HBs do indeed tend to be shortened while the weak ones tend to be elongated by NQEs. Specifically, at 300 K, in the A-T base pair, the stronger NH–N bond becomes shorter in the PIMD simulations, and the weaker NH–O bond barely changes. In the C-G base pair at 300 K, the weakest NH–O bond becomes longer in the PIMD simulations, while the other two HBs become shorter. At 100 K, the correlation also holds, and overall, the HB lengths change in a similar manner to what is observed at 300 K. It is clear, therefore, that NQEs impact the intermolecular separation, and this could be probed experimentally through, for example, isotopic substitution measurements of secondary geometric isotope effects. However, from the structural data alone, it remains unclear how NQEs have affected the interaction strength between the dimers in each complex.

Figure 2

(a) Differences between the heavy atom separation distances from PIMD and MD simulations. Positive changes mean that the N(H)–O or N(H)–N bonds are longer in the PIMD than those in the AIMD simulations; and negative values mean that they are shorter in PIMD. The five different HBs in the base pairs are arranged from left to right in order of decreasing strength, with strength characterized by the harmonic frequency of the N–H stretch in the HB divided by the harmonic frequency of the N–H stretch in the monomers.[18] A snapshot of the AT base pair taken from a PIMD simulation is also shown in the inset; each sphere is a “bead” in the PIMD simulation. (b) Plot of the binding free energy change due to NQEs (eq ) in the AT (blue) and CG (red) base pairs obtained from PIMD. A negative binding free energy change means that NQEs strengthen the binding, while a positive binding free energy change means that NQEs weaken the binding. Also shown with the dashed lines are the predictions of each base pair obtained within the harmonic approximation. The error bars in (a) and (b) have been calculated using block averaging.[55]

In order to unambiguously determine how NQEs alter the interaction strength within the complexes, we computed how the binding free energy changes upon going from classical to quantum nuclei. To this end, we employed a thermodynamic integration scheme, previously used in calculations of isotope effects.[60−64] Full details of the specific approach used here are given in section I of the SI. The key point is that we obtain a binding free energy change by performing a thermodynamic integration in which the mass of the nuclei is switched from classical to quantum (c → q). Specifically, the classical to quantum change in the binding free energy, ΔFbc→q, is given byHere ⟨K⟩ is the ensemble average of the quantum kinetic energy, which can be directly obtained from PIMD simulations, and g is a mass-dependent integration variable. Separate PIMD simulations must be performed for the two isolated molecules, M1 and M2, and the hydrogen-bonded complex M1:M2, and in order to obtain accurate values for the integrand, simulations must be performed for several values of g (7 in this study). In total, to calculate ΔFbc→q for a single system at a single temperature, trajectories equivalent to ∼2 ns of AIMD simulations must be accumulated. This enormous computational cost is a key reason that binding free energies have rarely been computed with ab initio PIMD.

The results obtained from the thermodynamic integration of our ab initio PIMD simulations are shown in Figure b at the two temperatures considered. A negative ΔFbc→q means that NQEs strengthen the binding, while a positive value means that NQEs weaken it. At room temperature, the binding of both the AT and CG complexes is strengthened by ∼20 meV (0.5 kcal/mol) when NQEs are accounted for, specifically, 24 ± 4 meV (0.55 ± 0.1 kcal/mol) for AT and 17 ± 4 meV (0.39 ± 0.1 kcal/mol) for n class="Chemical">CG. On an absolute scale of chemical bonding interactions, these are, of course, small energies. However, as we know in biology, energies tend to be finely balanced and very small changes in energy can be critical. For example, 20 meV is on the same scale as thermal energy at room temperature and similar to the cost to unzip double-stranded DNA in solution (estimates range from 20 to 100 meV at room temperature).[36]

Upon considering how NQEs alter the energetics when the temperature decreases from 300 K to a cryogenic temperature, one would expect the influence to be greater than that at room temperature. However, the opposite is the case, with the contribution of NQEs to the binding being smaller at 100 K for both base pair complexes. In the AT complex, NQEs strengthen the binding by 10 ± 4 meV (0.23 ± 0.1 kcal/mol), about half the value at 300 K. In the CG complex, the contribution of NQEs is even less (7 ± 4 meV or 0.16 ± 0.1 kcal/mol), and moreover, NQEs now act to destabilize the complex. Thus, although simple assumptions about the temperature dependence of NQEs have been useful in understanding, for example, structural properties of liquid n class="Chemical">water,[65] the same cannot be done when it comes to binding free energies.

Can we understand the free energy changes obtained? In the PIMD simulations, the free energies arise from thermal sampling of the quantum kinetic energy through a complex interplay of vibrational modes, which is not particularly straightforward to interpret. It is possible to project the quantum kinetic energy on to particular modes (see, e.g., refs (25) and (30)); however, here we opt to perform an analysis within the harmonic approximation, which provides a relatively straightforward means of establishing how specific groups of vibrational modes contribute to the observed changes in free energies. Within the harmonic approximation, the quantum kinetic energies in eq are calculated fromwhere ω are the 3N harmonic vibrational frequencies (including the zero frequency translation and rotation modes) of the ground-state geometry-optimized complexes, ℏ is the reduced Planck’s constant, and β is the inverse temperature.[66] Results from the harmonic approximation are shown in Figure . Clearly, the harmonic approximation does not reproduce PIMD exactly; anharmonic effects are certainly important in these systems. Nonetheless, however, the harmonic approximation does reasonably well at a qualitative level; for both the AT and CG complexes, the harmonic approximation gets the correct sign of the change and the correct temperature dependence. With the picture of competing quantum effects in mind, we explored if it could also be used to explain the observed changes in binding free energies. To this end, we define a separation between high-frequency and low-frequency vibrational modes; 2000 cm–1 is chosen as this represents a threshold between high-frequency covalent bond stretching modes and relatively low frequency bond bending and collective intermolecular stretching modes (see Figure S3 in the SI). As shown on the right of each plot in Figure , the high- and low-frequency vibrational modes have quite large (∼40–90 meV) but opposing contributions to the overall binding free energy change; the low-frequency modes reduce the binding free energy, whereas the high-frequency modes increase it. In the SI (section IV), we show the integration curves from which the histograms in Figure have been obtained, which show precisely how the observed changes in the quantum kinetic energy arise from competing contributions in the two vibrational regimes. In simple terms, this behavior arises because the high-frequency modes tend to be softened upon HB formation, which reduces the zero-point energy, whereas the low-frequency modes are hardened or new ones are created, which tends to increase the zero-point energy and therefore reduce the effective attraction. Overall, it is clear that the net impact of NQEs on DNA binding results from a significant cancellation of two larger quantum contributions and that the picture of competing quantum effects can be used to quantitatively understand how NQEs alter HB strengths.

Figure 3

Competing quantum effects and explanation of the anomalous temperature dependence. The binding free energy change (ΔFbc→q) obtained from the PIMD simulations is compared with results from the harmonic approximation (Harm.). The binding free energy changes within the harmonic approximation are also decomposed into high- (ωhigh) and low-frequency (ωlow) contributions, revealing that the net change in binding free energy arises from a significant cancellation of contributions from these two regimes. The unusual temperature dependence simply arises because of a greater cancellation of terms at 100 K (black bars) than at 300 K (red bars). The change with temperature is more pronounced for the contribution from the low-frequency modes than it is for the high-frequency modes.

Let us now move on to the seemingly anomalous temperature dependence of NQEs on the binding free energies. Having established that the overall change in binding free energy arises from a cancellation of two opposing effects, one can recognize that there is a greater cancellation of terms at 100 K than there is at 300 K (Figure ). Looking at this figure more closely, we can see that as the temperature increases from 100 to 300 K, the contributions to the binding free energy differences from both the low-frequency and high-frequency modes decrease. This is to be expected and is consistent with conventional understanding that quantum effects are less prominent at high temperatures. However, the change with temperature is more pronounced for the low-frequency modes than it is for the high-frequency modes. These changes are governed by occupation of the vibrational modes through the relation kBT/ℏω, and within the temperature regime being considered, the population of the high-frequency modes changes much less than that of the low-frequency modes. Hence, it is the underlying competition coupled with the differing temperature dependence of the high- and low-frequency modes that makes the net impact of NQEs on the binding free energies more significant at 300 K than that at 100 K. It is interesting to note that differences in the temperature dependence of the intermolecular and intramolecular modes have also been used to provide qualitatively the same explanation for isotope fractionation in water, in particular, to explain an interesting inversion at high temperatures in the liquid n class="Chemical">water/vapor fractionation ratio.[25,30]

The temperature dependence of NQEs is of relevance beyond the HBs in DNA, and we now show how a minimalistic model of HB formation can be used to make predictions about the temperature dependence of NQEs in HBs in general. In the model, we mimic the essence of the competing quantum effects with only two variables, ωlow and ωhigh, which represent the total red shift of the high-frequency modes and the total blue shift of the low-frequency modes, respectively. (No explicit HB parameters, i.e., bond length or model potentials, are introduced here.) Upon using these two variables to compute the change in quantum kinetic energy (see section V of the SI), one can predict whether NQEs strengthen or weaken the binding of a hydrogen-bonded system at a given temperature. In Figure , we show how the transition from NQEs strengthening HBs to NQEs weakening HBs depends on the interplay of these modes at room temperature (solid line) and at 100 K (dashed line); a broader range of temperatures is reported in the SI (section V). Also included in Figure are the results for some specific n class="Chemical">hydrogen-bonded dimers in which all of the modes are taken into account. The particular dimers considered include the two DNA base pair complexes as well as water, HF, formamide, and formic acid dimers. The model presented here is incredibly simple; for example, it neglects changes in modes at intermediate frequencies.[31] However, it produces qualitatively correct results. It correctly places the base pairs just within the strengthening regime at room temperature and (in agreement with full harmonic frequency calculations) shows that the water, HF, and formamide dimers are weakened at room temperature. It also correctly captures the behavior observed for the CG dimer wherein NQEs switch from strengthening to weakening upon lowering the temperature. As well as CG, the formic acid dimer also behaves in a similar manner, revealing that CG is not any sort of a special case and that other hydrogen-bonded systems could exhibit similar behavior. Although ab initio PIMD simulations are becoming increasingly affordable computationally,[49] on the whole, they remain expensive and are far from routine. However, the conceptual framework presented here allows for ballpark predictions to be made of the role of NQEs in HBs on the basis of (harmonic) vibrational frequencies. This is data that can be obtained from experiment or relatively cheap and easy vibrational frequency calculations.

Figure 4

Temperature at which NQEs switch from weakening to strengthening the binding for the model hydrogen-bonded system, plotted as a function of total high (ωhigh) and total low (ωlow) frequency mode shifts. The solid line marks the 300 K transition, whereas the dashed line marks the 100 K transition. Actual frequency changes corresponding to six specific dimers are also indicated on the figure; these data points correspond to the average changes per HB for frequencies computed within the harmonic approximation. At room temperature, the AT and CG base pairs and the formic acid dimer are in the strengthening regime, and the water dimer, HF dimer, and formamide dimer are in the weakening regime.

Temperature at which NQEs switch from weakening to strengthening the binding for the model n class="Chemical">hydrogen-bonded system, plotted as a function of total high (ωhigh) and total low (ωlow) frequency mode shifts. The solid line marks the 300 K transition, whereas the dashed line marks the 100 K transition. Actual frequency changes corresponding to six specific dimers are also indicated on the figure; these data points correspond to the average changes per HB for frequencies computed within the harmonic approximation. At room temperature, the AT and CG base pairs and the formic acid dimer are in the strengthening regime, and the water dimer, HF dimer, and formamide dimer are in the weakening regime.

Before concluding, it is useful to put the results and model obtained in the current study into context. The model presented here focuses on energetics, and therefore, it complements other successful HB models, such as the recently proposed diabatic model of McKenzie and co-workers.[21−23,67] We note that these models have been very useful in predicting geometric isotopic effects, proton transport, and frequency shifts, and therefore, together with the current model, a broader picture of the structure, dynamics, and energetics of HBs is emerging. Here, we have obtained free energy changes upon moving from the classical to quantum regime and have established the absolute contribution from NQEs. This is not a measure that is readily accessible in experiments where indications about the importance of quantum nuclear effects are obtained indirectly through isotopic substitution measurements. With this in mind, we note that we have also calculated the H/D isotope effect on the binding free energy in the harmonic approximation. From this, we find that at room temperature, the binding free energy of the AT base pair decreases by 3 meV upon deuteration and the binding of the CG base pair increases by 6 meV upon deuteration. Thus, as expected, the difference between H and D does not capture the full contribution of NQEs. Finally, we note again that biological environments will be more complex than the gas-phase systems considered here and that the presence of a solvent will no doubt have an impact on the nature of the HBs. For example, experiments have shown that the proton-transfer rate for base pair complexes varies with the dielectric constant of the solvent.[5,68,69] Our model does not directly address measurements such as these, but we do expect the physical insight obtained here to hold in more complex environments, and, to first order, estimates of the influence of a solvent could be made by examining how the solvent alters the vibrational frequencies within the n class="Chemical">hydrogen-bonded complexes being considered.

To conclude, we have reported what is, to the best of our knowledge, the first determinations of the quantum contribution to the binding free energy of DNA base pairs. We have found that NQEs strengthen the binding of both AT and CG complexes at room temperature. At a lower temperature, however, NQEs have a smaller impact on the binding free energies; analysis of the quantum kinetic energies in each system reveals that this seemingly anomalous temperature dependence arises from a balance between competing quantum effects associated with low-frequency and high-frequency modes of vibration. Upon forming a HB, the high-frequency (covalent bond stretching) modes are softened; hence, quantum kinetic energy is lost, and the system is stabilized. This stabilization, however, is offset by the quantum kinetic energy gained when low-frequency modes are hardened or created upon forming the HB. This shows that the picture of competing quantum effects can be applied to understand how NQEs alter the energetics of n class="Chemical">hydrogen bonding, and with this insight, we have presented a simple model to estimate the temperature dependence of NQEs in hydrogen-bonded systems in general. Of course, real DNA is much more complex than the simple dimers considered here, but at the very least, the current study demonstrates that the role played by quantum effects could be more significant than previously anticipated and deserves further study.