The influence of base pair tautomerism on single point mutations in aqueous DNA

4.1. Ensemble-based classical molecular dynamics

The reproducibility of MD remains a topic of debate for many theoretical chemists and biologists [60–62]. Because of the extreme sensitivity to the initial conditions in any chaotic MD simulation, one-off simulations are not reliable. Instead, we can obtain statistically robust results by performing ensemble-based simulations, that is, a collection of n replicas each differing from the other solely in terms of the initial velocities assigned to all the atoms, drawn from a Maxwell–Boltzmann distribution at the temperature of interest [62]. Furthermore, ensemble-based simulations provide a reliable means of quantifying uncertainty in general [63].

The parmbsc1 AMBER force field is used since it performs well for solvated double-helix DNA MD simulations [64,65]. The B-DNA structure of the Drew–Dickerson dodecamer (PDB ID. 1BNA) was neutralized (22 Na⁺) and solvated in a water box (dimensions: 71.15 Å × 73.13 Å × 85.94 Å) with the TIP3P water model [59,66]. Ensemble classical MD was then performed under periodic boundary conditions. The cut-off for the interaction distance for both electrostatic and van der Waals calculations was set to 10 Å. To prevent discontinuities in electrostatic and van der Waals energies, pairs of atoms greater than 11.5 Å apart were excluded. Periodic boundary electrostatics were calculated using the particle mesh Ewald method with a grid spacing of 1 Å. The following equilibration and simulation steps are repeated 10 times using the parmbsc1 AMBER force field, generating a total of 100 ns simulation for the ensemble.

Equilibration: The coordinates of the entire DNA were restrained as the geometry of the rest of the system was minimized using the conjugate gradient and line search algorithm for 10 000 steps. Next, the temperature of the system was incrementally raised from 50 K to 300 K over a time period of 30 ps with a time step of 1 fs. The temperature of the system was maintained constant using a Berendsen barostat at a pressure of 1 atm. The restraints were then systematically removed over 0.5 ns, followed by an unrestrained 0.5 ns run.

Simulation: The production runs were for 10 ns, simulated at a constant temperature and pressure (300 K, 1 atm). All bonds between heavy atoms and hydrogen were constrained to their nominal length during integration and a time step of 2 fs was employed.

All MD simulations were performed using NAMD 2.12 on the UCL high-performance computing (HPC) facility Grace.

4.2. Choosing the quantum-mechanical method

A group of QM methods were compared on their ability to accurately describe base pair geometries and their interaction energies. Various QM methods that are known to accurately describe proton transfer barriers were chosen [67]. The QM methods assessed include B3LYP, CAM-B3LYP, LC-ωPBE, M06-2X and MP2, in conjunction with different basis sets and dispersion corrections. Both double zeta and triple zeta basis sets were assessed in conjunction with a double diffuse function to properly evaluate the energies of hydrogen bonds. To accurately capture weak dispersion interactions between base pairs, the Grimme-D3 correction and the exchange-hole dipole model (XDM) dispersion correction were individually assessed [68–70]. All QM calculations were performed using NWChem 6.6 on the Blue Waters supercomputer at the National Center for Supercomputing Applications, USA.¹

The geometry of a gas phase G:C base pair was optimized using each QM method and the guanine–cytosine interaction energies were calculated. The energies were then compared with highly accurate coupled-cluster values obtained by Hobza and co-workers [71] (as shown in figure 4). The RMSD of the selected base pair geometries were then compared with an MP2/aug-cc-pvdz optimized G:C base pair (further detailed in electronic supplementary material, S2). Although it is not feasible to use a triple zeta basis set, owing to its computational cost, both the double and triple zeta basis sets were considered in this benchmark for comparison purposes. In general, the triple zeta basis sets are expected to outperform the double zeta basis sets. However, in figure 4, it is shown that this is not the case for the Dunning basis sets in conjunction with the MP2 and CAM-B3LYP method. The same trend can be seen for the Pople basis sets, whereby 6-31++G** outperforms 6-311++G** when used with the B3LYP, LC-ωPBE and MP2 methods. This may be due to the cancellation of errors.

• Download figure
• Open in new tab
• Download PowerPoint

