By treating DNA as a vibrating nonlinear lattice, an activated kinetic theory for DNA melting is developed to capture the breakage of the hydrogen bonds and subsequent softening of torsional and bending vibration modes. With a coarse-grained lattice model, we identify a key bending mode with GHz frequency that replaces the hydrogen vibration modes as the dominant out-of-phase phonon vibration at the transition state. By associating its bending modulus to a universal in-phase bending vibration modulus at equilibrium, we can hence estimate the entropic change in the out-of-phase vibration from near-equilibrium all-atom simulations. This and estimates of torsional and bending entropy changes lead to the first predictive and sequence-dependent theory with good quantitative agreement with experimental data for the activation energy of melting of short DNA molecules without intermediate hairpin structures.

## I. INTRODUCTION

Biological systems are characterized by their metabolism and self-replication. In DNA, self-replication initiates with the local formation of denaturation bubbles by overcoming energy barriers, making denaturation, and hybridization of DNA among the most important processes in biology. These reactions, in response to physical or chemical stimuli, are also the cornerstones of many laboratory techniques such as polymerase chain reaction (PCR) and DNA sequencing, as well as many oligo-based biosensors where a single-stranded DNA binds selectively to its complementary strand.^{1}

Thermodynamics and kinetics are two fundamental aspects of DNA denaturation and hybridization. Through extensive benchmarking against melting temperature measurements, nearest-neighbor models have proven to accurately represent the thermodynamics of these reactions, at least for sequences of up to 24 base pairs (bps).^{2,3} Kinetic properties are harder to capture both through simulations and theories, as the reaction pathway is not well defined^{4} and the time scales of these processes are large enough to prevent extensive studies through all-atom molecular dynamics (MD) simulations.^{5–7}

As most relevant melting and hybridization conditions occur away from the thermodynamic melting temperature and are under non-equilibrium external forces like shear and electric stress,^{8,9} irreversible kinetic rates are often more important than thermodynamic equilibrium constants. Irreversible melting has been shown to produce higher selectivity for base-pairs with comparable thermodynamic equilibrium dissociation constants.^{9} Coarse-grained models have been developed to meet this need for melting and hybridization rates. They have the ability to predict trends in the kinetics but typically produce rates many orders of magnitude faster than experimental values.^{10–12}

Irreversible DNA melting kinetics are known experimentally to obey the Eyring-Polanyi activated rate^{5–7,13–15} $k=kBThexp(\u2212\Delta G\u2021kBT)$, where *k* is the melting rate, *k*_{B} is the Boltzmann constant, *T* is the temperature, *h* is the Planck constant, and $\Delta G\u2021$ is the activation free energy, which can be decomposed into its enthalpy and entropy terms $\Delta G\u2021=\Delta H\u2021\u2212T\Delta S\u2021$. As the separation of the strands during the melting is small enough such that the rate is, within reasonable values, both concentration and salt content independent,^{6,16} the determination of the free energy of activation is sufficient to provide good estimates for the kinetic rate of melting. Nevertheless, we are not aware of any current theory on the activation energy of melting, even from the coarse-grained models.

One major obstacle is the estimate of the activation entropy $\Delta S\u2021$, which is expected to involve cooperative vibration dynamics with both in-phase and out-of-phase phonon fluctuations with respect to the hydrogen bonds at the transition state. Such long-wavelength and slow vibration dynamics are in violation of energy equi-partition as some of the vibration modes will dominate due to the long-range coupling.^{17} Dynamic tracking of this coupling and the evolution of the slow and long-range vibrations are beyond all-atom simulations except for dedicated efforts with large computing machines.^{4,11,18}

There are, however, coarse-grained models that may allow us to capture these cooperative vibrations at the transition state without directly simulating the dynamics. The most extensive all-atom molecular dynamics simulations of the DNA melting process to date present a mechanistic description where the duplex untwists until reaching a conformation which resembles a planar ladder. Then, the bases fray and peel until both strands are held together by a single base pair (see Fig. 1). Finally, this base pair breaks, leading to two fully separated strands.^{4,18} The melting process appears to be quite stochastic, sampling different routes with different fluxes and residence times, and thus making the definition of a simple reaction coordinate quite elusive.^{4} For sequences longer than 10 bps, loops or cruciform secondary structures form during denaturation, further complicating the dynamics and the energy landscape,^{19} with typically lower activation energies.^{20–22} Nevertheless, for a transition state theory, only the activation entropy and enthalpy are needed and the above dissociation of two linear structures should involve certain universal continuum vibration modes that can be quantitatively estimated without dynamic simulation.

In this work, we will model the melting process of short DNA molecules as a thermally activated escape from a multi-dimensional metastable state with specific entropy contributions not present in the usual transition state theory.^{23} The rate can be described instead by the Kramers-Langer theory with a prefactor related to the entropy^{23–26}

