Modelling proton tunnelling in the adenine–thymine base pair

Abstract

The energies of the canonical (standard, amino-keto) and tautomeric (non-standard, imino-enol) charge-neutral forms of the adenine–thymine base pair (A–T and A*–T*, respectively) are calculated using density functional theory. The reaction pathway is then computed using a transition state search to provide the asymmetric double-well potential minima along with the barrier height and shape, which are combined to create the potential energy surface using a polynomial fit. The influence of quantum tunnelling on proton transfer within a base pair H-bond (modelled as the DFT deduced double-well potential) is then investigated by solving the time-dependent master equation for the density matrix. The effect on a quantum system by its surrounding water molecules is explored via the inclusion of a dissipative Lindblad term in the master equation, in which the environment is modelled as a heat bath of harmonic oscillators. It is found that quantum tunnelling, due to transitions to higher energy eigenstates with significant amplitudes in the shallow (tautomeric) side of the potential, is unlikely to be a significant mechanism for the creation of adenine–thymine tautomers within DNA, with thermally assisted coupling of the environment only able to boost the tunnelling probability to a maximum of 2 × 10⁻⁹. This is barely increased for different choices of the starting wave function or when the geometry of the potential energy surface is varied.

---

1 Introduction

In their seminal work,¹ Watson and Crick proposed that the genetic code is stored in the form of hydrogen-bonds between the canonical purine and pyrimidine nucleic acid bases: adenine–thymine (A–T) and cytosine–guanine (C–G), which form the chains of the double helix structure of the DNA molecule. They recognised that tautomerization alters the hydrogen-bonding patterns and therefore could lead to mismatches in the canonical base pairs. The rare tautomer hypothesis of spontaneous mutagenesis thus states that mutations can arise through the spontaneous formation of high energy tautomers. However, accurate quantitative studies of the probability of proton transfer along the hydrogen bonds between the bases have only been possible in recent years once detailed knowledge of the potential energy surface felt by the protons could be accurately mapped. It is generally accepted that this surface is in the shape of an asymmetric double-well potential, allowing for the proton to be transferred either via classical over-the-barrier hopping or via quantum tunnelling through the barrier.

The likelihood of quantum tunnelling occurring within DNA and the question of whether this is a significant contributor to spontaneous point mutations has been the subject of much research over the past few decades.^2–15 It has become more prominent in recent years as more sophisticated computational methods for investigating the problem have become feasible. In addition, it is now known that the external environment, such as the presence of water molecules, also plays a role in the stability and structure of DNA. The question of interest in this study is the extent to which the environment might play a role in promoting or inhibiting proton transfer in the H-bonds in A–T and C–G base pairs.

Many studies have focussed on the hydrogen transfer itself without considering whether tunnelling is the cause. Some have claimed that tunnelling in DNA^7,16,17 is either not possible or so unlikely as to be statistically negligible. Others claim it is a reasonable possibility.¹⁸ In the case of DNA base pairs, Löwdin² correctly declared in 1963 that once a DNA replication event had occurred, the protons in the connecting hydrogen bonds would be in one of several quantum states, some or all of which could lead to potential tunnelling events, affecting any future replication events.

The canonical forms of the DNA base pairs are very stable, as would be required of the molecules responsible for the storage of the genetic information of life. However, one early study by Parker and Van Everv⁴ calculated the energy levels of both standard DNA base units (A–T and G–C) and a non-standard DNA base unit (G–T). They estimate that if the proton is initially in the ground state of the double well potential, then tunnelling across from the deeper well to the shallow well is extremely unlikely due to the large asymmetry. By using a simple semi-classical WKB approximation they conclude that the only way to promote tunnelling is through incoming radiation that can excite the proton to higher energy levels.

This conclusion was supported and expanded upon by Villani in a series of papers between 2005 and 2012, investigating whether the tautomeric Watson–Crick configurations are stable enough to exist for reasonable lengths of time.^18–21 He identifies that mapping a potential energy surface is a suitable approach to model hydrogen atom transfer, and that a double hydrogen transfer that preserves electro-neutrality, as is the case with A–T and the tautomeric A*–T*, requires a double well with relatively large asymmetry. Florián et al. had earlier studied the stability of the canonical and tautomeric structures in 1994⁷ and concluded that the difference in energy between the potential energy barrier and the A*–T* configuration is very small and thus A*–T* reverts to A–T very quickly.

In his 2007 paper, Villani shows using a two-dimensional model that the probability of the A*–T* tautomer existing at any given time is around 4 × 10⁻⁴, although the probability of other tautomers that exhibit charge separation occur with a slightly higher probability.¹⁹ In the four-dimensional model (with the four degrees of freedom being the two lengths of the hydrogen bonds, and the two out-of-plane positions of the hydrogen atoms in the hydrogen bridges), the probability is closer to 3 × 10⁻³, but this does not change his main conclusion that A*–T* does not appear often enough to be considered “a mutation point in the DNA chain”.

One of the most comprehensive studies of the suitability of various theoretical frameworks for solving the adenine–thymine system was published by Bende in 2009.²² Various types of DFT calculations using different exchange–correlation functionals were compared to results from ab initio MP2 (second-order Møller–Plesset perturbation theory),^23,24 which is more precise but more computationally expensive than DFT. An important conclusion of this work is that DFT produces very good results for this type of system when using certain functionals such as BLYP,^25,26 B3LYP,²⁷ and BHLYP,²⁸ with the latter two performing the best. Lower-level functionals such as PBE^29,30 and KMLYP³¹ did not perform so well, particularly in the case of the hydrogen bond lengths and interaction energies. The choice of exchange–correlation functional is very important when using DFT to find energies and optimise molecular geometries.

