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.

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 defined4 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 rate5–7,13–15k=kBThexp(ΔGkBT), where k is the melting rate, kB is the Boltzmann constant, T is the temperature, h is the Planck constant, and ΔG is the activation free energy, which can be decomposed into its enthalpy and entropy terms ΔG=ΔHTΔS. 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 ΔS, 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.

FIG. 1.

Conformations of the DNA duplex CAAAAAG sampled during an enhanced melting simulation.

FIG. 1.

Conformations of the DNA duplex CAAAAAG sampled during an enhanced melting simulation.

Close modal

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 entropy23–26 

k=Γ0exp(ΔHkBT).
(1)

For an N-dimensional potential, the prefactor Γ0 is given by the expression

Γ0=ωb(0)ω0(1)ω0(2)ω0(Ñ)2πγωb(2)ωb(3)ωb(Ñ),
(2)

where γ is the damping rate and ω0(i) and ωb(j) characterize the curvatures of the Hamiltonian at the equilibrium and at the barrier (transition state) through which the system escapes. The frequency ω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

ΔS=kBln(hkBTΓ0),
(3)

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 

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,

ΔH(N1N)ΔH0,
(4)

where ΔH0 is the thermodynamic enthalpy from the nearest-neighbor approach2,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±3 kcal/mol and 67.2±0.4 kcal/mol, respectively, while the thermodynamic enthalpies ΔH0 are 47.7 kcal/mol and 73.9 kcal/mol, giving rise to estimated values of ΔH 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 ΔH.

In the absence of other experimental measurements of ΔH, we shall use Eq. (4) in our theory for the activation enthalpy. We list the thermodynamic enthalpy ΔH0 and the activation enthalpy ΔH 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.

TABLE I.

Average intrinsic bending frequencies and length per bp (with standard deviation) from equilibrium MD simulations. Thermodynamic enthalpies (ΔH0) estimated using the nearest-neighbor approach.2,3 Activation enthalpies (ΔH) estimated from Eq. (4). References for measured melting rates of each sequence are listed.