For an *N*-dimensional potential, the prefactor Γ_{0} is given by the expression

where *γ* is the damping rate and $\omega 0(i)$ and $\omega b(j)$ characterize the curvatures of the Hamiltonian at the equilibrium and at the barrier (transition state) through which the system escapes. The frequency $\omega b(0)$ is the curvature of the potential at the barrier associated with the negative eigenvalue of an unstable normal mode corresponding to the theoretical reaction coordinate, indicating flow away from the barrier towards the stable states.

The consideration of this prefactor leads to an Eyring-Polanyi equation with an entropy of activation

where assuming the dissociating strands remain linked at the transition state, the thermodynamic entropy of mixing is not included within the activation entropy. Instead, we shall assign these entropies to certain continuum lattice modes with long wavelengths and low frequencies as in phonon melting theories for solids.^{26}

## II. THEORETICAL MODEL

### A. Activation enthalpy

Comparing measured activation enthalpies from Refs. 7 and 13 to thermodynamic enthalpies, we note that the activation enthalpy is comparable to the change in the thermodynamic enthalpy due to all but one base pair,

where Δ*H*^{0} is the thermodynamic enthalpy from the nearest-neighbor approach^{2,3} and *N* is the number of base pairs of the molecule. More specifically, for sequences CAAAAAG and CAGGTCACAG, the measured activation enthalpies are $40.1\xb13$ kcal/mol and $67.2\xb10.4$ kcal/mol, respectively, while the thermodynamic enthalpies Δ*H*^{0} are 47.7 kcal/mol and 73.9 kcal/mol, giving rise to estimated values of $\Delta H\u2021$ of 40.9 kcal/mol and 66.5 kcal/mol from Eq. (4). This agreement suggests that the transition state involves one remaining hinged base-pair as in Fig. 1, and Eq. (4) can be used to accurately estimate $\Delta H\u2021$.

In the absence of other experimental measurements of $\Delta H\u2021$, we shall use Eq. (4) in our theory for the activation enthalpy. We list the thermodynamic enthalpy Δ*H*^{0} and the activation enthalpy $\Delta H\u2021$ estimated from the nearest-neighbor theory for 8 specific sequences in Table I. These are sequences with accurate measured melting rates that will be compared to our kinetic theory. The key to realizing a kinetic theory for DNA melting is then the proper accounting of the entropy change at the barrier.

Sequence . | ω_{bd}/2π (cm^{−1})
. | l (Å)
. | S.D. (Å) . | ΔH^{0} (kcal/mol)
. | $\Delta H\u2021$ (kcal/mol) . | Reference . |
---|---|---|---|---|---|---|

CAAAAAG | 11.995 | 3.19 | 0.11 | 47.7 | 40.9 | 7 |

GGTGAAT | 12.021 | 3.32 | 0.12 | 45.8 | 39.25 | 14 |

CACAGCAC | 11.205 | 3.56 | 0.14 | 59.7 | 52.2 | 42 |

CACGGCTC | 11.229 | 3.37 | 0.10 | 61.1 | 53.5 | 42 |

CAGGAGCA | 11.865 | 3.47 | 0.11 | 56.2 | 49.2 | 43 |

TACGTGGA | 11.186 | 3.44 | 0.09 | 54.7 | 47.9 | 14 |

GGTGAATG | 11.841 | 3.29 | 0.09 | 56.5 | 49.4 | 14 |

CAGGTCACAG | 11.593 | 3.63 | 0.09 | 73.9 | 66.5 | 13 |

GCATCTGGGC | 11.613 | 3.60 | 0.08 | 75.6 | … | … |

Sequence . | ω_{bd}/2π (cm^{−1})
. | l (Å)
. | S.D. (Å) . | ΔH^{0} (kcal/mol)
. | $\Delta H\u2021$ (kcal/mol) . | Reference . |
---|---|---|---|---|---|---|

CAAAAAG | 11.995 | 3.19 | 0.11 | 47.7 | 40.9 | 7 |

GGTGAAT | 12.021 | 3.32 | 0.12 | 45.8 | 39.25 | 14 |

CACAGCAC | 11.205 | 3.56 | 0.14 | 59.7 | 52.2 | 42 |

CACGGCTC | 11.229 | 3.37 | 0.10 | 61.1 | 53.5 | 42 |

CAGGAGCA | 11.865 | 3.47 | 0.11 | 56.2 | 49.2 | 43 |

TACGTGGA | 11.186 | 3.44 | 0.09 | 54.7 | 47.9 | 14 |

GGTGAATG | 11.841 | 3.29 | 0.09 | 56.5 | 49.4 | 14 |

CAGGTCACAG | 11.593 | 3.63 | 0.09 | 73.9 | 66.5 | 13 |

GCATCTGGGC | 11.613 | 3.60 | 0.08 | 75.6 | … | … |