A further illuminating study was carried out by Brovarets and Hovorun in 2013,¹⁷ which also compares MP2 and various types of DFT calculation, but in addition considers tunnelling and its potential impact on the appearance of adenine–thymine tautomers within DNA strands. They showed that the tautomerisation of A–T to A*–T* via double hydrogen transfer is both concerted and asynchronous; that is, the reaction pathway (found using the Hessian-based predictor-corrector integration algorithm^32–34) does not feature an intermediate step, but the hydrogen atoms move one at a time.

Curiously, they demonstrate a step-wise double hydrogen transfer where the hydrogen atom on the N–N bridge moves across first, followed by that of the O–N bridge, counter to what Villani had found three years earlier.²⁰ It is possible that there are various reaction pathways and that different models lead to different pathways being more prominent. Crucially, however, Brovarets and Hovorun also considered the effect of the surrounding environment: when modelling the hydrophobic interfaces of protein–nucleic acid interactions using a low dielectric constant (ε = 4), the stability and lifetime of A*–T* was not improved significantly enough for A*–T* to occur very often. This is counter to various studies^35–38 that put the instance rate of A*–T* tautomers at between 10⁻³ and 10⁻⁶. Hence, they conclude that tautomerisation cannot be a source of spontaneous mutations during the DNA replication process, in agreement with Florián et al.⁷ In addition, they claim that the appearance of A*–T* tautomers does not occur via proton tunnelling but instead by above-barrier vibrations.

In 2010, Pérez et al. studied the stability of the Watson–Crick base pair tautomers and found that, although previous studies based on classical mechanics had declared them stable, when taking into account quantum effects, they were found to be dynamically metastable (i.e. stable enough to have reasonable lifetimes).³⁹ In that work, several versions of density functional theory were used, as well as methods such as MP2 and Hartree–Fock. They showed that Hartree–Fock performs poorly for DNA base pairs and, in agreement with Bende,²² DFT using hybrid functionals such as PBE0 and B3LYP produce good results.

2 The H-bonds in A–T

A double hydrogen bond exists between A and T, with C being joined to G via a triple hydrogen bond. This makes the A–T bond somewhat easier to model (Fig. 1(a)). When DNA undergoes replication, the weak hydrogen bonds between the bases of the nucleotides break and are free to recombine with new nucleotides in the nucleoplasm, assisted by the enzyme DNA polymerase.⁴⁰ However, if the canonical forms of the nucleotides become tautomeric due to a proton tunnelling event in the hydrogen bond, the base pair is said to consist of new nucleotides, known as A* and T* (Fig. 1(b)). If a DNA strand splits, A* will no longer combine with T and it is far more energetically favourable for it to combine with C instead. Likewise, T* will almost certainly combine with G. If this error makes it past the error correction processes such that neither new DNA strand matches the original, a mutation is said to have occurred.

Fig. 1 Double hydrogen transfer in an adenine–thymine (A–T) pair.

Such a hydrogen bond can be modelled as a superposition of quantum states inside a double well potential. To get from a canonical to a tautomeric state a proton transfer would have to take place. Various studies attempt to describe this process as either a double proton transfer^2,5,7,11 or a single proton transfer,^6,8,10 in order to explain these mutations. Once in a tautomeric state, DNA replication of this base pair could result in a codon error (see Fig. 2).

Fig. 2 Adapted from Löwdin.² Canonical A–T base pair undergoing tautomerisation and then replication. Neither new base pair is of the usual A–T or C–G forms.

In this work, we employ density functional theory to model the 1-D potential energy surface felt by one of the two protons in the double H-bond of the A–T base pair. A full two-dimensional map of the reaction pathway would be ideal but also computationally prohibitive; fortunately, a one-dimensional map appears sufficient for such a proton transfer.³⁹ Once this static double-well potential is found it is used to solve the time-dependent generalised Liouville equation, including a dissipative term to model the action of the environment (e.g. H₂O molecules) surrounding the base pair. We assume that if one of the two protons moves across then so will the second one, whether together or sequentially in a two-step process, to create the tautomeric form, A*–T*. It has been found that using ab initio constrained molecular dynamics calculations⁴¹ that double potential transfer in the A–T base pair is a stepwise and most likely asynchronous process. However, whether the double proton transfer occurs in a concerted or step-wise fashion is of little importance when utilising a one-dimensional reaction pathway because the “secondary” hydrogen atom will always move to preserve electroneutrality. Thus, the movement of a single proton is the rate-determining action that should be investigated.^20,17,39,42

3 Mapping the potential energy surface

The differences between the various density functional theory (DFT) methods hinge on their choice of exchange–correlation functional. In the absence of a perfect method for calculating geometries and energies, the most suitable functional in terms of performance and computational complexity must be chosen. For this work, we have used B3LYP,²⁷ a hybrid generalised gradient approximation (GGA) functional derived from the GGA functional BLYP, which itself combines Becke's exchange functional²⁵ with the correlation functional by Lee, Yang, and Parr.²⁶ There is evidence in the literature to suggest that hybrid GGA functionals like B3LYP are well suited to a system like adenine–thymine^{19–22,39,42,43} (and indeed match well with results from higher level methods such as MP2). It should be noted, however, that B3LYP has been demonstrated to underestimate hydrogen bond strengths for A–T.⁴⁴ In addition, DFT functionals in general tend to underestimate barrier energies compared to other methods.^45–47

