Sub-picosecond proton tunnelling in deformed DNA hydrogen bonds under an asymmetric double-oscillator model

We present a model of proton tunnelling across DNA hydrogen bonds, compute the characteristic tunnelling time (CTT) from donor to acceptor and discuss its biological implications. The model is a double oscillator characterised by three geometry parameters describing planar deformations of the H bond, and a symmetry parameter representing the energy ratio between ground states in the individual oscillators. If the symmetry parameter takes its maximum value of 1, then we recover a known model which produced CTTs too large to be biologically relevant; but this is reduced by up to 40 orders of magnitude as the symmetry parameter is decreased. We discover that unless the symmetry parameter is close to 1 or 0, the proton's CTT under any planar deformation is guaranteed to be below one picosecond, which is a biologically relevant time-scale. This supports theories of links between proton tunnelling and biological processes such as spontaneous mutation.


Introduction
In the DNA double helix, the two strands of nucleobases are held together by hydrogen bonds, each consisting of a proton being covalently bonded with a donor atom from a donor molecule, and electrostatically attracted to an acceptor atom from an acceptor molecule [1,2]. Löwdin proposed that the proton in an H bond may break away from the donor atom and form a new covalent bond with the acceptor atom, by the mechanism of quantum tunnelling across the potential barrier between the donor and acceptor, and that this process may cause spontaneous mutation [3]. McFadden and Al-Khalili later demonstrated that quantum coherence between the tunnelling proton and its environment can be maintained for biological time-scales, which validates modelling the proton's dynamics as being entirely quantum mechanical [4].
In a normal H bond, all atoms in the donor and acceptor molecules are co-planar, and the donor and acceptor atoms are co-linear with the proton. A planar deformation of the normal H bond is some combination of translations and rotations in the donor-acceptor molecular plane [5][6][7]. It has been theorised that planar deformations of the H bond can have significant effects on the characteristic time-scale of proton tunnelling, and Krasilnikov studied these effects by modelling the potential in the H bond as a double harmonic oscillator which, when the bond is normal, is symmetric about the potential barrier [8]. It was found that the characteristic tunnelling time (CTT) of the proton was extremely sensitive to bond deformation, taking values up to O(10 27 )s, which was not a biologically relevant time-scale.
We propose a generalisation to Krasilnikov's model, in which we associate the symmetry of the double-well potential in a normal H bond with a parameter, γ, whose value equals the energy ratio between a ground-state proton covalently bonded with the donor and one covalently bonded with the acceptor. When γ takes its maximum value of 1, we recover Krasilnikov's model; when 0 < γ < 1, the two local wells in the H bond potential are not equivalent, and the proton has a preferred equilibrium state near the donor rather than acceptor. We further encode the planar deformation of the H bond in three other parameters, d x , d y representing relative shifts between the donor and acceptor, and θ representing the relative rotation, all of which are defined in detail in Section 2. We then derive an analytical expression for the proton's CTT. Fixing all other parameters such as proton mass and covalent bond lengths at values appropriate to DNA H bonds, the CTT is a function of γ, d x , d y and θ. We discover that moderate values of γ guarantee sub-picosecond proton tunnelling, regardless of bond deformation. In Section 3, we discuss the biological implications of our results.