The interaction energies of stacked DNA base pairs were evaluated using B3LYP in combination with different basis sets and dispersion correction schemes. The same protocol as outlined in Šponer and co-workers [72] was followed and the interaction energies for 10 different combinations of stacked base pairs were calculated. A detailed definition of how the interaction energies are defined in addition to the computational methodology applied is offered in electronic supplementary material, S2.2. The results of the benchmark were then compared with the highly accurate DLPNO-CCSD(T)/CBS reference values (table 1). It is demonstrated that the B3LYP+XDM/aug-cc-pvdz binding energies compare well to the reference values for both the two-body and four-body stacking interaction energies. Despite its remarkable accuracy in predicting binding energies of hydrogen-bonded base pairs, CAM-B3LYP fails to accurately predict the binding energies of stacked base pairs.

Therefore, the B3LYP/aug-cc-pvdz with XDM dispersion correction was selected as the QM method in this study, based on its ability to accurately reproduce the binding energies and geometries of both the hydrogen bonded and stacked base pairs (as shown in electronic supplementary material, S2).

4.3. Ensemble quantum mechanics/molecular mechanics

Previous base pair proton transfer studies that use ab initio MD are limited by the expense (and therefore the accuracy) of the QM method, by the size of the explicit base pair environment and by the accessible time scale of simulation [44,73,74]. The impracticality and cost of performing multiple ab initio MD simulations is obviated by the use of a multiscale ensemble based on the QM/MM method [75]. In this way, the statistical relevance of various proton transfer pathways can be assessed.

Despite numerous studies of DNA using QM/MM methods, there is no single agreed approach to modelling it [10,48–52]. Previous work has shown that the effects of adjacent base pair stacking can modify both the strengths and lengths of base pair hydrogen bonds, which ultimately influence the proton transfer energy profile [30,44,50]. However, a more recent QM/MM study (performed by Das et al. [52] in which 10 snapshots from a single MD simulation were studied) showed that adjacent base pair stacking has a relatively small influence on the base pair proton transfer energetic profile of the order of 1 kcal mol⁻¹. The authors investigated the effect of neighbouring base pair polarization by systematically increasing the QM region size from one base pair through to five nucleosides. Of the 10 QM/MM calculations performed, the authors reported the electronic reaction energy for the single proton transfer in the G⁺:C base pair to vary between $7\text{ and }14\hspace{0.17em}\text{kcal}\hspace{0.17em}{\text{mol}}^{-1}$. The relatively large variance they observed (approx. ±3 kcal mol⁻¹) indicates that the mean reaction energy is largely independent of (i) the QM region size and (ii) the polarization effects of the adjacent base pairs. In other words, the effect of using a larger QM region lies within the uncertainty of the single proton transfer itself.

We also independently considered the benefits that might accrue from increasing the QM region size to include the adjacent nucleotide base pairs. Our findings corroborate those of Das et al. In particular, while barrier heights from our one-off QM/MM calculations are lowered slightly, they lie within the errors arising from ensemble averaging of single base pair QM regions.

The initial configurations for the subsequent QM/MM ensemble were selected from the prior classical MD simulations based on the distance between the base pairs (the distribution of these distances is shown in figure 5). By sampling from the average of this base pair distance distribution, n configurations are selected as starting points for the ensemble QM/MM study. The number of QM/MM replicas (n) is determined in §4.4.

• Download figure
• Open in new tab
• Download PowerPoint

The QM/MM simulations were performed using ChemShell 3.7 to link NWChem 6.6 (QM) with DL-POLY (MM). All of the QM/MM tasks and optimizations were performed using the DL-FIND module as implemented in ChemShell [76]. The periodicity of each initial QM/MM configuration was removed, which resulted in a solvation sphere of 15 Å containing approximately 9000 atoms. The QM region consisted of a single base pair with hydrogen linker atoms placed between the deoxyribose C1′ and the corresponding terminal N of the nucleobase. For this work, the electrostatic interactions between the QM and MM regions were modelled using electrostatic embedding. In principle, the use of a polarizable force field might be expected to provide a more realistic charge description at the price of introducing further fitting parameters. The use of such force fields in QM/MM calculations alters energetic calculations by up to ca 1 kcal mol⁻¹ [77], while established electrostatic embedding techniques are more robust and computationally efficient than the best currently available counterparts [78]. In accordance with a previous ChemShell QM/MM DNA study [79], our ‘active’ MM region consists of all residues within 15 Å of the QM base pair; the remaining residues beyond this distance were frozen. This large active MM region was selected to permit structural changes in the phosphate backbone and the inner solvation sphere during geometry optimization. The proton transfer reaction pathways were calculated between the QM/MM optimized base pair and the proton transfer product using the adiabatic climbing image nudged elastic band (CI-NEB) method [80]. Depending on the curvature of the reaction pathway, the replica was then categorized as either a stepwise or concerted process. Transition states on the reaction pathway were further optimized using the dimer method [81] and verified by a single imaginary frequency in the Hessian. When appropriate, the geometry of the reaction pathway intermediate was optimized to a local minimum. The Hessian for the QM region of each geometry in the optimized proton transfer reaction pathway was then calculated to include thermal corrections of 300 K. Next, the rate coefficients for proton transfer were calculated using conventional transition state theory (equation (3.8)) as implemented in the ChemShell rate module [82]. All ChemShell QM/MM calculations were performed using the UK national supercomputer ARCHER.²