The structure of the canonical form of A–T is listed in the S22 database,⁴⁸ which was calculated using MP2 with the cc-pVTZ basis set.⁴⁹ The structure data was converted to a format readable by CASTEP⁵⁰ and then a geometry optimisation was performed using the “precise” preset for plane-wave energy cut off, with BLYP norm-conserving pseudopotentials.⁵¹ Since the tautomeric form, A*–T*, is not available in the S22 database, it was necessary to estimate its form before a geometry optimisation could be performed. To do this, the B3LYP optimised geometry for A–T was adjusted to create a starting geometry for A*–T*. The lengths and the angles of the hydrogen bond ions were estimated based on the average of values reported in literature, then those hydrogen bonds were fixed in place whilst the rest of the geometry was optimised. Lastly, that result was used as the starting structure for a final pass during which all ions were free to move. Both the canonical and tautomeric optimised geometries are shown in Fig. 3 and the difference in energy between the two forms, ΔE = E_A*–T* − E_A–T was found to be 4613.5 cm⁻¹ (0.572 eV or 13.19 kcal mol⁻¹). This is compared to existing results in literature in Table 1.

Fig. 3 Diagrams of optimised geometries generated using the software package J mol,^52,53 with bond lengths displayed.

Table 1 Summary of results for A–T base pair energies. ΔE is the energy difference between the canonical and tautomeric states; E_B is the barrier height compared to the canonical state. See Fig. 4

Source	ΔE (cm^−1)	E B (cm^−1)
This work	4597	8066
Villani (2005)^18	4517	—
Bende (2009)^22	4759–5243	—
Yue-Jie Ai et al. (2010)^43	4597	4759–4920
Florián et al. (1994)^7	3307	3388
	5243–6614	6614–7017
Guallar et al. (1999)^54	7872	7864
Gorb et al. (2004)^16	4194–5727	4033–6049
Cerón-Carrasco et al. (2009)^55	3780–5070	—
Pérez et al. (2010)^39	2984–3468	3307–4355
Gobbo et al. (2012)^42	5001	5650–18 600
Herrera and Torro-Labbe (2004)^56	4200–7000	15 700–17 100

---

It is important to note that the z-axis positions for all of the ions are essentially static for both states (unsurprising given that the molecule is very “flat”). The relative change in y-axis positions for the protons in the hydrogen bonds is ∼5% that of the change in x-axis positions. This suggests that a one-dimensional potential, varying along the x-axis, can accurately describe the transition state between these two molecular forms, in agreement with other authors.^18,39

Once minimum energies and optimised geometries were known for stable versions of both A–T and A*–T*, it was necessary to find the height and shape of the barrier between them to complete the double well potential picture. There are a few different methods to perform a transition state search, in order to calculate a minimum energy pathway or potential energy surface of a geometry transition. The most powerful one offered by CASTEP is the complete LST/QST,⁵⁷ which combines the linear synchronous transit (LST) and quadratic synchronous transit (QST) methods. Initially, a single LST interpolation is performed, bracketing the maximum energy between the reactants and product, which is followed by continual energy minimisations in “directions conjugate to the reaction pathway” until a true energy minimum is reached. This approximation is then used to perform QST “constrained minimisations” to much more accurately refine the transition state. Once this is complete the cycle is repeated beginning with another LST maximisation until a stationary point is found, which represents the transition state barrier.

Next, a fourth-order polynomial fit was carried out to map the full shape of the potential with five boundary conditions: the two well minima energies, the barrier maximum, and the two well minima ζ-axis positions, where ζ is the dimensionless reduced position co-ordinate of the proton along the x-axis defined as ζ = x/a, where a is half the distance between the two well minima positioned such that they are equidistant from x = 0, placing them at ζ = ±1. (Note that, due to the asymmetry in the potential, the peak of the barrier may not be exactly at ζ = 0.) The nature of a fourth order polynomial is similar to that of a quartic well and produces the high well walls necessary for outer-box boundary conditions. The resulting polynomial is described in eqn (1) and appears alongside the energies found during the CASTEP transition state search in Fig. 4. Note that the 8065.73 factor is to convert eV to cm⁻¹.

V(ζ) = 8065.73(−4464.49 + 0.429ζ − 1.126ζ² − 0.143ζ³ + 0.563ζ⁴) (1)

Fig. 4 Polynomial fit applied to the results of CASTEP's transition state search for stable and optimised forms of A–T & A*–T*.

In Fig. 4 the potential energy scale has been adjusted such that the energy of the canonical state E_A–T = 0 since only the relative energies of the barrier and minima, and the well shape, are important.

4 Time evolution of the density matrix

The density matrix operator describing a single proton in a 1-D double-well potential is defined aswhere the time-dependent state vector for the system, |ψ(t)〉, is expanded in the well's energy eigenstate basis

Ĥ₀|ϕ_i〉 = E_i|ϕ_i〉, (3)

whereand m is the mass of the proton.

In order to model an open (dissipative) system, coupling to the environment (in the limit of weak coupling to a Markovian bath) can be included in the generalised Liouville equation that include a dissipative (Lindblad) term on the right hand side^58,59