Model and Results
In this Section, we firstly describe the geometry of an H bond under planar deformation, then define our double-well potential within this H bond, before solving the Schrödinger equation under this potential to obtain the proton's wavefunction. From this wavefunction, we derive the proton's CTT. We make the following assumptions and approximations in our model. Firstly, we consider only stationary bonds, meaning that the bond is not actively undergoing deformation whilst proton dynamics is taking place. Secondly, we assume that the lengths and relative angles of all covalent bonds in the donor and acceptor molecules are unaffected by the deformation. In other words, we only consider translations and rotations of the donor molecule as a whole and, independently, of the acceptor molecule as a whole. Finally, even though the proton's global equilibrium is in a covalent bond with the donor atom, we assume that the proton can exist with a higher energy in a locally-stable state of being covalently bonded to the acceptor atom. That there are two local potential minima for the proton in the H bond is the foundation of our double-oscillator model. Since the H bond is planar, it suffices to model the potential for the proton as a function of two spatial dimensions. The geometry of the deformed H bond is shown in Figure 1. Thick lines marked A and D represent, respectively, the acceptor and donor molecules in a deformed bond, whilst the dotted line D marks where the donor molecule would be in a normal bond. N 1 and N 2 mark the acceptor and donor atoms, respectively, and N 2 marks where the donor atom would be in a normal bond. We set up three Cartesian coordinate systems as follows.  Figure 1 as L, D 0 , d x , d y , θ, and defined as follows. L is the distance between N 1 and O 1 , which we assume to be the same as the distance between N 2 and O 2 , as well as the distance bweteen N 2 and O 2 , since we have assumed that no deformation affects the lengths of covalent bonds. D 0 is the distance between O 1 and O 2 , in a normal bond. d x and d y are, respectively, the shifts in the x 1 and y 1 directions of the donor molecule from its normal position, so that, for instance, d x < 0 represents a shift of the donor molecule towards the acceptor molecule. Finally, θ is the anticlockwise angle by which the donor molecule is rotated from its normal orientation, about the point N 2 . We emphasise that the shifts are independent from the rotation, which means that the order in which d x , d y and θ act on the system does not affect its final configuration.
By comparing the coordinates of an arbitrary point in the three systems, O 1 x 1 y 1 , O 2 x 2 y 2 and Oxy, we write down the following coordinate transformation equations.
where θ 1 is the anticlockwise angle from x 1 to x, θ 2 is the anticlockwise angle from x 2 to x, D θ is the distance between O 1 and O 2 in the deformed bond, and λD θ where 0 < λ < 1 is the distance between O 1 and O in the deformed bond. We express θ 1 , θ 2 , D θ and λ in terms of L, D 0 , d x , d y and θ as follows.
which imply cos θ 2 = cos θ 1 cos θ + sin θ 1 sin θ, sin θ 2 = sin θ 1 cos θ − cos θ 1 sin θ, and λ is dependent upon the form of the potential function over the (x, y) plane. For our asymmetric double-oscillator model, we consider a potential function and the proton wavefunction, Ψ(x, y, t), evolves in time according to the Schrödinger equation, In eqs. (4) and (5), m is proton mass, ω 1 and ω 2 respectively are natural angular frequencies of the single oscillators U 1 and U 2 , and g > 0 is an isotropy parameter which we assume to be the same for U 1 and U 2 . We define the symmetry parameter, so that if γ < 1 then there is a lower ground state in U 2 than in U 1 , and this represents the fact that the proton's preferred equilibrium is in U 2 . V is a double oscillator which is identical to U 1 to the left of the line x = 0 and identical to U 2 to the right of x = 0. Thus, there is a potential barrier along the line x = 0 where, in general, we have U 1 = U 2 , so that there is a discontinuity in V .
With the potential function in place, we now calculate λ. In the O 1 x 1 y 1 frame, the local potential well's equipotential curve through the point O is an ellipse, with equation where U 0 is the potential energy at O. One could write a similar ellipse equation, in terms of (x 2 , y 2 ), for the equipotential curve through O in U 2 . Instead, using eqs. (1) and (2), we write both ellipse equations in the Oxy frame, as follows.
Since the ellipses intersect at O, we set (x, y) = (0, 0) in eqs. (7a) and (7b), to obtain from which it follows that We proceed to compute the characteristic time-scale of proton tunnelling from being localised in U 2 to being maximally localised in U 1 , using the Rayleigh-Ritz ansatz [9], in which the ground state wavefunction of the proton is approximately where α 1,2 are complex coefficients, and φ 1,2 are normalised ground state wavefunctions that the proton would have if it existed in the single-well potential U 1 or U 2 , with their domains extended to the infinite plane. We note that if a proton were in the single oscillator U 1 or U 2 , then its ground state energy would be so that the symmetry parameter, γ, equals the energy ratio E 2 /E 1 . Scaling length by we have φ 1 and φ 2 in the following dimensionless forms, in terms of coordinates ξ 1,2 := x 1,2 /x 0 and η 1,2 := y 1,2 /x 0 .
where u 1,2 and φ 1,2 are expressed in terms of coordinates (ξ, η) as follows. Defining and using the dimensionless version of eq. (1), we obtain, for j = 1, 2, where We note that λ [cf. eq. (9)] can now be written from which it follows that r 1 = r 2 . We therefore define For φ j with j = 1, 2, we have, for −∞ < ξ, η < ∞, where Defining the inner product f |g := ∞ −∞ dξ ∞ −∞ dη f * g, we take the inner product of eq. (14) with φ 1 | and φ 2 | respectively to obtain where the overdot denotes differentiation with respect to τ , and We note that since φ 1 , φ 2 are positve, square normalised functions, and since φ 1 ≡ φ 2 , we have 0 < S < 1. Next, using eq. (17), we deduce where By invoking the change of variable ξ → −ξ where necessary, we write, for j = 1, 2 and k = 1, 2, where a = a 2 − a 1 , A jk = A j + A k , and analogous definitions hold for b, B jk , c, C jk , p, P jk , q, Q jk and R jk . Each transition integral I jk can be evaluated exactly, as can the overlap integral, S. We present closed-form expressions for these integrals in the Appendix.
In order for our model to represent tunnelling, rather than scattering, we must have the height U 0 [cf. eq. (8)] of the saddle point in the double-well potential surface being greater than the ground-state energy of φ 2 [cf. eq. (11)]; that is, we must have  Figure 2 shows d crit x as a function of γ. As γ decreases towards 0, greater values of d x would be needed in order to guarantee that every combination of (d y , θ) produces a valid tunnelling model. This is because γ is positively correlated with the steepness of the local potential well U 2 (x 2 , y 2 ). The smaller γ is, the further away from (x 2 , y 2 ) = (0, 0) one needs to go before U 2 reaches the required height, namely the ground-state energy of φ 2 ; thus, in order to ensure that the saddle point between U 1 and U 2 is sufficiently high, U 1 and U 2 must be far enough apart, hence the large d crit x . Meanwhile, as γ → 1, we observe that d crit x → −0.44Å. For 0.01 ≤ γ ≤ 1, we vary d x , d y , θ as follows. −0.45Å ≤ d x , d y ≤ 0.45Å, −90 • ≤ θ ≤ 90 • , and we only consider combinations of (γ, d x , d y , θ) such that eq. (42) holds. We find that for each γ, t p falls in a range between some t min p (γ) and some t max p (γ), and in Figure 3 we present these extremal values as functions of γ. Crucially, our results show that for 0.01 ≤ γ ≤ 0.99, we always have 8.5fs ≤ t p (γ, d x , d y , θ) ≤ 770fs. We also observe that t max p (γ) increases steeply both as γ → 0 and as γ → 1. Indeed, when γ = 1, t max p (γ) becomes ∼ O(10 27 )s; and even though t min p (γ) is still ∼ O(10 −14 )s, t p increases rapidly as (d x , d y , θ) moves away from the combination which minimises t p . Moreover, t min p (γ) is slowly varying with γ, and there is a range of values of γ, namely 0.2 γ 0.4, for which t min p (γ) becomes close to t max p (γ). In this case, varying (d x , d y , θ) has little effect on t p , which contrasts strongly with the large-γ and small-γ cases where t p is very sensitive to (d x , d y , θ). We have defined t ave p (γ) as the mean t p , given a fixed γ, over all combinations of (d x , d y , θ) which satisfy eq. (42), and we have presented t ave p (γ) for 0.1 ≤ γ ≤ 0.6 in the small box in Figure 3. As γ → 1, we have t ave p (γ) ∼ t max p (γ), and for intermediate values of γ, namely γ ≈ 0.3, we have t ave p (γ) ∼ t min p (γ), but as γ → 0, t ave p (γ) is asymptotic to neither t min p (γ) nor t max p (γ).  Furthermore, our results show that for every (γ, d x ), we have This is because a deformation consisting of a shift of d y and rotation of θ is intrinsically identical to one consisting of a shift and rotation of the same magnitudes but both in the opposite direction. Figure 4 shows variations in t p as θ varies between −90 • and 90 • , whilst (γ, d x , d y ) are fixed at certain values. For every combination of (γ, d x ), we have presented only results relating to d y ≥ 0, since one can simply reflect these curves about θ = 0 to obtain results for d y < 0. For fixed (γ, d x ) with d y = 0, the graph of t p (θ) is symmetric about θ = 0, where the graph has a local mimimum under some (γ, d x ) and a local maximum under others; we find from our results that for every γ there is one value of d x at which the graph transitions from having a local minimum to having a local maximum at θ = 0, and that this value of d x increases with γ. For fixed (γ, d x ) with d y = 0, the symmetry of t p (θ) about θ = 0 is broken, and as d y increases, the local extremum which was at θ = 0 when d y = 0 moves towards larger θ. There are cases where this local extremum ceases to exist when d y becomes large, for instance the case of (γ, d x ) = (0.55, −0.3Å), as we can see in Figure 4a: there is a local minimum at θ = 0 if d y = 0 and at θ = 5 • if d y = 0.2Å, but if d y = 0.4Å then this local mimimum disappears. For any fixed (γ, d x , d y ), we always have t p tending to some value as θ tends to ±90 • , typically with several local extrema between θ = 0 and θ = ±90 • ; the value of this limit at ±90 • is dependent only on γ. Calling this limit t 90 p (γ), we have t 90 p (0.55) = 17.1fs, and t 90 p (0.85) = 51.2fs. As γ → 1 and as γ → 0, we have t 90 p (γ) ∼ t max p (γ), and for 0.02 γ 0.4, we have t 90 p (γ) ∼ t min p (γ).  : t p as surfaces over the parameter subspace (d x , d y ), given various combinations of (γ, θ). In each case, the range of d x is d crit We further observe by comparing Figures 4a and 4b that, when γ = 0.85, there is a larger overall variation in t p as a result of varying (d x , d y , θ), compared to when γ = 0.55. This agrees with our observation about Figure 3 that the gap between t min p (γ) and t max p (γ) increases as γ → 1. Indeed, this gap also increases as γ → 0. Moreover, for fixed γ, the larger d x is, the less t p varies with θ or with d y . As we see in Figures 5b, 5c, 5e and 5f, if γ is far from 0, then for fixed (γ, θ), t p as a surface over (d x , d y ) is almost constant given sufficiently large d x . As d x → ∞, t p always tends to some limit, whose value is independent of d y . Meanwhile, we see in Figures 5a to 5c that if θ = 0, then for fixed (γ, θ), t p as a surface over (d x , d y ) is symmetric about the line d y = 0. This is due to eq. (43). If θ = 0 and γ is moderate, such as 0.55, then for each d y sufficiently to 0 we have some small value of d x which maximises t p , as we can see in Figure 5b. This shows that increasing d x , which represents moving the donor away from the acceptor in the H bond, does not necessarily prolong the proton tunnelling. If θ = 0, then the symmetry about d y = 0 is broken, and reflecting a surface for θ > 0 about the line d y = 0 produces corresponding results for θ < 0.