### B. Activation entropy

Estimates of the configurational entropy can be computationally challenging, owing to the difficulty of selecting a proper reaction coordinate and sufficiently sampling the phase space^{27} because of the long length and time scales of the cooperative phonon vibrations at the transition state. We suggest that, through Kramers-Langers formalism, this entropy can be associated with the change in vibration frequencies between the barrier and the double-stranded conformation. Such vibration frequencies can be extracted from the power spectrum or velocity density of states (VDOS) from equilibrium molecular dynamics simulations.^{17,28–30}

Considering a system composed only of harmonic oscillators, the autocorrelation of each oscillator’s elongations is a superposition of all the oscillations within the different eigenfrequencies. Therefore, the Fourier transform of the position autocorrelation functions $\u2211j=1N\u222b\u2212\u221e+\u221e\u27e8rj(t),rj(0)\u27e9e\u2212i\omega tdt$, where **r**_{j} is the displacement of the *j*th oscillator and *j* runs over all atoms of the system, leads to a spectrum with a peak at each eigenfrequency.^{29} The intensities of these peaks are determined by the corresponding oscillator amplitudes, which are inversely proportional to the eigenfrequency *ω* and to the square root of the reduced mass of each oscillator *m*_{j}. Hence, as long as the system is well equilibrated (i.e., the equipartition theorem is fulfilled), a spectrum with the same intensity for each eigenfrequency can be obtained through^{17}

From properties of the Fourier transform, this formula is equivalent to the Fourier transform of the velocity autocorrelation functions

leading to the definition of the power spectrum (or VDOS) of the system. For phonons in an ideal crystal lattice, the peaks of the power spectrum are associated with the different branches of the dispersion relations^{31} (see Appendix A).

In general, low-frequency acoustic vibration modes appear as collective motions of large effective masses and include the motions of many atoms over large distances such that their wavenumber and frequency can approach zero. This makes them quite sensitive to the overall structure of the molecule and thus dominates the entropy change in many DNA transitions.^{32–35} As the amplitude of the normal modes and their consequent importance to the overall motions diminish at higher frequencies, our estimation of the vibrational entropy will focus on low-frequency vibrations.

Figure 2 shows the VDOS calculated directly from equilibrium MD trajectories of a double stranded molecule. Simulation details are presented in Appendix B. By neglecting combination bands, four main branches can be identified from the VDOS of the DNA at equilibrium: the torsional, compressional and longitudinal-bending phonon branches, which correspond to acoustic, in-phase branches, and the transversal branch that corresponds to an optical, out-of-phase branch associated with the breathing of the hydrogen bonds.^{31,36,37} These branches are designated by dispersion relationships from lattice models with various one-dimensional modes that reflect the topology of the DNA. The in-phase acoustic modes represent compressional, torsional, and bending phonons and do not involve direct separation of the two strands, while the out-of-phase optical modes represent the actual separation of the two strands. Nevertheless, both sets of modes will exhibit large changes in vibration frequencies during the melting process, and hence big changes in entropy. The dispersion relationships for the acoustic modes exhibit frequency maxima with respect to the wavenumber, where the group velocities are zero. Wave packet vibrations driven by thermal noise would typically have frequency spectra peaks at values corresponding to zero group velocity.^{31} We hence assigned the most pronounced spectral peaks to the maxima of these dispersion relationships with zero group velocity in Fig. 2. Allocating these dispersion branches to the different stiffnesses of the molecule, it is not only more convenient but also physically reasonable to characterize low-frequency motions through the use of continuous or coarse-grained lattice models.^{35}

### C. Coarse-grained lattice model

In this work, we will model the DNA molecule through a modified version of the Peyrard-Bishop-Dauxois (PBD) model. The PBD model was originally motivated by the study of dynamic properties of DNA such as the breathing of the hydrogen bonds.^{38–41} The potential energy associated with this model can be written as *H* = *H*_{bend} + *H*_{hydr}, where

correspond, respectively, to the elastic bending energy and hydrogen-bonding interactions of the molecule.^{38} The variables *u*_{n} and $vn$ represent the transversal displacements of the complementary nucleotides from their equilibrium position.^{41} Hydrogen bonding interactions are modeled by Morse potentials whose inverse width *a* can be selected to match the vibration frequencies of the hydrogen bonds, while *k*_{b} is a nonlinear parameter related to the bending stiffness of the molecule, which varies through different conformations.

The double lattice structure of the PBD model captures the in-phase and out-of-phase components of the bending dynamics of the two strands. The in-phase “snaking” dynamics do not contribute to the separation of the two strands (melting) but correspond to a distinct phonon branch that can be differentiated from other phonon modes with near-equilibrium simulations. The out-of-phase component does contribute to melting but it is much weaker than the other out-of-phase dynamics [the hydrogen bond vibrations in Eq. (8)] near equilibrium. At the transition state, however, the out-of-phase component becomes dominant and this transition from hydrogen bond dominant out-of-phase vibration at equilibrium to bending dominant out-of-phase vibration at the transition state is a key contribution to the entropy change.