where the extra term is generally written in the form:^58–62whereand W_ij are environment induced transition rates between well states |ϕ_i〉 and 〈ϕ_j|.

Substituting the above back into eqn (5) leads to diagonal and off-diagonal density matrix elements

There are a number of ways of calculating the transition matrix elements, W_ij. They are derived here using the microscopic theory of Meyer and Ernst.⁶³ We begin with the full system + bath Hamiltonian,

Ĥ = Ĥ₀ + Ĥ_b + ΔĤ, (9)

where Ĥ₀ is the internal Hamiltonian of the A–T system defined in eqn (3) and Ĥ_b is the bath Hamiltonian defined as a sum of harmonic oscillatorswhere m is the set of bath oscillators, p_m are their momenta, q_m are their spatial positions and ω_m their frequency. Finally, is the interaction between the system and bath, with coupling constant, f_m.

The transition probability W_ij between states i and j is defined as

where ω_ij is a transition frequency depending on the energy of the eigenstates i and j, The correlation functions, C_ij(t), required in eqn (11) are calculated from an appropriately chosen power spectral density function of the active bath displacement:where T is temperature, ω_p is a characteristic phonon frequency of the heat bath and ΔV_R is the rearrangement energy gained by the bath oscillators upon displacement from Q_m = 0 to their optimal values for configurations ζ = ±1 of the tunnelling system near its potential minima.⁶³ In summary, this standard definition of the power spectral density function^62–64 is related to the chosen model for the bath oscillators, using Debye theory.

The power spectral density function can be transformed to give two correlation functions, one for the active bath displacement

and another for the active bath momentum (the force of the heat bath on the quantum system) A final correlation function, C^tun, for the active bath rearrangement, is given by

C^tun(t; d²) = e^{d2ℏ−2(C_ss(t)−C_ss(0))}, (17)

where d is the transfer distance between any two energy eigenstates (or doublets):

d = |ζ_ii′ − ζ_jj′|. (18)

Here ζ_ij′ is obtained through a transformation of the ζ matrix, which comprises elements The diagonal elements of this matrix are simply the position expectation values of the energy eigenfunctions of the quantum system, 〈ϕ_i|ζ|ϕ_i〉. The transformation used to obtain the ζ′ matrix is defined to be that which diagonalises ζ for pairs of states that are considered tunnel doublets. These involve states of nearly-equal local vibrational excitation in each of the two wells, where the only direct method to shift from one to the other is to tunnel through the barrier.

Finally, the correlation function

C_ij(t) = (ζ_ij′)²C_rr(t)C^tun(t; d²) (20)

controls the rate of transition between states i and j and can be entered into eqn (11) to find transition probabilities for each eigenstate pair. Note that these correlation functions are explicit for baths consisting of simple harmonic oscillators within this model; they can be thought of as defining how the bath oscillators behave when interacting with the quantum system. However, the power spectral density function given by eqn (14), J_rr, can be thought of as describing the fundamental properties of the heat bath and these are typically altered to fit experimental data. Changes to the power spectral density would affect the system-bath interaction but the formula used here is based on the general form typically found in textbooks.^62–64

Unfortunately, transition probabilities are not trivial to calculate for this system because the difference between each energy eigenvalue is relatively large, and as the difference in energy between eigenstates i and j increases, the smaller the values for C_rr(t) and C_ss(t) become around t = 0. For numerical integration, many of these smaller values must be summed, and they are less prominent compared with the already small values at more extreme values of t. Fortunately, an analytical approach for calculating these correlation functions using the residue theorem can be used to solve the troublesome integrals, C_ss(t) and C_rr(t).

The analytic solutions of C_ss(t) and C_rr(t) are shown below:

where , c = ℏ/k_BT, , and γ_n = 2πn. With analytic solutions for C_rr(t) and C_ss(t), the expression for C_ij(t) also becomes analytic. However, there is a single case, at C_rr(t = 0), where the infinite sum is divergent, leading to numerical difficulty when solving the Fourier transform integral over time in eqn (11). However, there is a way around this. The infinite sum in eqn (21) and (22) only contributes to the real part of C_rr(t). With this in mind, it is useful to find a means of solving the FT integral for W_ij without requiring the real component of C_ij(t).

To this end, W_ij is split into its real and imaginary constituents as follows:

Since W_ij is a parameter, this can be simplified towhere C_ij = C^Re_ij + iC^Im_ij and both C^Re_ij and C^Im_ij are real functions of t.

Eqn (26) still includes the real component of C_ij(t) and so doesn't solve the divergent-sum problem entirely. However, there are several useful identities and relations that can be employed to solve this. Since ω_ij = ℏ⁻¹(E_i − E_j) = −ω_ij, then C_ij(t) = C_ij(t). Also, the transition probabilities between energy eigenstates must conform to the principle of detailed balance,⁶³

due to the fact that the Hamiltonian for the system-bath interaction is not Hermitian but has Hermitian-like symmetry, i.e. (ΔH_kj)^† = ΔH_ij.⁶³ If we rewrite W_ji using the above relations, we can see that it is nearly identical to W_ij,with the only difference between eqn (28) and (26) being that the two terms in the integral are subtracted instead of added. This is helpful because by substituting for W_ij and W_ji into eqn (27) we can rearrange to write W_ij in terms of the known quantity, β_ij, This means we do not need to know the real part of C_ij(t) and therefore the troublesome real part of C_rr(t)). The above integral is easy to solve numerically and converges with reasonable limits (t = ±1 ps).