4.4. Uncertainty quantification

To establish an appropriate number of replicas n constituting the QM/MM ensemble, the bootstrap statistical method was employed [63]. By applying the QM/MM methodology as described in §4.3, the reaction energy (ΔE_rxn) for the G:C → G*:C* tautomerism was calculated per replica. Figure 6 shows that the mean bootstrap standard deviation (σ) for the reaction energy (ΔE_rxn) decreases as the number of replicas (n) increases. After 25 QM/MM replicas are included, the change in σ is negligible and, for this reason, n is set to 25.

• Download figure
• Open in new tab
• Download PowerPoint

5.1. G:C tautomerism

The stepwise process is the most probable double proton transfer mechanism in G:C, occurring in 84% of the replicas from the QM/MM ensemble. This process proceeds via two transition states [(G:C)${\hspace{0.17em}}_{\mathrm{a}}^{‡}$, (G:C)${\hspace{0.17em}}_{\mathrm{b}}^{‡}$] and an intermediate (G:C)_Int. The less frequent double proton transfer mechanism in G:C is concerted, occurring only 12% of the time. This process proceeds via a single transition state (G:C)^‡ and has a larger electronic energy barrier of approximately 2 kcal mol⁻¹ compared with the stepwise process. The optimized double proton transfer reaction pathways are displayed in figure 7. Tolosa et al. [83] were among the rare groups to report more than one possible reaction pathway for the double proton transfer reactions in DNA base pairs. They modelled a microhydrated G:C base pair using M06-2X/6-311++G** and steered MD. They reported stepwise, concerted and water-assisted mechanisms. The forward barrier height for the first step of the stepwise mechanism is 17.98 kcal mol⁻¹, which is well out of the range we report. In addition, the barrier height they reported for the concerted mechanism was more than twice as high as our reported values.

• Download figure
• Open in new tab
• Download PowerPoint

The energetics, kinetics and thermodynamics for the stepwise double proton transfer mechanism are shown in the left column of table 3. The standard deviation for the mean electronic and Gibbs free energies are consistently approximately 1 kcal mol⁻¹ for each step of the stepwise reaction coordinate. As shown in figure 7 (left), the energy of each part in the reaction coordinate lies within the standard deviation error of all the other points.

On the other hand, only three replicas from the QM/MM ensemble displayed the concerted G:C$\hspace{0.17em}\rightleftharpoons \hspace{0.17em}$G*:C* tautomerization pathway. The energetics, kinetics and thermodynamics for this pathway are shown in the middle column of table 3. Only three replicas participated in the concerted pathway, indicating the rarity of this phenomenon. As a result, the standard deviation of the energy barrier ($\mathrm{\Delta }{E}^{‡}$) for the concerted process is twice that of the stepwise.