A transformation to the center-of-mass coordinates that represent the in-phase and out-of-phase transversal motions of the two strands ($xn=un+vn2$, $yn=un\u2212vn2$) yields two decoupled equations of motion^{41}

where *m* is the average mass of a base pair (650 amu). Linearizing Eq. (10) around the equilibrium position, the out-of-phase breathing and in-phase bending branches from the dispersion relation can be determined explicitly (see Fig. 2) by

where *l* is the distance between adjacent nucleotides on the same strand, *q* is the wavenumber ($\pi n(N+1)l$ with $n=1,\u2026,N$), and *ω*_{bd} is a frequency related to the bending stiffness of the molecule ($kb/m$).

From our MD simulations of 9 different sequences, at equilibrium, the VDOS spectra show a sequence and length independent parameter *k*_{b} such that $\omega bd/2\pi \u224811.5\xb10.5$ cm^{−1} (see Table I). This frequency will be taken as a universal value in the simplest version of our theory. Similarly, the separation between two nucleotides does not seem to vary far from the universal value of 0.34 nm and this value will be used as a universal value for *l*. This gives rise to a sequence-independent frequency of the first fundamental mode (*n* = 1) for the double stranded bending motions at equilibrium,

that takes values *ω*_{bend}/2*π* of $2.24\u2009\xb1\u20090.1$, $2.00\u2009\xb1\u20090.1$, and $1.64\u2009\xb1\u20090.1$ cm^{−1} for sequences of lengths 7, 8, and 10 bps, respectively. We also obtained a near sequence and length independent vibration frequency *ω*_{H} for the hydrogen bonds such that $\omega H/2\pi \u224892\xb11$ cm^{−1}, in consistency with other experimental and simulation studies.^{44,45}

Despite its success in the description of thermodynamic properties of the molecule, the usual choices of parameters of the PBD model have proven to overestimate denaturation rates by several orders of magnitude.^{10,19,39} In this article, we will modify the PBD model by adding two potentials, *H*_{tors} and *H*_{compr}, corresponding to the torsional and compressional energies of the molecule,^{46,47}

where the variables *θ*_{n,n−1} are related to the angular displacement of two consecutive hydrogen bond vectors (being *θ*_{0} the angle near the equilibrium position), the variables *z*_{n,n−1} are related to the separation of the centers of mass of two consecutive Watson-Crick pairs in the direction of the axis of the molecule (being *z*_{0} the separation at the equilibrium position) and *k*_{θ} and *k*_{c} are nonlinear parameters related to the torsional and compressional stiffnesses of the molecule (see Fig. 3). The addition of these potentials to a similar model by Sicard *et al.*^{46} has proven to be able to reproduce opening and closure times of DNA denaturation bubbles, suggesting this addition may lead to an improvement of the PBD model towards properly describing kinetic properties of the system.

Note that, from these potentials, the torsional and compressional branches can also be determined explicitly through,

where *ω*_{t} and *ω*_{c} are frequencies related to the torsional ($k\theta /IM$) and compressional ($kc/m$) stiffnesses of the molecule, with *I*_{M} being the moment of inertia per bps around the twisting axis.

Two main phenomena have been suggested to give rise to the difference in thermodynamic entropy between the double and single stranded configurations: the drop in the elastic stiffness due to the decrease of the overlap of the *π*-electrons of adjacent bases, which is mostly responsible of the stacking interactions, and the increment in angular configurations that the bases may sample upon the unwinding of the helix.^{38,40,46,48,49} In our kinetic theory with similar elasticity changes, we will introduce the coupling between base-stacking and bending elasticity by varying the bending persistence length, while the change in torsional stiffness will be coupled to the angular sampling space of the bases. As the compressional modulus of the molecule does not change significantly upon melting,^{50} we will consider *k*_{c} as constant. As a result, compressional frequency change does not occur for thermal melting of short DNAs and will hence be omitted in this study. However, for shear or force induced melting, the change in compressional entropy may become important and we will include it here for completeness.

### D. Change in frequencies at the barrier

where $\Delta Sbend\u2021$, $\Delta Stors\u2021$, $\Delta Scompr\u2021$, and $\Delta Soop\u2021$ correspond to the change in bending, torsional, compressional, and out-of-phase breathing vibration frequencies, while $\Delta Src\u2021$ is a correction to the activation entropy due to the change in curvature along the reaction coordinate.