5 Results

The energy eigenstates for the adenine–thymine double potential well, shown in Fig. 4, were calculated using the DSTEV function included in the LAPACK library.⁶⁵ The first eight eigenstates are shown in Fig. 5.

Fig. 5 Energy eigenstates of adenine–thymine base pair 1-D model.

The most important feature to note here is that the first five eigenstates are wholly localised in the deeper well on the left (i.e. the canonical state). The 6th eigenstate is the first that is localised in the shallower well (the tautomeric state). The 5th and 6th eigenstates have sufficiently close energy eigenvalues to be considered a tunnel doublet. The same is true for the 7th and 8th eigenstates. The 12th energy eigenstate (not shown in figure), at 7407 cm⁻¹ (or 0.92 eV), is the first one above the barrier, which means there are six eigenstates below the barrier with energies high enough to allow for tunnelling between the wells.

With the eigenstates known, the time evolution of the density matrix can be computed using eqn (5). The proton is chosen to be initially in the ground state and the initial wave function is thus (|ψ(t = 0)〉 = |ϕ₁〉). This is a reasonable choice since the ground state probability amplitude in such an asymmetric double well is almost entirely in the deep well anyway, where we would wish to start the proton in. Of course it should be noted here that without a dissipative Lindblad term in eqn (5) there would be no tunnelling since this is a stationary state. The effect of the dissipative term is therefore to allow transitions between eigenstates and it is only the coupling to the higher states (6th to 11th) that leads to any significant probability of finding the proton in the shallow well.

The values of several parameters in the Lindblad term are now required. The temperature of the environment, appearing in the spectral density function, was fixed at 320 K. The other two parameters, ΔV_R and ω_p, relate to the properties of the heat bath of harmonic oscillators. However, without experimental data for this system it is difficult to choose appropriate values and it is instructive to carry out computations over a range of values. Fortunately, since the double hydrogen-bond in the A–T system shares its properties with that of a much simpler and better studied molecule, the benzoic acid dimer,^66–68 it was decided to use parameters from the latter (ΔV_R = 44 cm⁻¹, ℏω_p = 186 cm⁻¹)⁶⁹ as a benchmark around which to explore.