It is important to note that the overall shape of the double proton transfer reaction profiles, shown in figure 7, changes substantially between electronic and Gibbs free energy. The transition state energies are stabilized by approximately 3.5 kcal mol⁻¹ in terms of the Gibbs free energy. In the case of the stepwise pathway (figure 7, left), this results in the transition states no longer being the maximum relative Gibbs free energies. Consequently, the overall instability of the G*:C* tautomer is reinforced as it occupies the highest Gibbs free energy state on the reaction coordinate. This is due to the zero-point energy contribution stabilizing both stepwise transition states ($\mathrm{\Delta }{G}_a^{‡}$ and $\mathrm{\Delta }{G}_b^{‡}$) by approximately 2 kcal mol⁻¹ more than the intermediate and products. The computed barrier heights using B3LYP with a dispersion correction are expected to be slightly underestimated (although not to the extent of using B3LYP alone) [84]. The reverse barrier heights for the concerted mechanism in our study are 1.0 ± 1.5 kcal mol⁻¹. These are in good agreement with the values of approximately 1 kcal mol⁻¹ computed using BP86/6-311++G** [42]. Other studies, however, reported reverse free energy barriers less than 0 kcal mol⁻¹ using MP2/aug-cc-pVTZ and B3LYP/6-311++G** [13,15].

The range of the rate coefficients for each step in the G:C base pair stepwise double proton transfer is displayed as a normalized histogram in figure 8. The rate coefficients fit well to the normal distribution of the data, reinforcing the calculated means and standard deviations. The forward rate coefficient for the first step (k_a), shown in orange, is approximately eight orders of magnitude lower than that of the reverse (k′_a), shown in red. The reverse of the first step (G:C)_Int → G:C is favoured, establishing an equilibrium constant (K_a) of the order of 10⁻⁸. The values for the forward (k_b) and reverse (k^′_b) rate coefficients, coloured blue and green, respectively (figure 8b), are shown to partially overlap around 10¹³ s⁻¹. The second step of this process (G:C)${\hspace{0.17em}}_{\mathrm{Int}}\hspace{0.17em}\rightleftharpoons \hspace{0.17em}$G*:C* lies close to equilibrium, with a K_b of approximately 0.6. By using equation (3.13), the overall stepwise double proton transfer equilibrium constant (K) is calculated to be of the order of 10⁻⁹. Recent computational studies have predicted the G:C double proton transfer equilibrium constant to be between ${10}^{-6}\text{ }\mathrm{and}\hspace{0.17em}{10}^{-9}$ but this is largely dependent on the different QM approximations used [13–16,83,85]. Our results for K (approx. 10⁻⁹), are in agreement with Céron-Carrasco & Jacquemin (4.22 × 10⁻⁹) [85], who also used QM/MM methods to simulate proton transfer within a d(GGG)₂ codon. They modelled the central base pair using M06-2X/6-311G**, the stacking base pairs using M06-2X/6-31G* and the sugar–phosphate backbone with the semi-empirical PM6 method. Despite predicting a similar equilibrium constant, they calculated the respective forward and reverse Gibbs free energy barrier to be 19.54 kcal mol⁻¹ and 8.12 kcal mol⁻¹, which are both approximately 8 kcal mol⁻¹ larger than our results. Céron-Carrasco et al. [50] have shown that the inclusion of base pair stacking and explicit solvation causes the double proton transfer equilibrium constant to be smaller than that in the gas phase. In corroboration with their findings, we report equilibrium constants to be one order of magnitude smaller than recent QM-only gas phase studies that have used B3LYP/6-311++G** (1.68 × 10⁻⁸) and MP2/aug-cc-pvTZ (7.5 × 10⁻⁸) methods [13,15].

• Download figure
• Open in new tab
• Download PowerPoint

The rate coefficient of the first forward step (k_a′) is the slowest by several orders of magnitude and, therefore, is the rate-determining step of the stepwise mechanism. For the stepwise process, the half-life of the G*:C* tautomer is calculated to be approximately 20 fs. The equilibrium constant (K) of the concerted mechanism is of the same order of magnitude as the stepwise mechanism (10⁻⁹). However, the reverse rate coefficient (k′_r) is roughly two orders of magnitude smaller (approx. 10¹² s⁻¹) and, consequently, the half-life of the concerted G*:C* tautomer is calculated to be approximately 3 ps.