Optical Kerr-effect spectroscopies show that the out-of-phase branch associated with the breathing of the hydrogen bonds disappears upon melting,^{44} while MD studies show that denaturation bubbles have much smaller bending and torsional moduli than the double helix.^{46} As the compressional modulus of the molecule does not change significantly upon melting, we will neglect the entropy due to the change in the compressional vibrations ($\Delta Scompr\u2021\u22480$) and will focus on the change in bending, torsional, and breathing vibrations.

#### 1. Bending entropy

At the barrier, we assume that the molecule is held together by a single base pair. Hence, it is expected that both strands undergo independent bending dynamics and the bending persistence length is then comparable to that of a single stranded DNA. As there are *N* bending modes, each with a frequency proportional to the square root of the persistence length,^{51,52} we obtain an entropy term

where $lpds$ is the persistence length associated with the double stranded configuration and $lpss$ is the persistence length associated with the single stranded one. At these short scales, the persistence length associated with the double stranded molecule can be calculated from the VDOS by using the theory of vibrations of a homogeneous rod,^{51,52} as

where *M* = *mN* is the total mass of the molecule, *L* = (*N* − 1)*l* is the contour length, *ω*_{bend} is the frequency of the first fundamental mode from Eq. (13), and *p*_{1} is an inverse length such that *p*_{1}*L* = 4.73. This treatment of the persistence length yields average values of $9\xb12$, $13\xb13$, and $25\xb12$ nm for duplexes of lengths 7, 8, and 10 bps, respectively [see Fig. 4(a)], which are comparable to previous simulated and experimental studies.^{53,54}

The same treatment applied to calculate the persistence length of single stranded molecules yields unreasonably stiff values, as in reality these molecules form strong non-bonded interactions which make them more folded or collapsed than ideal rods. There are many inconsistencies in experimental and theoretical calculations of the persistence length of single stranded molecules, which fall within a rather wide range of values (0.7–7.8 nm).^{56} Estimates from the end-to-end distance^{53} of our molecular dynamics simulations give rise to results comparable to other simulations that work in similar length scales,^{55} yielding persistence lengths of short single stranded molecules close to $1\xb10.5$ nm [see Fig. 4(b)].

#### 2. Torsional entropy

The torsional stiffness associated with the double stranded molecule can also be calculated from the VDOS by using the theory of vibrations of a homogeneous rod;^{51,52} however, to our knowledge, there is no consensual definition of the twisting stiffness of a single stranded molecule. At the barrier, the torsional modulus can be estimated by measuring the twist angle between consecutive base pairs in the molecule, which has been suggested by Sicard *et al.*^{46} to follow a power-law relation $k\theta \u221d(N\u22121)\u2212\alpha $, where *N* − 1 is the number of melted bases and $\alpha =2.2\xb10.1$. Assuming the change in inertia is also included in this formula, this estimate yields a torsional entropy

#### 3. Out-of-phase entropy

When the DNA molecule is in its double stranded configuration, transverse optical phonons associated with the vibrations of the hydrogen bonds are supported, having each optical phonon a finite frequency at zero wave vector. Upon melting, these transverse waves are no longer supported and the zero wave vector frequency of this phonon has to become zero. Therefore, as seen in Eq. (12), the out-of-phase optical phonon modes associated with the breathing of the hydrogen bonds become out-of-phase acoustic phonon modes in a phenomenon which resembles the soft phonon phase transition observed in crystal structures.^{57,58}

This physical effect can be observed in the PBD model [see Eq. (10)]. Once the hydrogen bonds break, out-of-phase bending effects previously withheld by the much stronger hydrogen bonds come into play, as the equations of motion become

Out-of-phase dynamics and in-phase dynamics share the same bending stiffness *k*_{b}; hence, knowledge of the in-phase bending dynamics can lead to predictions of the out-of-phase vibration frequencies.

We will designate base pair *j* to be the single connecting (hinge) base pair at the transition state, as assumed in our activation enthalpy estimation [Eq. (4)]. The exact value of *j* does not affect our estimate of the out-of-phase bending vibration frequency. The vibration frequencies of the out-of-phase normal modes obey the equation $A\omega bd\u20322v=\omega v$, where *A* is the $N\xd7N$ matrix

whose *j*th column is zero, $\omega bd\u2032$ is the average intrinsic bending frequency at the barrier ($kb/m$ evaluated at the barrier), and the eigenvector **v** corresponds to the normal mode *ω*. As is true for the bending frequencies of Eq. (19), we correct for the change in the bending modulus *k*_{b} when the DNA has transformed from the double-stranded configuration at equilibrium to the configuration at the transition state. This is achieved by approximating the value of *k*_{b} at the barrier through the product of *k*_{b} at the equilibrium position and $lpss/lpds$. The frequency $\omega bd\u2032$ can hence be estimated from the double stranded bending frequency at equilibrium through

leading to values $\omega bd\u2032/2\pi $ of $3.32\u2009\xb1\u20092$, $2.76\u2009\xb1\u20091.7$, and $2.13\xb1\u20091.0$ cm^{−1} for sequences of lengths 7, 8, and 10 bps, respectively.