Discussions and Conclusions
We have studied the quantum mechanical tunnelling of a proton across the potential barrier between the donor and acceptor of a planar hydrogen bond in DNA, and computed an analytical expression for the proton's characteristic tunnelling time (CTT) as a function of four parameters describing the geometry of the bond. Three of these parameters, d x , d y and θ, represent the deformation of the H bond from its normal alignment, under the assumption that any deformation consists of planar translations and rotations of the donor and acceptor molecules as independent units. With the acceptor molecule treated without loss of generality as fixed, d x and d y respectively represent the longitudinal and lateral displacements of the donor molecule from its normal position, while θ represents the rotation of the donor molecule about the donor atom from its normal orientation. The fourth parameter, γ, taking values 0 < γ ≤ 1, represents the intrinsic symmetry that the potential in the H bond possesses when the bond is in its normal alignment. When γ = 1, we recover a model previously studied in [8], whose potential function in the normal H bond was symmetric about the potential barrier, so that the local potential wells near the donor and acceptor are equivalent to each other. This symmetry is broken only if some of (d x , d y , θ) is non-zero. For 0 < γ < 1, the symmetry is broken even if d x = d y = θ = 0, in the sense that the local potential well near the donor has a less energetic ground state than the one near the acceptor, and this gives a better representation of the physical property of the H bond than γ = 1. In addition, setting any of d x , d y and θ to non-zero values further distorts the symmetry between the two local potential wells.
The most important difference that generalising from γ = 1 to 0 < γ ≤ 1 has made is that, for most γ values in 0 < γ < 1, the proton CTT is sub-picosecond regardless of (d x , d y , θ). Compared to the γ = 1 case in which some (d x , d y , θ) give CTTs of O(10 27 )s, the sub-picosecond time-scale is much more biologically relevant. Moreover, if γ is such that the CTT is guaranteed to be sub-picosecond, then it varies by no more than 2 orders of magnitude as the H bond deforms. This means that the tunnelling process is much more stable with respect to bond deformation compared to the γ = 1 case, under which the CTT varies by over 30 orders of magnitude as the H bond deforms. Overall, our model under moderate γ-values produces CTTs on a biological time-scale with strong stability against bond deformation, and therefore it supports the theory that proton tunnelling across DNA hydrogen bonds may be a mechanism responsible for biological processes such as spontaneous mutation.
The author is grateful to Dr. Emma Coutts and Dr. Bernard Piette for their kind support.
with K 0,jk = A jk B jk − C 2 jk , K 1,jk = and erfc being the cumulative error function, defined for all real X by erfc(X) = (2/ √ π) ∞ X e −z 2 dz.
The parameters a, b, c, p, q, A jk , B jk , C jk , P jk , Q jk , R jk were defined in the main text.