Previous gas phase base pair double proton transfer studies have typically observed only one of the G:C mechanisms shown in table 2 for a given QM method [10,12,13,17–19,30]. In doing so, the forward rate coefficient for the G:C double proton transfer has been estimated to be between ${10}^2\hspace{0.17em}\mathrm{and}\hspace{0.17em}{10}^6\hspace{0.17em}{\text{s}}^{-1}$, depending on the approximations made and the QM method used [14–16]. Because of the nature of these QM-only models, the associated errors in measuring the rate coefficient are not reported. A recent nuclear magnetic resonance (NMR) experimental study has estimated the lower bound for the rates of double proton transfer within the (G:T*$\hspace{0.17em}\rightleftharpoons \hspace{0.17em}$G*:T) enol tautomerism to be approximately 10⁵ s⁻¹ or larger [33]. This G:T* base pair mismatch corresponds to guanine connected by three hydrogen bonds to a thymine enol tautomer and is comparatively similar to the G:C base pair. Similarly, the mean forward rate coefficients obtained from our ensemble QM/MM model are 10⁶ s⁻¹ via the first step in the stepwise pathway (k_a) and 10⁵ s⁻¹ via the concerted pathway (k_r). Our results are statistically robust as well as consistent with these experimental data, thus providing substantial improvements on the previous QM calculations.

5.2. A:T tautomerism

For the case of proton transfer occurring in A:T, the greatest percentage of QM/MM replicas (28%) display the concerted single proton transfer reaction. This concerted process involves the transfer of the proton in thymine to adenine to form the A⁺:T⁻ zwitterion (as shown in figure 1b). The preponderance of single proton transfer over double proton transfer is in agreement with a previous study, which demonstrated that the inclusion of explicit water molecules stabilizes the A⁺:T⁻ zwitterion over the A*:T* tautomer [42].

It should be noted that two of the QM/MM replicas performed have shown that the A:T → A*:T* double proton transfer reaction pathway can occur, although locating the exact transition states for these replicas was not possible. This is because the electronic energy levels of the approximate transition state and product (A*:T*) were within close proximity to each other, which resulted in an essentially flat reaction coordinate. Therefore, the transition state geometry optimization was unable to achieve satisfactory gradient convergence for the A*:T* tautomerism and the rate coefficients of the reaction were not determined. The A:T → A*:T* reaction coordinates are displayed in electronic supplementary material, S4. These two replicas were therefore removed from the QM/MM ensemble and replaced by another two, in accordance with the method described in §4.3.

The energetics, kinetics and equilibrium constants for the concerted single proton transfer (A:T → A⁺:T⁻) are presented in the right-hand column of table 3. The mean values and their standard deviations are calculated from a total of seven QM/MM replicas. The reaction coordinate for the concerted single proton transfer in A:T is shown in figure 9. The relative energies of the single proton transfer reaction in A:T (ΔE_rxn, ΔG_rxn) and the transition state ($\mathrm{\Delta }{E}_{}^{‡}$, $\mathrm{\Delta }{G}_{}^{‡}$) are within the standard deviation of one another, similar to that of G:C. The calculated values of the relative transition state energy ($\mathrm{\Delta }{E}_{}^{‡}$) are unevenly distributed between 5.00 and 11.75 kcal mol⁻¹, resulting in a large standard deviation (approx. 2 kcal mol⁻¹). The Gibbs free energy modifies the single proton transfer reaction coordinate shape (figure 9), by stabilizing the transition state approximately 3 kcal mol⁻¹ more than the A⁺:T⁻ zwitterion product (ΔG_rxn). Several studies have calculated the reverse Gibbs free energy barrier to be 0.5 and 3 kcal mol⁻¹ when using the M05-2X/6-311++G** and the M06-2X/6-311++G** methods, respectively, and in conjunction with the PCM solvent approximation [16,74]. Our work is partially in agreement with an earlier study by Céron-Carrasco et al. [30], who used the BP86/6-311++G** method and a micro-hydrated A:T model to predict the reverse free energy barrier to be approximately −1 kcal mol⁻¹. However, the same study predicts the forward Gibbs free energy barrier to be approximately 2 kcal mol⁻¹ smaller than ours (4.4 kcal mol⁻¹). Since there is a distinct lack of QM/MM proton transfer models of A:T in the literature, we expect the difference in our results to be a consequence of the stacking and hydration effects we include. Overall, we find a negative reverse Gibbs free energy barrier for the single proton transfer reaction, which indicates that the A⁺:T⁻ zwitterion product is energetically unstable relative to the transition state (A:T)^‡.

• Download figure
• Open in new tab
• Download PowerPoint