The matrix *A* is singular because of the zero eigenvalue associated with the zero *j*th column, corresponding to the hinge base pair. However, this zero eigenvalue can be factored out from the determinant of *A* such that the product of all the non-zero eigenvalues, corresponding to the product of the frequencies of the out-of-phase mode, can be estimated explicitly to be $\omega bd\u2032N\u22121$, independent of the position of the hinge. We note that the geometric average of these out-of-phase vibration frequencies $\omega bd\u2032/2\pi $ (about 50–100 GHz) is lower than the hydrogen bond vibration frequencies at equilibrium *ω*_{H}/2*π* (at about 2.75 THz) by almost two orders of magnitude.

The breakage of the hydrogen bond oscillators followed by the manifestation of these out-of-phase bending modes leads to an entropy term

where the out-of-phase branch at the equilibrium position is approximated by a constant *ω*_{H} as $\omega H\u226b\omega bd$. This represents the entropy change due to the softening of the out-of-phase mode in Eq. (10), as the hydrogen bonds at equilibrium melt (except for those at the hinge base pair) and yield to the out-of-phase bending modes.

#### 4. Correction due to curvature effects

Out-of-phase dynamics involve the breathing of the two strands, dominating the melting process. Therefore, by working near the barrier where all except the *j*th base pairs have melted, the curvature along the reaction coordinate can be associated with the curvature along a one-dimensional reaction coordinate which involves the out-of-phase dynamics of this last base pair. Approximating the frequency of the hydrogen bond breathing modes to the average vibration frequency of each hydrogen bond oscillator *ω*_{H} and making a selection of the location of the transition state and its negative curvature, the correction factor can be estimated as

where the frequency at the equilibrium position along this reaction coordinate is taken to be *ω*_{H} and at the barrier will be assumed to take the most negative frequency of the Morse potential, $\omega b0=\u2212\omega H/2$. The damping coefficient will be chosen as *γ*/2*π* = 0.5 ps^{−1} in concordance with experimental results.^{44} This rather arbitrary designation of the position of the transition state along the reaction coordinate does not introduce much error, as this curvature effect will be shown later to contribute to less than 2% of the activation entropy.

In summary, our theory yields the following formula for the activation entropy:

where near-equilibrium all-atom simulations have provided universal values for the hydrogen bond vibration frequencies $\omega H/2\pi \u2009\u2248\u200992\u2009\xb1\u20091$ cm^{−1} and the characteristic bending frequency $\omega bd/2\pi \u2009\u2248\u200911.5\u2009\xb1\u20090.5$ cm^{−1}. The damping factor is roughly *γ*/2*π* = 0.5 ps^{−1}, leading to a correction term $T\Delta Src\xb1$ of barely 0.5 *k*_{B}*T*. For our sequences, the persistence length ratio between double stranded and single-stranded DNA is estimated by

where the average base pair mass *m* = 650 amu and nucleotide separation *l* = 0.34 nm from near-equilibrium simulations in Table I are used. The values of the persistence length ratios obtained are $lpds/lpss=12.7\xb17$, $17.3\xb110.0$, and $33.3\xb117$ for sequences of lengths 7, 8, and 10, respectively. The softening of the bending modulus increases the entropy change of both the out-of-phase breathing vibration and the bending mode. However, the former is further affected by the change in frequency due to the breakdown of the hydrogen bonds to give rise to the bending out-of-phase modes. The entropy is hence only a function of the base number *N*. With the sequence dependent activation enthalpy of Eq. (4), we can hence estimate the activated melting rate of short DNAs.

## III. COMPARISON TO MEASURED RATE DATA AND DISCUSSION

Figure 5(a) shows the comparison of experimentally measured free energies of melting for eight different sequences to our theoretical estimates (see Appendix C for details on the extraction of the enthalpies and entropies from measured rates and the estimate of their error bars). Our theoretical estimates were calculated by taking the average persistence length ratios. An analysis of how the choice of different parameters changes these results is presented in Appendix D. By using the estimated thermodynamic enthalpy of each of these sequences and using our previous estimates for the entropy, we were able to predict the experimental free energies of activation of the eight sequences at different temperatures within 10% (2 kcal/mol) without fitting, which is comparable to the experimental error of the measurements^{7,14,42} [see Fig. 5(a)]. The experimental kinetic rates vary over 4 orders of magnitude, as seen in the supplementary material.

Figure 5(b) shows the contribution of each entropy term from Eq. (18) to the activation entropy. We find that the torsional entropy change from Eq. (21), usually neglected in the PBD model, is just as important as the bending entropy change from Eq. (19). We also find that, using changes in stiffness comparable to those suggested by Sicard *et al.*,^{46} our theoretical calculations yield free energy barriers of similar orders to those obtained for denaturation bubbles through metadynamics calculations.