Sequenceωbd/2π (cm−1)l (Å)S.D. (Å)ΔH0 (kcal/mol)ΔH (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. (Å)ΔH0 (kcal/mol)ΔH (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 … … 

Estimates of the configurational entropy can be computationally challenging, owing to the difficulty of selecting a proper reaction coordinate and sufficiently sampling the phase space27 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 j=1N+rj(t),rj(0)eiωtdt, where rj is the displacement of the jth 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 mj. 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 through17 

F(ω)=ω2j=1Nmj+rj(t),rj(0)eiωtdt.
(5)

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

F(ω)=j=1Nmj+rj̇(t),rj̇(0)eiωtdt,
(6)

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 relations31 (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 

FIG. 2.

Velocity density of states (VDOS) of the DNA duplex CAAAAAG at equilibrium along with theoretical phonon branches. The peaks of the dispersion branches correspond, from lower to higher frequency, to ωbd/2π, ωt/2π, ωc/2π, and ωH2+ωbd2/2π from Eqs. (11), (16), (17), and (12). The VDOS is normalized to unity in the output frequency range and multiplied by a factor of 500.

FIG. 2.

Velocity density of states (VDOS) of the DNA duplex CAAAAAG at equilibrium along with theoretical phonon branches. The peaks of the dispersion branches correspond, from lower to higher frequency, to ωbd/2π, ωt/2π, ωc/2π, and ωH2+ωbd2/2π from Eqs. (11), (16), (17), and (12). The VDOS is normalized to unity in the output frequency range and multiplied by a factor of 500.

Close modal

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 = Hbend + Hhydr, where

Hbend=n=2Nkb2[(unun1)2+(vnvn1)2],
(7)
Hhydr=n=1NΔH0N(ea(unvn)/21)2,
(8)

correspond, respectively, to the elastic bending energy and hydrogen-bonding interactions of the molecule.38 The variables un 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 kb 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=unvn2) yields two decoupled equations of motion41 

mẍn=kb(xn+1+xn12xn),
(9)
mÿn=kb(yn+1+yn12yn)+2aΔH0N(eayn1)eayn,
(10)

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

ωbend2(q)=ωbd2sin2(ql2),
(11)
ωoop2(q)=ωH2+ωbend2(q),
(12)

where l is the distance between adjacent nucleotides on the same strand, q is the wavenumber (πn(N+1)l with n=1,,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 kb such that ωbd/2π11.5±0.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,

ωbend=ωbdsin(π2(N+1)),
(13)

that takes values ωbend/2π of 2.24±0.1, 2.00±0.1, and 1.64±0.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 ωH/2π92±1 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, Htors and Hcompr, corresponding to the torsional and compressional energies of the molecule,46,47

Htors=n=2Nkθ2(θn,n1θ0)2,
(14)
Hcompr=n=2Nkc2(zn,n1z0)2,
(15)

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 zn,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 z0 the separation at the equilibrium position) and kθ and kc 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.

FIG. 3.

Schematic view of the model considered.

FIG. 3.

Schematic view of the model considered.

Close modal

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

ωtors2=ωt2sin2(ql2),
(16)
ωcompr2=ωc2sin2(ql2),
(17)

where ωt and ωc are frequencies related to the torsional (kθ/IM) and compressional (kc/m) stiffnesses of the molecule, with IM 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 kc 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.

From Eqs. (2) and (3), the change in vibrational entropy can be written as

ΔS=kBln(ω0(2)ω0(Ñ)ωb(2)ωb(Ñ))+kBln(hkBT|ωb(0)|ω0(1)2πγ)(ΔSbend+ΔStors+ΔScompr+ΔSoop)+ΔSrc,
(18)

where ΔSbend, ΔStors, ΔScompr, and ΔSoop correspond to the change in bending, torsional, compressional, and out-of-phase breathing vibration frequencies, while ΔSrc 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 (ΔScompr0) 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

ΔSbend=N2kBln(lpdslpss),
(19)

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

lpds=ωbend2ML3(p1L)4kBT,
(20)

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 p1 is an inverse length such that p1L = 4.73. This treatment of the persistence length yields average values of 9±2, 13±3, and 25±2 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

FIG. 4.

(a) Persistence lengths of the double stranded molecules. The persistence length was calculated using Eq. (20) (circle) and those parameters from Table I. Experimental persistence lengths (diamond) were calculated in Ref. 53, where an estimate of the persistence length from simulations for the duplex GCATCTGGGC is also provided (red square). Dashed lines correspond to persistence lengths calculated from the WLC model (blue from Ref. 53, green from Ref. 54). (b) Persistence lengths of the single stranded molecules calculated from molecular dynamics simulations using the end-to-end (ETE) distance.53 Dashed lines correspond to persistence lengths calculated in Ref. 55.

FIG. 4.

(a) Persistence lengths of the double stranded molecules. The persistence length was calculated using Eq. (20) (circle) and those parameters from Table I. Experimental persistence lengths (diamond) were calculated in Ref. 53, where an estimate of the persistence length from simulations for the duplex GCATCTGGGC is also provided (red square). Dashed lines correspond to persistence lengths calculated from the WLC model (blue from Ref. 53, green from Ref. 54). (b) Persistence lengths of the single stranded molecules calculated from molecular dynamics simulations using the end-to-end (ETE) distance.53 Dashed lines correspond to persistence lengths calculated in Ref. 55.

Close modal

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 distance53 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±0.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θ(N1)α, where N − 1 is the number of melted bases and α=2.2±0.1. Assuming the change in inertia is also included in this formula, this estimate yields a torsional entropy

ΔStors=N2kBαln(N1).
(21)

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

mẍn=kb(xn+1+xn12xn),
(22)
mÿn=kb(yn+1+yn12yn).
(23)

Out-of-phase dynamics and in-phase dynamics share the same bending stiffness kb; 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ωbd2v=ωv, where A is the N×N matrix

A=1100121001200001,
(24)

whose jth column is zero, ωbd 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 kb 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 kb at the barrier through the product of kb at the equilibrium position and lpss/lpds. The frequency ωbd can hence be estimated from the double stranded bending frequency at equilibrium through

ωbd=ωbdlpss/lpds,
(25)

leading to values ωbd/2π of 3.32±2, 2.76±1.7, and 2.13±1.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 jth 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 ωbdN1, independent of the position of the hinge. We note that the geometric average of these out-of-phase vibration frequencies ωbd/2π (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

ΔSoop(N1)kBln(ωHωbd),
(26)

where the out-of-phase branch at the equilibrium position is approximated by a constant ωH as ωHω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 jth 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

ΔSrc=kBln(hkBTωH|ωb(0)|2πγ),
(27)

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, ωb0=ω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:

ΔS=ΔSoop+ΔStors+ΔSbend+ΔSrc(N1)kBln(ηωHωbd)+1.1NkBln(N1)+NkBln(η)+kBln(hkBTωH222πγ),
(28)

where near-equilibrium all-atom simulations have provided universal values for the hydrogen bond vibration frequencies ωH/2π92±1 cm−1 and the characteristic bending frequency ωbd/2π11.5±0.5 cm−1. The damping factor is roughly γ/2π = 0.5 ps−1, leading to a correction term TΔSrc± of barely 0.5 kBT. For our sequences, the persistence length ratio between double stranded and single-stranded DNA is estimated by

η2=lpdslpss=(mωbd2l3lpss×kBT)(N4500.5)sin2(π2(N+1)),
(29)

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±7, 17.3±10.0, and 33.3±17 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.

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 measurements7,14,42 [see Fig. 5(a)]. The experimental kinetic rates vary over 4 orders of magnitude, as seen in the supplementary material.

FIG. 5.

(a) Comparison of activation free energies in kcal/mol from experiments and theory for different DNA sequences at different temperatures. Equation (29) was used to estimate the ratio between the persistence lengths. The numbers next to each data point mark the temperature of each measurement in K. The black dashed diagonal marks the ideal collapse of experiments and theory. Red dashed lines mark the cutoff of 2 kcal/mol from the center line. (b) Contribution of each entropy term from Eq. (18) to the activation entropy of the DNA duplex CAAAAAG. From left to right, out-of-phase entropy (black), torsional entropy (green), bending entropy (blue), and correction due to curvature effects (red). (1 kBT0.5922 kcal/mol at 298.15 K.)

FIG. 5.

(a) Comparison of activation free energies in kcal/mol from experiments and theory for different DNA sequences at different temperatures. Equation (29) was used to estimate the ratio between the persistence lengths. The numbers next to each data point mark the temperature of each measurement in K. The black dashed diagonal marks the ideal collapse of experiments and theory. Red dashed lines mark the cutoff of 2 kcal/mol from the center line. (b) Contribution of each entropy term from Eq. (18) to the activation entropy of the DNA duplex CAAAAAG. From left to right, out-of-phase entropy (black), torsional entropy (green), bending entropy (blue), and correction due to curvature effects (red). (1 kBT0.5922 kcal/mol at 298.15 K.)

Close modal

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 kc may be a strong function of the secondary structure and we shall include nonlinear corrections to our Hamiltonian Hcompr, as was done for the hydrogen bond Hamiltonian Hhydr.

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

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.

From the Wiener-Khintchine theorem,59,60 the power spectrum from Eq. (6) can be written as

F(ω)=j=1Nmj|aj(ω)|2,
(A1)

where aj(ω)=+ṙj(t)eiωtdt. The displacement of the kth coordinate of the atom j can be decomposed by means of the normal modes λ through

rkj(t)=λAλkjeiωλt,
(A2)

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

F(ω)=λn=13N|Aλm|2ωλ2δ(ω+ωλ),
(A3)

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λm|ωλ2=kBT. Hence, as the density of states is given by Ref. 61 ρ(ω)=λδ(ω+ωλ), the power spectrum can be written as F(ω) = 3NkB(ω).

For phonons in a crystal lattice, the density of states presents singularities where the group velocity is zero62 (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.

DNA structures were generated using the NAB module of Amber 14.63 Each structure was solvated in a rectangular box with explicit TIP3P water molecules64 (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) method66 and a real space cutoff radius of 9 Å.

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 Å was applied to short-range, non-bonded interactions, and 9.0 Å 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 analysis67 (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.

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=kBThexp(ΔGkBT), where ΔG  =  ΔHTΔS. The logarithmic form of this equation yields

ln(k)=ln(kBTh)+ΔSkBΔHkBT.
(C1)

Given the minimal temperature dependence of ln(kBTh), a plot of ln(k) versus 1/T yields ΔH and ΔS 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,

k=konKD,
(C2)

where kon is the hybridization rate and KD 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).

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 (ωH/2π[91,93] cm−1, ωbd/2π[11.0,12.0] cm−1, and η2[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 ΔS0 as the activation entropy evaluated with our suggested set of parameters and ΔS1 as the activation entropy evaluated with a different set, from Eq. (28) we may write

ΔS1ΔS0(2N1)kBln(η1η0)+(N+1)kBln(ωH1ωH0)(N1)kBln(ωbd1ωbd0),
(D1)

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

ΔSdiffkBln(η12η02)+kBln(ωH1ωH0)kBln(ωbd1ωbd0).
(D2)

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

Defining the parameter μ=η12/η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 kBT per base pair (net differences of up to 9 kBT 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±0.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 kBT per base pair, leading to net deviations from our suggested free energies close to 3.15 kBT, 3.6 kBT, and 4.5 kBT 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) (1kBT0.5922 kcal/mol at 298.15 K).

1.
A.
Egatz-Gomez
,
C.
Wang
,
F.
Klacsmann
,
Z.
Pan
,
S.
Marczak
,
Y.
Wang
,
G.
Sun
,
S.
Senapati
, and
H.-C.
Chang
, “
Future microfluidic and nanofluidic modular platforms for nucleic acid liquid biopsy in precision medicine
,”
Biomicrofluidics
10
,
032902
(
2016
).
2.
J.
SantaLucia
, Jr.
and
D.
Hicks
, “
The thermodynamics of DNA structural motifs
,”
Annu. Rev. Biophys. Biomol. Struct.
33
,
415
440
(
2004
).
3.
J.
SantaLucia
, “
A unified view of polymer, dumbbell, and oligonucleotide DNA nearest-neighbor thermodynamics
,”
Proc. Natl. Acad. Sci. U. S. A.
95
,
1460
1465
(
1998
).
4.
A.
Perez
and
M.
Orozco
, “
Real-time atomistic description of DNA unfolding
,”
Angew. Chem., Int. Ed.
49
,
4805
4808
(
2010
).
5.
L. E.
Morrison
and
L. M.
Stols
, “
Sensitive fluorescence-based thermodynamic and kinetic measurements of DNA hybridization in solution
,”
Biochemistry
32
,
3095
3104
(
1993
).
6.
A. P.
Williams
,
C. E.
Longfellow
,
S. M.
Freier
,
R.
Kierzek
, and
D. H.
Turner
, “
Laser temperature-jump, spectroscopic, and thermodynamic study of salt effects on duplex formation by dGCATGC
,”
Biochemistry
28
,
4283
4291
(
1989
).
7.
J. W.
Nelson
and
I.
Tinoco
, Jr.
, “
Comparison of the kinetics of ribooligonucleotide, deoxyribooligonucleotide, and hybrid oligonucleotide double-strand formation by temperature-jump kinetics
,”
Biochemistry
21
,
5289
5295
(
1982
).
8.
I.-F.
Cheng
,
S.
Senapati
,
X.
Cheng
,
S.
Basuray
,
H.-C.
Chang
, and
H.-C.
Chang
, “
A rapid field-use assay for mismatch number and location of hybridized DNAs
,”
Lab Chip
10
,
828
831
(
2010
).
9.
S.
Marczak
,
S.
Senapati
,
Z.
Slouka
, and
H.-C.
Chang
, “
Induced nanoparticle aggregation for short nucleic acid quantification by depletion isotachophoresis
,”
Biosens. Bioelectron.
86
,
840
848
(
2016
).
10.
O.
Dahlen
and
T. S.
van Erp
, “
Mesoscopic modeling of DNA denaturation rates: Sequence dependence and experimental comparison
,”
J. Chem. Phys.
142
,
235101
(
2015
).
11.
T. E.
Ouldridge
,
P.
Šulc
,
F.
Romano
,
J. P.
Doye
, and
A. A.
Louis
, “
DNA hybridization kinetics: Zippering, internal displacement and sequence dependence
,”
Nucleic Acids Res.
41
,
8886
8895
(
2013
).
12.
D. M.
Hinckley
,
G. S.
Freeman
,
J. K.
Whitmer
, and
J. J.
de Pablo
, “
An experimentally-informed coarse-grained 3-site-per-nucleotide model of DNA: Structure, thermodynamics, and dynamics of hybridization
,”
J. Chem. Phys.
139
,
144903
(
2013
).
13.
X.-B.
Gu
,
S.-i.
Nakano
, and
N.
Sugimoto
, “
Consecutive gc base pairs determine the energy barrier of DNA duplex formation under molecularly crowded conditions
,”
Chem. Commun.
26
,
2750
2752
(
2007
).
14.
S.
Howorka
,
L.
Movileanu
,
O.
Braha
, and
H.
Bayley
, “
Kinetics of duplex formation for individual DNA strands within a single protein nanopore
,”
Proc. Natl. Acad. Sci. U. S. A.
98
,
12996
13001
(
2001
).
15.
T.
Ohmichi
,
H.
Nakamuta
,
K.
Yasuda
, and
N.
Sugimoto
, “
Kinetic property of bulged helix formation: Analysis of kinetic behavior using nearest-neighbor parameters
,”
J. Am. Chem. Soc.
122
,
11286
11294
(
2000
).
16.
A.
Madhumalar
and
M.
Bansal
, “
Structural insights into the effect of hydration and ions on a-tract DNA: A molecular dynamics study
,”
Biophys. J.
85
,
1805
1816
(
2003
).
17.
M.
Thomas
,
M.
Brehm
,
R.
Fligg
,
P.
Vöhringer
, and
B.
Kirchner
, “
Computing vibrational spectra from ab initio molecular dynamics
,”
Phys. Chem. Chem. Phys.
15
,
6608
6622
(
2013
).
18.
K.-Y.
Wong
and
B. M.
Pettitt
, “
The pathway of oligomeric DNA melting investigated by molecular dynamics simulations
,”
Biophys. J.
95
,
5618
5626
(
2008
).
19.
T. S.
Van Erp
,
S.
Cuesta-Lopez
, and
M.
Peyrard
, “
Bubbles and denaturation in DNA
,”
Eur. Phys. J. E
20
,
421
434
(
2006
).
20.
J. S.
Schreck
,
T. E.
Ouldridge
,
F.
Romano
,
P.
Šulc
,
L. P.
Shaw
,
A. A.
Louis
, and
J. P.
Doye
, “
DNA hairpins destabilize duplexes primarily by promoting melting rather than by inhibiting hybridization
,”
Nucleic Acids Res.
43
,
6181
6190
(
2015
).
21.
D. Z.
Avizonis
and
D. R.
Kearns
, “
Structural characterization of d(CAACCCGTTG) and d(CAACGGGTTG) mini-hairpin loops by heteronuclear NMR: The effects of purines versus pyrimidines in DNA hairpins
,”
Nucleic Acids Res.
23
,
1260
1268
(
1995
).
22.
S.
Roy
,
S.
Weinstein
,
B.
Borah
,
J.
Nickol
,
E.
Appella
,
J. L.
Sussman
,
M.
Miller
,
H.
Shindo
, and
J. S.
Cohen
, “
Mechanism of oligonucleotide loop formation in solution
,”
Biochemistry
25
,
7417
7423
(
1986
).
23.
R.
Landauer
and
J.
Swanson
, “
Frequency factors in the thermally activated process
,”
Phys. Rev.
121
,
1668
(
1961
).
24.
L. I.
McCann
,
M.
Dykman
, and
B.
Golding
, “
Thermally activated transitions in a bistable three-dimensional optical trap
,”
Nature
402
,
785
787
(
1999
).
25.
P.
Hanggi
, “
Escape from a metastable state
,”
J. Stat. Phys.
42
,
105
148
(
1986
).
26.
G. H.
Vineyard
, “
Frequency factors and isotope effects in solid state rate processes
,”
J. Phys. Chem. Solids
3
,
121
127
(
1957
).
27.
H.
Hu
and
W.
Yang
, “
Free energies of chemical reactions in solution and in enzymes with ab initio quantum mechanics/molecular mechanics methods
,”
Annu. Rev. Phys. Chem.
59
,
573
601
(
2008
).
28.
K.
Wendler
,
M.
Brehm
,
F.
Malberg
,
B.
Kirchner
, and
L.
Delle Site
, “
Short time dynamics of ionic liquids in AIMD-based power spectra
,”
J. Chem. Theory Comput.
8
,
1570
1579
(
2012
).
29.
M.-P.
Gaigeot
,
M.
Martinez
, and
R.
Vuilleumier
, “
Infrared spectroscopy in the gas and liquid phase from first principle molecular dynamics simulations: Application to small peptides
,”
Mol. Phys.
105
,
2857
2878
(
2007
).
30.
N.
Papanicolaou
,
I.
Lagaris
, and
G.
Evangelakis
, “
Modification of phonon spectral densities of the (001) copper surface due to copper adatoms by molecular dynamics simulation
,”
Surface Sci.
337
,
L819
L824
(
1995
).
31.
M.
Krisch
,
A.
Mermet
,
H.
Grimm
,
V.
Forsyth
, and
A.
Rupprecht
, “
Phonon dispersion of oriented DNA by inelastic x-ray scattering
,”
Phys. Rev. E
73
,
061909
(
2006
).
32.
D.
Lin
,
A.
Matsumoto
, and
N.
Go
, “
Normal mode analysis of a double-stranded DNA dodecamer d(CGCGAATTCGCG)
,”
J. Chem. Phys.
107
,
3684
3690
(
1997
).
33.
K.
Irikura
,
B.
Tidor
,
B.
Brooks
, and
M.
Karplus
, “
Transition from B to Z DNA: Contribution of internal fluctuations to the configurational entropy difference
,”
Science
229
,
571
573
(
1985
).
34.
K.-C.
Chou
and
Y.-S.
Kiang
, “
The biological functions of low-frequency vibrations (phonons): 5. A phenomenological theory
,”
Biophys. Chem.
22
,
219
235
(
1985
).
35.
K.-C.
Chou
, “
Biological functions of low-frequency vibrations (phonons). III. Helical structures and microenvironment
,”
Biophys. J.
45
,
881
889
(
1984
).
36.
F.
Merzel
,
F.
Fontaine-Vive
,
M.
Johnson
, and
G.
Kearley
, “
Atomistic model of DNA: Phonons and base-pair opening
,”
Phys. Rev. E
76
,
031917
(
2007
).
37.
T.
Dauxois
, “
Dynamics of breather modes in a nonlinear ‘helicoidal’ model of DNA
,”
Phys. Lett. A
159
,
390
395
(
1991
).
38.
N.
Theodorakopoulos
, “
Melting of genomic DNA: Predictive modeling by nonlinear lattice dynamics
,”
Phys. Rev. E
82
,
021905
(
2010
).
39.
M.
Peyrard
and
G.
James
, “
Intrinsic localized modes in nonlinear models inspired by DNA
,”
Nonlinear Theory Its Appl., IEICE
3
,
27
51
(
2012
).
40.
T.
Dauxois
,
M.
Peyrard
, and
A. R.
Bishop
, “
Entropy-driven DNA denaturation
,”
Phys. Rev. E
47
,
R44
(
1993
).
41.
S.
Zdravković
, “
Helicoidal Peyrard–Bishop model of DNA dynamics
,”
J. Nonlinear Math. Phys.
18
,
463
484
(
2011
).
42.
B.
Rauzan
,
E.
McMichael
,
R.
Cave
,
L. R.
Sevcik
,
K.
Ostrosky
,
E.
Whitman
,
R.
Stegemann
,
A. L.
Sinclair
,
M. J.
Serra
, and
A. A.
Deckert
, “
Kinetics and thermodynamics of DNA, RNA, and hybrid duplex formation
,”
Biochemistry
52
,
765
772
(
2013
).
43.
U.
Christensen
,
N.
Jacobsen
,
V. K.
Rajwanshi
,
J.
Wengel
, and
K.
Troels
, “
Stopped-flow kinetics of locked nucleic acid (LNA)–oligonucleotide duplex formation: Studies of LNA–DNA and DNA–DNA interactions
,”
Biochem. J.
354
,
481
484
(
2001
).
44.
M.
González-Jiménez
,
G.
Ramakrishnan
,
T.
Harwood
,
A. J.
Lapthorn
,
S. M.
Kelly
,
E. M.
Ellis
, and
K.
Wynne
, “
Observation of coherent delocalized phonon-like modes in DNA under physiological conditions
,”
Nat. Commun.
7
,
11799
(
2016
).
45.
F.
Merzel
and
M. R.
Johnson
, “
Low-frequency vibrations of DNA and base pair opening
,”
Acta Chim. Slov.
58
,
442
447
(
2011
).
46.
F.
Sicard
,
N.
Destainville
, and
M.
Manghi
, “
DNA denaturation bubbles: Free-energy landscape and nucleation/closure rates
,”
J. Chem. Phys.
142
,
034903
(
2015
).
47.
A. K.
Dasanna
,
N.
Destainville
,
J.
Palmeri
, and
M.
Manghi
, “
Slow closure of denaturation bubbles in DNA: Twist matters
,”
Phys. Rev. E
87
,
052703
(
2013
).
48.
M.
Zoli
, “
Twist versus nonlinear stacking in short DNA molecules
,”
J. Theor. Biol.
354
,
95
104
(
2014
).
49.
M.
Barbi
,
S.
Cocco
, and
M.
Peyrard
, “
Helicoidal model for DNA opening
,”
Phys. Lett. A
253
,
358
369
(
1999
).
50.
S. B.
Smith
and
L.
Finzi
, “
Direct mechanical measurements of the elasticity of single DNA molecules by using magnetic beads
,”
Science
258
,
1122
1126
(
1992
).
51.
A.
Matsumoto
and
W. K.
Olson
, “
Sequence-dependent motions of DNA: A normal mode analysis at the base-pair level
,”
Biophys. J.
83
,
22
41
(
2002
).
52.
A.
Matsumoto
and
N.
Go
, “
Dynamic properties of double-stranded DNA by normal mode analysis
,”
J. Chem. Phys.
110
,
11070
11075
(
1999
).
53.
Y.-Y.
Wu
,
L.
Bao
,
X.
Zhang
, and
Z.-J.
Tan
, “
Flexibility of short DNA helices with finite-length effect: From base pairs to tens of base pairs
,”
J. Chem. Phys.
142
,
125103
(
2015
).
54.
C.
Yuan
,
H.
Chen
,
X. W.
Lou
, and
L. A.
Archer
, “
DNA bending stiffness on small length scales
,”
Phys. Rev. Lett.
100
,
018102
(
2008
).
55.
F.-H.
Wang
,
Y.-Y.
Wu
, and
Z.-J.
Tan
, “
Salt contribution to the flexibility of single-stranded nucleic acid of finite length
,”
Biopolymers
99
,
370
381
(
2013
).
56.
A. Y.
Sim
,
J.
Lipfert
,
D.
Herschlag
, and
S.
Doniach
, “
Salt dependence of the radius of gyration and flexibility of single-stranded DNA in solution probed by small-angle x-ray scattering
,”
Phys. Rev. E
86
,
021901
(
2012
).
57.
N.
Nakanishi
,
A.
Nagasawa
, and
Y.
Murakami
, “
Lattice stability and soft modes
,”
J. Phys. Colloq.
43
,
C4-35
(
1982
).
58.
G.
Ackland
, “
Elastic properties as a pointer to phase transitions
,”
RIKEN Rev.
29
,
34
38
(
2000
).
59.
N.
Wiener
, “
Generalized harmonic analysis
,”
Acta Math.
55
,
117
258
(
1930
).
60.
A.
Khintchine
, “
Korrelationstheorie der stationären stochastischen prozesse
,”
Math. Ann.
109
,
604
615
(
1934
).
61.
H. M.
Flores Ruiz
, “
Modos vibracionales de baja frecuencia y su impacto en la formacion de vidrios
,” Ph.D. thesis,
Universidad Nacional Autonoma de Mexico
,
2012
.
62.
C.
Kittel
,
Introduction to Solid State Physics
(
Wiley
,
2005
).
63.
D. A.
Case
,
V.
Babin
,
J.
Berryman
,
R.
Betz
,
Q.
Cai
,
D.
Cerutti
,
T.
Cheatham
 III
,
T.
Darden
,
R.
Duke
,
H.
Gohlke
 et al, “Amber 14,”
2014
.
64.
W. L.
Jorgensen
,
J.
Chandrasekhar
,
J. D.
Madura
,
R. W.
Impey
, and
M. L.
Klein
, “
Comparison of simple potential functions for simulating liquid water
,”
J. Chem. Phys.
79
,
926
935
(
1983
).
65.
A.
Pérez
,
I.
Marchán
,
D.
Svozil
,
J.
Sponer
,
T. E.
Cheatham
,
C. A.
Laughton
, and
M.
Orozco
, “
Refinement of the amber force field for nucleic acids: Improving the description of α/γ conformers
,”
Biophys. J.
92
,
3817
3829
(
2007
).
66.
T.
Darden
,
D.
York
, and
L.
Pedersen
, “
Particle mesh Ewald: An N·log(N) method for Ewald sums in large systems
,”
J. Chem. Phys.
98
,
10089
10092
(
1993
).
67.
D. R.
Roe
and
T. E.
Cheatham
 III
, “
PTRAJ and CPPTRAJ: Software for processing and analysis of molecular dynamics trajectory data
,”
J. Chem. Theory Comput.
9
,
3084
3095
(
2013
).
68.
W.
Humphrey
,
A.
Dalke
, and
K.
Schulten
, “
VMD: Visual molecular dynamics
,”
J. Mol. Graphics
14
,
33
38
(
1996
).

Supplementary Material