The distribution of the forward (k_r) and reverse (k_r′) rate coefficients for single proton transfer in the A:T base pair is displayed in figure 10. The distribution of k_r (yellow) has a much larger spread than k_r′ (red). Because of this, the per cent standard deviation for the forward rate coefficient is much larger than that for the reverse reaction. The mean concerted single proton transfer rate coefficients k_r and k′_r are within the order of 10¹⁰ and 10¹⁴ s⁻¹, respectively. The reverse rate coefficient is four orders of magnitude larger than the forward one, resulting in an equilibrium constant (K) for the single proton transfer reaction in A:T of the order of 10⁻⁴. From the reverse rate coefficient, the half-life of the A⁺:T⁻ zwitterion is estimated to be 6.3 ± 4.0 fs.

• Download figure
• Open in new tab
• Download PowerPoint

5.3. Relation between proton transfer and single point mutation rates

An upper bound on the spontaneous mutation rates in human DNA can be estimated under the following assumptions.

If the above assumptions are satisfied, the proton transfer equilibrium constant (K) is proportional to the number of tautomers at equilibrium and, in turn, the maximum rate of spontaneous mutations per genome. Double proton transfer is assessed independently to single proton transfer, since the two mechanisms and their effect on spontaneous mutation rates in vivo are mutually exclusive.

The influence of the G:C double proton transfer reaction on spontaneous mutation rates is detailed in the left and centre columns in table 4. Despite occurring less frequently, the concerted pathway produces 33% more G*:C* tautomers at equilibrium than the stepwise pathway. It is estimated that there are a total of 20 G*:C* tautomers (N_taut) within the human genome at equilibrium. It should be noted that the stepwise G*:C* tautomer product has a very short associated half-life (approx. 20 fs; table 3), as well as, on average, a barrierless reverse reaction (−1.0 ± 0.3 kcal mol⁻¹ and $-0.8\pm 0.6\hspace{0.17em}\text{kcal}\hspace{0.17em}{\text{mol}}^{-1}$ for the respective first and second steps). Therefore, in accordance with assumption (4), the stepwise double proton transfer process in G:C is excluded from our single point mutation rate estimate. Given that the concerted tautomerism of G:C has a mean reverse Gibbs energy barrier of 1.0 ± 1.5 kcal mol⁻¹, we are not able to say with certainty whether or not the reverse pathway is actually larger than 0 kcal mol⁻¹, and thus there is no guarantee that the criterion of assumption (4) is met. Also, keeping in mind that the barrier heights in our study are slightly underestimated, as discussed above, and given that the 3 kcal mol⁻¹ threshold in assumption (5) is only an approximation, we are not able to say with certainty whether or not assumption (5) is met. Therefore, we can neither assert that the probabilities of the mutations are null nor affirm that the concerted pathway leads to permanent mutations, especially since the time scale of DNA opening during replication is around one billion times larger [86] than the half-life of the G*:C* tautomer (approx. 3 ps). Therefore, it is reasonable to conclude that only a small fraction of the already small number of G*:C* tautomers (as a consequence of the concerted double proton transfer reaction) may lead to a permanent G:C → A:T mutation during DNA replication; in practical terms, it is truly negligible at less than one base pair per human genome replication.

By contrast, the formation of the double proton transfer A*:T* tautomer is rare and was not observed once within the entire QM/MM ensemble. On the other hand, the A:T → A⁺:T⁻ process was found to be significantly more thermodynamically and kinetically favourable than that of G:C → G*:C*. Specifically, there are on average five orders of magnitude more A⁺:T⁻ zwitterions than G*:C* tautomers at equilibrium, which contributes towards approximately $0.06\text{\%}$ of the A:T content in the human genome. Despite this abundance of A⁺:T⁻ content, the zwitterions have a very short half-life (approx. 6 fs), as well as an average barrierless reverse reaction. Therefore, according to assumption (4), the zwitterion is not considered to influence mutation rates during DNA replication.

That G:C base pairs are more likely to spontaneously mutate is consistent with the universal G:C → A:T transition mutation bias observed in vivo during the replication process [39,87]. Despite the debate concerning the influence of external agents on this mutation bias, recent studies have reported that spontaneous mutations could be the driving force [88–90]. This bias is also observed in humans, whereby the mutation rate per genome replication for G:C sites is two orders of magnitude faster than for non-G:C sites [38,39]. Similar to recent studies by Brovarets & Hovorun [12,13], our work suggests that the double proton transfer tautomerism reactions in DNA are extremely unlikely to contribute towards the overall rates of transition mutations during the DNA replication process.