For the sequences analyzed, our approach introduces the sequence dependence through the nearest-neighbor enthalpy terms. For example, CAGGAGCA and CACAGCAC possess equal percentage of GC base pairs and activation entropies, but very different activation free energies. Further analysis should be performed to guarantee the sequence independence of the entropy; however, each entropy term can be estimated explicitly from equilibrium simulations.

This approach represents the first kinetic theory for the activation free energy of DNA melting, which captures the important elastic torsional, bending, and acoustic vibration entropies. We are currently extending this theory to longer molecules with secondary intermediated structures, whose dynamics can also be captured by a local version of the current theory. Other than such topological changes, which will obviously affect the all-important vibration spectra and the activation entropy, we expect protein and intercalator interaction will also likewise affect the vibration spectra and significantly change the melting kinetics and thermodynamics. For longer DNAs under shear stress, the change in compressional entropy may become important. However, the compression modulus *k*_{c} may be a strong function of the secondary structure and we shall include nonlinear corrections to our Hamiltonian *H*_{compr}, as was done for the hydrogen bond Hamiltonian *H*_{hydr}.

## SUPPLEMENTARY MATERIAL

See supplementary material for information about how the experimental activation energies are obtained from the kinetic measurements and how their error bars estimated.

## ACKNOWLEDGMENTS

We thank J. K. Whitmer and S. A. Corcelli for useful discussions. S.S. and H.C.C. are partially supported by NIHR21CA206904. Z.P. acknowledges the support from the National Science Foundation (No. 1706436-CBET) and the Indiana Clinical and Translational Sciences Institute (No. UL1TR001108) from the National Institutes of Health, National Center for Advancing Translational Sciences, and Clinical and Translational Sciences Award.

### APPENDIX A: RELATION BETWEEN POWER SPECTRUM AND DISPERSION BRANCHES

where $aj(\omega )=\u222b\u2212\u221e+\u221er\u0307j(t)ei\omega tdt$. The displacement of the *k*th coordinate of the atom *j* can be decomposed by means of the normal modes λ through

where *ω*_{λ} is the eigenfrequency of the mode λ and *A*_{λkj} corresponds to the coordinates in the base of the eigenvectors defined by these modes. Considering the derivative of the previous equation and using the definition of the delta function, we may write

where *n* runs over all atoms (and coordinates) of the system and the masses have been incorporated into the amplitudes along the normal modes. If the equipartition theorem holds, then all normal modes have the same kinetic energy and $|A\lambda m|\omega \lambda 2=kBT$. Hence, as the density of states is given by Ref. 61 $\rho (\omega )=\u2211\lambda \delta (\omega +\omega \lambda )$, the power spectrum can be written as *F*(*ω*) = 3*Nk*_{B}*Tρ*(*ω*).

For phonons in a crystal lattice, the density of states presents singularities where the group velocity is zero^{62} (i.e., at the maximum of each dispersion branch). Hence, the maxima of the dispersion branches are expected to correlate with the highest peaks of the power spectrum.

### APPENDIX B: SIMULATION DETAILS

DNA structures were generated using the NAB module of Amber 14.^{63} Each structure was solvated in a rectangular box with explicit TIP3P water molecules^{64} (a minimum distance of 10 Å between box border and any atom of the solute), neutralized with sodium ions by means of the Amber 14 leap module, and modeled using the ff99+parmbsc0 force field for DNA.^{65} Long range electrostatic interactions were calculated with the Particle Mesh Ewald (PME) method^{66} and a real space cutoff radius of 9 $\xc5$.

Simulations were performed using the SANDER module of Amber 14. For minimization, we ran 500 steps of steepest descent followed by 500 steps of conjugate gradient minimization, with a 500 kcal/mol restraint force on the DNA molecule, and 1000 steps of steepest descent followed by 1500 steps of conjugate gradient minimization with no restraints.

During 20 ps, the system was heated up at a constant volume from 0 K to a temperature of 275 K while the DNA molecule was weakly restrained with 25 kcal/mol. Langevin temperature controls, SHAKE constraints on hydrogen atoms, and a 2 fs time step were applied. Subsequently, positional restraints were gradually removed during 150 ps at 275 K and constant pressure (1 bar). The system was then equilibrated for 50 ps at NVT and later for 1 ns at NVE.

During the rest of the simulation, temperature was kept constant using the Langevin thermostat with a damping factor of 0.5 ps^{−1}. Pressure was maintained at 1 bar using the Nose-Hoover Langevin piston. A distance cutoff of 11.0 $\xc5$ was applied to short-range, non-bonded interactions, and 9.0 $\xc5$ for the smothering functions.