The time evolution was computed over a range of 10 ns and repeated for ΔV_R = 10, 20, 44, 100, 400 cm⁻¹ and ℏω_p = 20, 100, 500 cm⁻¹. The probability of tunnelling is plotted over time in Fig. 6. Each sub-figure is for a specific value of ω_p with curves for different values of ΔV_R. The probability of the wave function being in the shallow well at each time interval required for outputted results was calculated from the diagonal elements of the density matrix in the position basis, |ψ(ζ)|², then normalising this wave function and integrating over the relevant range along the ζ-axis that covers the width of the shallow well. We feel that this is a more informative quantity to calculate in our approach than the equivalent tunnelling times (for example using a semi-classical WKB approach) or transition rates (using Fermi's golden rule).

Fig. 6 Tunnelling probability in adenine–thymine nucleobase pair for |ψ(t = 0)〉 = |ϕ₁〉, T = 320 K, and varying Lindblad parameters. ΔV_R has units cm⁻¹.

The key feature of these plots is that the probability of finding the proton in the shallow well increases over time (very slowly in comparison with rate of tunnelling itself – several nanoseconds compared typically with around a tenth of a picosecond), but this probability reaches a plateau of around 1.6 × 10⁻⁹ for all parameters used; in fact, the only parameter that affects the maximum tunnelling probability is temperature. We showed in a previous study⁶⁸ that for the case of the benzoic acid dimer the tunnelling rate is reduced at lower temperatures since the weaker coupling to the environment means it is less likely that there will be strong transitions between eigenstates, particularly those at energies above the bottom of the shallow well. Secondly, both ω_p and ΔV_R affect how quickly the system reaches equilibrium.

Next we investigated the sensitivity of the tunnelling probability to the A–T well geometry. Firstly, the height of the barrier above to the deeper well, E_B, was reduced by 30%. This adjusted potential and its energy eigenstates are shown in Fig. 7. The first 7 eigenstates are now wholly localised in the deeper well, with the 7th & 8th and 9th & 10th states forming tunnel doublets. The 8th eigenstate is the lowest in the shallow well, but even this state is partially in the deep well because the barrier is so close in energy to the shallow well minimum. The tunnelling probability in this potential are shown in Fig. 8.

Fig. 7 Energy eigenstates of adenine–thymine base pair 1-D model with E_B reduced by 30%.

Fig. 8 Tunnelling probability in adenine–thymine nucleobase pair for |ψ(t = 0)〉 = |ϕ₁〉, T = 320 K, and varying Lindblad parameters. The height of the potential barrier compared to the deeper well minimum has been reduced by 30%. ΔV_R has units cm⁻¹. This figure should be compared with Fig. 6(a).

There are two important differences to note here. Firstly, the time it takes for an equilibrium to be reached is much reduced in all cases, being rarely greater than 100 ps. Secondly, the maximum tunnelling probability has increased slightly to 2.4 × 10⁻⁹.

Next, the original adenine–thymine potential was adjusted by increasing the depth of the shallow well compared to the potential barrier, E_B − ΔE, by 50%. This adjusted potential and its energy eigenstates are shown in Fig. 9. There are now 3 pairs of tunnel doublets (states 4 & 5, 6 & 7, and 8 & 9) and there are more localised eigenstates in the shallow well, starting with the 4th. The tunnelling probability results for this new potential are shown in Fig. 10.

Fig. 9 Energy eigenstates of adenine–thymine base pair 1-D model with E_T–C reduced by 50%.

Fig. 10 Tunnelling probability in adenine–thymine nucleobase pair for |ψ(t = 0)〉 = |ϕ₁〉, T = 320 K, and varying Lindblad parameters. The depth of the shallow well compared to the potential barrier has been increased by 50%. ΔV_R has units cm⁻¹. This figure should be compared with Fig. 6(c).

Now, the time taken to reach equilibrium is dramatically increased compared to the standard adenine–thymine model, being approximately an order of magnitude greater. However, the tunnelling probability is increased significantly to almost 4 × 10⁻⁷. This still represents a very low chance of tunnelling having taken place, but is nearly 250 times higher than for the original adenine–thymine potential.

Finally, it was decided to alter the initial (t = 0) wave function from the ground state to a Boltzmann distribution in order to simulate an “average” starting wave function based on the density of states:⁷⁰

where E_i is the ith energy eigenvalue. Fig. 11 shows the tunnelling probability when the system is set up in this way without any environment interaction present.

Fig. 11 Tunnelling probability in adenine–thymine nucleobase pair when the starting wave function, ψ(t = 0), is set to a Boltzmann distribution with no environment coupling.

Now the tunnelling probability oscillates rapidly due to the complex nature of the wave function. The probability of tunnelling is still extremely low at ∼1.6 × 10⁻⁹, which is due to the fact that the first energy eigenstate (being localised to the deep well) is by far the largest contributor, having a coefficient of 0.9936. Including the effect of the environment makes no difference to the resulting tunnelling probability and is negligible compared to the complex oscillations due to interference between the system's eigenfunctions.

It is worth noting that even when starting the proton in the highest excited state below the bottom of the shallow well (|ψ(t = 0)〉 = |ϕ₅〉), mimicking similar tests carried out by Gobbo et al.,⁴² the transitions between eigenstates due to the environmental coupling bring the tunnelling probability down to 1.6 × 10⁻⁹. This takes approximately 70 ns with default parameters, during which time the tunnelling probability briefly reaches a peak of 8 × 10⁻⁵. Whilst this is much higher than when the proton is in the ground state (which is to be expected since the barrier is relatively lower), it still cannot account for a large percentage of occurring tautomers.

6 Summary and conclusions

An improved method for discovering transition probabilities was devised and applied to the adenine–thymine potential. These were calculated and applied in order to find that the tunnelling probability, even with an environment present, is very low in this model. It thus seems unlikely that the instance rate of A*–T* is due to proton tunnelling in any significant way, with other mechanisms being more important. The effect of changing the asymmetry and barrier height of the potential is small, although the former is the more potent factor of the two. In addition, altering the mass parameter to simulate a deuterated molecule has a minimal effect on the tunnelling probability and thus it seems unlikely that any change in the instance rate of A*–T* due to deuteration is because of tunnelling.

The main purpose of this study has been to add to the body of work on the influence of quantum tunnelling on proton transfer within DNA base pair H-bonds. Although previous studies have looked at the stability, lifetimes, and instance rates of tautomeric Watson–Crick adenine–thymine base pairs within DNA,^{2,7,16–21,35–39,42,71} few have look specifically at tunnelling, particularly in the absence of external radiation exciting the protons to higher energies where tunnelling would be more likely. Even the notion that base pair tautomers exist for any meaningful length of time within DNA is yet to be settled, with some indications that A*–T* is metastable and others that it reverts to the canonical A–T form too quickly to have any influence during DNA replication.

The adenine–thymine potential was created using DFT with the B3LYP functional. The results obtained compared favourably with those in the existing literature, particularly in the case of the energy difference between the canonical state and tautomeric state. An initial wave function was set up in both the ground state (|ψ(t = 0)〉 = |ϕ₁〉) and in a Boltzmann distribution, then transition probabilities were calculated using an improved, mostly analytical, method to include the dissipative effects of coupling to an external heat bath. The density matrix was then allowed to evolve in time and it was found that tunnelling was promoted by the environmental coupling, but by a very low amount (2 × 10⁻⁹). It was confirmed in an earlier work⁶⁸ that the effect of the dissipative Lindblad term is to induce thermally assisted tunnelling, but (crucially) not to encourage classical over the barrier hopping. This was determined in that work, which was applied to proton tunnelling in the benzoic acid dimer molecule, by comparing the probability of finding the proton in the shallow well in two eigenstate basis calculations. In the first, a large basis set of energy eigenstates was used, which included many eigenstates with energies above the barrier. In the second, a smaller basis set was used, restricted to just those eigenstates below the barrier, some of which could support tunnelling if their energies were above the base of the shallow well. Hardly any change to the tunnelling probability was found between the two calculations, implying that tunnelling is the dominant mechanism (namely, that the only way for the proton to find its way to the shallow well was via quantum tunnelling rather than classical over the barrier hopping. The reason for this is that the transition probabilities W_ij between eigenstates that link to the above barrier states were negligibly small. This is key because it demonstrates that the primary effect of dissipation is not to allow classical over-the-barrier hoping but to increase the likelihood of tunnelling below the barrier. We have made the assumption that this argument also holds for the case of proton transfer in A–T. However, it may be that further work is necessary to confirm this.

To test the sensitivity of the outcomes to the input parameters, various adjustments were made to the potential and environment parameters. The only environment parameter to actually affect the maximum tunnelling rate was temperature, thus confirming that future adjustments to the other parameters would not drastically affect the conclusions of this work. Despite this, it would be interesting to generate more accurate parameters for the nucleoplasm environment using experimental data and then repeat the procedure. In terms of the potential itself, both the depth of the shallow well and the height of the barrier were adjusted and it was found that the asymmetry was by far the most important factor (even more so than the particle mass) for determining the increase in tunnelling probability in this model. Therefore, the range of barrier heights in the literature, which is fairly large (3200–18 600 cm⁻¹ or 0.4–2.3 eV), is unlikely to be critical; even if the barrier height is dropped significantly (by 30%), the affect on tunnelling rates remains relatively small.

A comparison between the method used here, where first a static potential is generated via a DFT approach and then the time-dependent master equation for an open quantum system is solved, with a fully time-dependent DFT approach would be of interest and could prove useful in testing our conclusions. Finally, even if quantum tunnelling really does not play an important role in point mutations in DNA, there is still the potential for utilising the tools developed here in other biochemical processes, such as those involving proton transfer promoted by enzyme catalysis.

References

1. J. D. Watson and F. H. C. Crick, Nature, 1953, 171, 964–967.
2. P.-O. Löwdin, Rev. Mod. Phys., 1963, 35, 724–732.
3. P.-O. Löwdin, Advances in Quantum Chemistry, Academic Press, 1966, vol. 2, pp. 213–360.
4. B. R. Parker and J. V. Everv, Chem. Phys. Lett., 1971, 8, 94–99.
5. W. Sokalski, H. Romanowski and A. Jaworski, Adv. Mol. Relax. Interact. Processes, 1977, 11, 29–41.
6. A. O. Colson, B. Besler and M. D. Sevilla, J. Phys. Chem., 1992, 96, 9787–9794.
7. J. Florián, V. Hrouda and P. Hobza, J. Am. Chem. Soc., 1994, 116, 1457–1460.
8. M. Hutter and T. Clark, J. Am. Chem. Soc., 1996, 118, 7574–7577.
9. K. Brameld, S. Dasgupta and W. A. Goddard, J. Phys. Chem. B, 1997, 101, 4851–4859.
10. J. Bertran, A. Oliva, L. Rodríguez-Santiago and M. Sodupe, J. Am. Chem. Soc., 1998, 120, 8159–8167.
11. V. L. Golo and Y. S. Volkov, Int. J. Mod. Phys. C, 2003, 14, 133–156.
12. E. K. Ivanova, N. N. Turaeva and B. L. Oksengendler, 2011, arXiv:1105.8282.
13. P. Ball, Nature, 2011, 474, 272–274.
14. D. Jacquemin, J. Zúñiga, A. Requena and J. P. Cerón-Carrasco, Acc. Chem. Res., 2014, 2467–2474.
15. A. Kumar and M. D. Sevilla, Chem. Rev., 2010, 110, 7002–7023.
16. L. Gorb, Y. Podolyan, P. Dziekonski, W. A. Sokalski and J. Leszczynski, J. Am. Chem. Soc., 2004, 126, 10119–10129.
17. O. O. Brovarets and D. M. Hovorun, J. Biomol. Struct. Dyn., 2013, 32, 127–154.
18. G. Villani, Chem. Phys., 2005, 316, 1–8.
19. G. Villani, Chem. Phys., 2007, 336, 143–149.
20. G. Villani, Phys. Chem. Chem. Phys., 2010, 12, 2664–2669.
21. G. Villani, Chem. Phys., 2012, 394, 9–16.
22. A. Bende, Theor. Chem. Acc., 2010, 125, 253–268.
23. C. Møller and M. S. Plesset, Phys. Rev., 1934, 46, 618–622.
24. M. Head-Gordon, J. A. Pople and M. J. Frisch, Chem. Phys. Lett., 1988, 153, 503–506.
25. A. D. Becke, Phys. Rev. A: At., Mol., Opt. Phys., 1988, 38, 3098–3100.
26. C. Lee, W. Yang and R. G. Parr, Phys. Rev. B: Condens. Matter Mater. Phys., 1988, 37, 785–789.
27. A. D. Becke, J. Chem. Phys., 1993, 98, 5648–5652.
28. A. D. Becke, J. Chem. Phys., 1993, 98, 1372–1377.
29. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1996, 77, 3865–3868.
30. J. P. Perdew, K. Burke and M. Ernzerhof, Phys. Rev. Lett., 1997, 78, 1396.
31. J. K. Kang and C. B. Musgrave, J. Chem. Phys., 2001, 115, 11040–11051.
32. H. P. Hratchian and H. B. Schlegel, J. Chem. Phys., 2004, 120, 9918–9924.
33. C. Dykstra, G. Frenking, K. Kim and G. Scuseria, Theory and Applications of Computational Chemistry: The First Forty Years, Elsevier Science, 2011, pp. 195–249.
34. H. P. Hratchian and H. B. Schlegel, J. Chem. Theory Comput., 2005, 1, 61–69.
35. S. Basu, R. Majumdar, G. Das and D. Bhattacharyya, Indian J. Biochem. Biophys., 2005, 42, 378–385.
36. A. Shukla, Chemical Basis of Inheritance, Discovery Publishing House Pvt. Limited, 2009, p. 51.
37. M. Williams and L. Maher, Biophysics of DNA-Protein Interactions: From Single Molecules to Biological Systems, Springer, 2010.
38. T. A. Kunkel and K. Bebenek, Annu. Rev. Biochem., 2000, 69, 497–529.
39. A. Pérez, M. E. Tuckerman, H. P. Hjalmarson and O. A. von Lilienfeld, J. Am. Chem. Soc., 2010, 132, 11510–11515.
40. A. Stanford, Foundations of biophysics, Academic Press, 1975.
41. S. Xiao, L. Wang, Y. Liu, X. Lin and H. Liang, J. Chem. Phys., 2012, 137, 195101.
42. J. P. Gobbo, V. Saur, D. Roca-Sanjun, L. Serrano-Andrs, M. Merchn and A. C. Borin, J. Phys. Chem. B, 2012, 116, 4089–4097.
43. Y.-J. Ai, F. Zhang, G.-L. Cui, Y. Luo and W.-H. Fang, J. Chem. Phys., 2010, 133, 064302.
44. T. van der Wijst, C. F. Guerra, M. Swart and F. M. Bickelhaupt, Chem. Phys. Lett., 2006, 426, 415–421.
45. B. J. Lynch and D. G. Truhlar, J. Phys. Chem. A, 2001, 105, 2936–2941.
46. M. N. Glukhovtsev, R. D. Bach, A. Pross and L. Radom, Chem. Phys. Lett., 1996, 260, 558–564.
47. R. L. Bell, D. L. Taveras, T. N. Truong and J. Simons, Int. J. Quantum Chem., 1997, 63, 861–874.
48. B. Energy and G. Database, Adenine–Thymine Watson–Crick Complex, 2006.
49. P. Jurecka, J. Sponer, J. Cerny and P. Hobza, Phys. Chem. Chem. Phys., 2006, 8, 1985–1993.
50. S. J. Clark, M. D. Segall, C. J. Pickard, P. J. Hasnip, M. J. Probert, K. Refson and M. Payne, Z. Kristall., 2005, 220, 567–570.
51. R. Sabatini, Quantum Espresso - Pseudopotentials, 2014.
52. Jmol, Jmol: an open-source Java viewer for chemical structures in 3D, 2014.
53. A. Herráez, How to Use Jmol to Study and Present Molecular Structures, http://Lulu.com, 2007.
54. V. Guallar, A. Douhal, M. Moreno and J. M. Lluch, J. Phys. Chem. A, 1999, 103, 6251–6256.
55. J. P. Cerón-Carrasco, A. Requena, J. Ziga, C. Michaux, E. A. Perpte and D. Jacquemin, J. Phys. Chem. A, 2009, 113, 10549–10556.
56. B. Herrera and A. Toro-Labbe, J. Chem. Phys., 2004, 121, 7096–7102.
57. N. Govind, M. Petersen, G. Fitzgerald, D. King-Smith and J. Andzelm, Comput. Mater. Sci., 2003, 28, 250–258.
58. V. Gorini, A. Kossakowski and E. C. G. Sudarshan, J. Math. Phys., 1976, 17, 821–825.
59. G. Lindblad, Commun. Math. Phys., 1976, 48, 119–130.
60. C. Scheurer and P. Saalfrank, J. Chem. Phys., 1996, 104, 2869–2882.
61. K. Lendi and R. Alicki, Quantum Dynamical Semigroups and Applications (Lecture Notes in Physics), Society for Industrial and Applied Mathematics, 11th edn, 1987, vol. 286.
62. H. Breuer and F. Petruccione, The Theory of Open Quantum Systems, Oxford University Press, USA, 2002.
63. R. Meyer and R. R. Ernst, J. Chem. Phys., 1990, 93, 5518–5532.
64. P. Rebentrost, M. Mohseni, I. Kassal, S. Lloyd and A. Aspuru-Guzik, New J. Phys., 2009, 11, 033003.
65. E. Anderson, Z. Bai, C. Bischof, S. Blackford, J. Demmel, J. Dongarra, J. Du Croz, A. Greenbaum, S. Hammarling, A. McKenney and D. Sorensen, LAPACK Users' Guide, Society for Industrial and Applied Mathematics, Philadelphia, PA, 3rd edn, 1999.
66. C. Foces-Foces, L. Infantes, F. Aguilar-Parrilla, N. S. Golubev, H.-H. Limbach and J. Elguero, J. Chem. Soc., Perkin Trans. 2, 1996, 349–353.
67. M. J. Wjcik, K. Szczeponek and M. Boczar, Int. J. Mol. Sci., 2003, 4, 422–433.
68. A. D. Godbeer, J. S. Al-Khalili and P. D. Stevenson, Phys. Rev. A: At., Mol., Opt. Phys., 2014, 90, 012102.
69. A. Stockli, B. H. Meier, R. Kreis, R. Meyer and R. R. Ernst, J. Chem. Phys., 1990, 93, 1502–1520.
70. A. Durrant, Quantum Physics of Matter, Taylor & Francis, 2000.
71. E. S. Kryachko, Int. J. Quantum Chem., 2002, 90, 910–923.

---