In order to obtain more accurate results for the contour length of the molecules, 100 ns of equilibration was performed followed by 10 ns of production. The output of the production stage was generated every 100 fs. The cpptraj module of AMBER 14 was used for basic analysis^{67} (centering and imaging of the trajectories, computations of distances and angles, etc.). The specden plugin of the VMD Signal Processing Plugin Package was used to calculate the VDOS of the DNA molecule.^{68} Convergence was tested by comparing the VDOS at 5 and 10 ns.

Properties of the DNA molecule were determined after removing the solvent molecules and center of mass motions from the trajectory. The contour length was obtained by calculating the distance between consecutive centers of mass of the base pairs (excluding the bases to avoid overestimations due to potentially flipped out ones). VMD was used for visualization of different conformations and trajectories.^{68} Figure 1 was obtained by using an enhanced melting simulation method that will be published elsewhere.

### APPENDIX C: EXPERIMENTAL ACTIVATION ENERGIES

All experimental data considered have been previously published (see Table I for references) and the activation enthalpies and entropies for the dissociation of the DNA duplex were derived from the temperature dependence of the denaturation rates *k* by using the Eyring-Polanyi equation $k\u2009=\u2009kBThexp(\u2212\Delta G\u2021kBT)$, where $\Delta G\u2021\u2009\u2009=\u2009\u2009\Delta H\u2021\u2212T\Delta S\u2021$. The logarithmic form of this equation yields

Given the minimal temperature dependence of $ln(kBTh)$, a plot of $ln(k)$ versus 1/*T* yields $\Delta H\u2021$ and $\Delta S\u2021$ from the slope and intercept, respectively.^{14}

The experimental methods reported can be classified in three categories depending on the implemented technique: stopped flow analysis,^{42,43} temperature jump analysis,^{7,13} and single-channel current analysis on a protein nanopore.^{14} The two last techniques provide the denaturation rates explicitly while the stopped flow technique only provides the hybridization rate and estimates the denaturation rate by assuming detailed balance,

where *k*_{on} is the hybridization rate and *K*_{D} is the equilibrium dissociation constant, which is experimentally derived by using a van’t Hoff method on the basis of the dependence of the melting temperature on the concentration of DNA strands.^{14,42}

In supplementary material, we present the Eyring-Polanyi plots of various reported kinetic data with error bars. When errors on the activation parameters or rates were not given explicitly, these errors led to the error bars in Fig. 5(a).

### APPENDIX D: ROBUSTNESS OF THE THEORETICAL PREDICTIONS

Although our theory provides good estimates to the free energies of denaturation through using the suggested values, it is pertinent to give an idea of how the choice of these parameters could change the results. This will be studied by altering the parameters *ω*_{H}, *ω*_{bd}, and *η* within the suggested domains where these values are physically reasonable ($\omega H/2\pi \u2208[91,93]$ cm^{−1}, $\omega bd/2\pi \u2208[11.0,12.0]$ cm^{−1}, and $\eta 2\u2208[5,19],[7.3,27.3]$, and [16.3, 50.3] for sequences of lengths 7, 8, and 10 bps, respectively). These parameters are only related to the activation entropy, and hence we will focus on this term.

Defining $\Delta S0\u2021$ as the activation entropy evaluated with our suggested set of parameters and $\Delta S1\u2021$ as the activation entropy evaluated with a different set, from Eq. (28) we may write

where the numeric subindexes identify the set of parameters. To simplify the discussion, let us define a function $\Delta Sdiff\u2021=1N(\Delta S1\u2021\u2212\Delta S0\u2021)$, which can be approximated by

Multiplying these terms by the temperature *T*, we can see that the second and third terms add up to less than 0.055 *k*_{B}*T* (leading to net differences of less than 0.55 *k*_{B}*T* for sequences of 10 bps) and the biggest differences arise from the quotient of the persistence lengths.

Defining the parameter $\mu =\eta 12/\eta 02$, this parameter varies within [0.4, 1.5] for *N* = 7, [0.4, 2.2] for *N* = 8, and [0.5, 1.5] for *N* = 10, what leads to changes in activation entropies of up to 0.9 *k*_{B}*T* per base pair (net differences of up to 9 *k*_{B}*T* for sequences of 10 bps). This number is quite large mostly due to the error bars of the persistence length of the single stranded molecule ($1\xb10.5$ nm). If we only allow up to a 20% deviation on this persistence length,^{55} then the difference in entropy is close to 0.45 *k*_{B}*T* per base pair, leading to net deviations from our suggested free energies close to 3.15 *k*_{B}*T*, 3.6 *k*_{B}*T*, and 4.5 *k*_{B}*T* for sequences of lengths 7, 8, and 10 bps, respectively. These variations are all comparable to the 2 kcal/mol error we estimated for the experimental error and for our theory in Fig. 5(a) ($1\u2009kBT\u22480.5922$ kcal/mol at 298.15 K).