Comparative melting and healing of B-DNA and Z-DNA by an infrared laser pulse

Recently, free electron lasers have been used for melting biomolecular complexes.^1–5 Kawasaki and co-workers^3–5 have developed a mid-infrared free electron laser with highly specific oscillation characteristics having a high photon density, a picosecond pulse structure, and a range of tunable frequencies. To date, such a laser pulse with frequencies tuned to the amide I bands has been used in experiments to dissociate amyloid-like fibrils of lyzozyme into its native forms,³ convert insulin fibrils into soluble monomers,⁴ and dissociate a short human thyroid hormone peptide.⁵ A simulated laser pulse has also been used in atomistic molecular dynamics simulations to dissociate amyloid fibrils⁶ and viruses⁷ and to study the break-up of a peptide-based nanotube.⁸ Aside from investigating the dissociation of selected biomolecules, in this work we propose that the application of a simulated laser pulse could provide for new computational opportunities to probe the relative response and, quite possibly, the relative stability of similar structures, such as nucleic acids with the same sequence but with different secondary structures.

A crucial aspect for the understanding of nucleic acid structures and their function is the relative stability of different competing structures that are involved in cellular processes. The discernment of the relative structural stabilities of different nucleic acid structures, both in vitro and in silico, is rather challenging. A time honored technique to address this issue is to carry out a thermal denaturation experiments,⁹ where a sample is heated beyond the melting point. This causes a conformational transition in the molecule that is measured by recording a specific variable as a function of temperature. Typical experiments include the recording of the absorbance at 260 nm in a UV-visible spectrophotometer;¹⁰ circular dichroism spectra;^11,12 fluorescence spectra;¹³ Raman signals;^14,15 or NMR measurements.^16,17

Thermal denaturing experiments can also be carried out computationally, but fully atomistic melting simulations are extremely expensive. Instead, simplified models have been introduced for estimating the melting temperature of DNA duplexes, based on empirical or statistical thermodynamics models^18–34 that predict the stability of nucleic acid secondary structure, including RNA.^22,35–48 Unfortunately, these models are still empirical and can fail to capture subtle sequence-dependent effects, or unusual conformations different from A and B duplexes, or conformations with noncanonical or mismatched base pairs, etc. A state-of-the-art web-based tool for predicting fluorescent high-resolution DNA melting curves and denaturation profiles of PCR products is the software package uMELT.⁴⁹ This tool, however, can only be applied to standard DNA duplexes. Since secondary structures other than B-DNA can play crucial roles in processes such as gene expression, extending the computational predictive capacity to these unusual structures is a desirable goal.

In this work, we use fast melting by an infrared laser pulse of B-DNA and Z-DNA duplexes with sequence d(5′-CGCGCGCGCGCG-3′)₂ in order to compare the melting and healing response of the duplexes and to explore whether this response reflects their relative stability. The left-handed, double-helix Z-DNA with two antiparallel chains joined by Watson-Crick (WC) base pairs was first crystallized in 1979.⁵⁰ Z-DNA is favored by sequences that alternate purines and pyrimidines, mainly CG or GC. The left-handed helix is thinner and more rigid than B-DNA, and while the pyrimidine-purine base pairs are in anti-anti conformation in B-DNA, the guanine rotates around its glycosidic bond in Z-DNA, resulting in an anti-syn configuration.^51,52 As a consequence of these rotations, the phosphate groups become closer in Z-DNA; the sugar-phosphate backbone displays the characteristic zig-zag pattern that gave rise to the name Z-DNA; and the CpG and GpC steps stack differently, causing a dinucleotide step to be the repeating unit in Z-DNA. A high density of base sequences favoring Z-DNA is found near transcription start sites,⁵³ where Z-DNA is stabilized by negative supercoiling of DNA.^51,52,54 In the cell, Z-DNA is induced by a set of binding proteins near promoter regions, which boosts the transcription of downstream genes.⁵⁵ Z-DNA is highly immunogenic, and antibodies against it^56–58 are used to find locations prone to Z-DNA conformations. The current view is that Z-DNA formation plays a role in gene expression, regulation, and recombination.^{52,55,59–65}

B-DNA is more stable than Z-DNA under physiological ionic strength and pH conditions. While in the cell Z-DNA can be stabilized by the negative supercoiling^51,52,54 that, for instance, arises from the relocation of the RNA polymerase during transcription, and through binding with highly specific Z-DNA binding proteins such as ADAR1, E3L, and PKZ,^66–69 in vitro inducers and stabilizers are needed for Z-DNA, such as high ionic concentrations, base chemical modifications, organic solvents, divalent cations and transition metal complexes, and small molecular complexes. (these are reviewed in Ref. 70). Computationally, the determination of the relative stability of B-DNA and Z-DNA is not trivial. Existing predictive tools do not have provisions for this alternative structure, and therefore one has to resort to very long time simulations or careful theoretical and statistical modeling. Considerable insight has been obtained by studying the microscopics behind the B-Z DNA transition,^50,71–81 and the ion distribution around these nucleic acid structures.⁸² Recent molecular dynamics (MD) simulations indicate that the transition is governed by a complex free energy landscape which allows for the coexistence of several competing mechanisms so that the transition is better described in terms of a reaction path ensemble.⁸¹ 

In this work, we choose B-DNA and Z-DNA duplexes with sequence d(5′-CGCGCGCGCGCG-3′)₂ to test their behavior when exposed to fast melting by an infrared laser pulse. In addition to the intrinsic interest of the melting and healing process, we propose that this technique could be used as a relatively inexpensive tool to determine relative stability in fully solvated, atomically accurate nucleic acid structures. Naturally, a laser pulse implies a nonequilibrium process and the results cannot be used to construct an equilibrium melting curve, and therefore the ensuing free energy estimates cannot be obtained either. The perturbation, however, can map out the responses of the structures for the entire range of applied fields, covering extremely small perturbations all the way to complete melting. The response can be measured by a wide range of variables including normalized handedness, base stacking, Watson-Crick hydrogen bonds, etc. A systematic trend in all these variables (in this case that Z-DNA is more “melted” than B-DNA) for all strengths of the field is interpreted as being indicative of the greater stability of B-DNA under the solvation conditions in the simulations. A traditional equilibrium molecular dynamics simulation not only is orders of magnitude more expensive but also faces additional challenges. Since B-DNA is more prone to fray at the ends than Z-DNA, the melting of short B-DNA oligomers is more susceptible to length effects than Z-DNA, a difficulty that can be partially bypassed with the laser melting technique. In addition, Z-DNA is believed to be stabilized by higher temperatures,^73,83–85 which complicates the interpretation of the traditional melting experiments.

Simulations of B-DNA and Z-DNA duplexes with sequence (5′-CGCGCGCGCGCG-3′)₂ in TIP3P waters⁸⁶ were carried out using the GROMACS 4.5.5 package⁸⁷ with the AMBER99SB-ILDN force field and DNA parameters ff99BSCO.⁸⁸ The initial DNA configurations were taken from our previous studies.^81,82 For each DNA (B-DNA or Z-DNA) structure, we performed simulations at five different NaCl salt concentrations: 0, 1.0, 2.5, 3.0, and 4.0M, with ion parameters by Joung and Cheatham.⁸⁹ Specifically, 22 Na⁺ neutralizing ions were added, followed by an additional number of Na⁺ and Cl⁻.

Ions to make up the specified NaCl salt concentration. The DNA was placed in an octahedral box containing about 7050 water molecules with a 1 nm distance from the solute to the box boundary, resulting in a box size and volume of 6.6 nm and 200 nm³, respectively. Periodic boundary conditions were imposed with a 1 nm cutoff for the van der Waals and electrostatic interactions. The long-range electrostatic interactions were computed with the Particle-Mesh Ewald summation method,⁹⁰ and the non-bonded interaction pair list updated every 10 fs. The covalent bonds were constrained using the LINCS algorithm⁹¹ with a relative geometrical tolerance of 10⁻⁴. The simulations were primarily NPT and made use of the Berendsen pressure coupling method⁹² to model the barostat (pressure was kept at 1 atm). To ensure stability, a time step of 0.2 fs was used, and the equations of motion were integrated using the Leapfrog method.⁹³ To control the temperature, we used a V-rescaling temperature coupling scheme based on rescaling the velocity with a stochastic term.⁹⁴ Since a laser pulse was applied to DNA, only the solvent molecules were coupled to the heat bath which was kept at 310 K, with a temperature coupling constant of 0.1 ps. Experimentally, the DNA samples would be immersed in an effectively “infinite” reservoir of solvent (perhaps even with added flow). However, the simulations are limited in terms of the sizes the one can deal with, even with current supercomputers necessitating the use of a thermostat in order to damp out the thermal fluctuations of the solvent.

The laser pulse was modeled by the following equation:

with E₀ denoting the amplitude of the electric field, σ the pulse width, t the time, t₀ the time when the pulse is maximal, ω the frequency, and c the velocity of light. An important issue that needs consideration is the direction of the laser pulse with respect to the orientation of the DNA molecule. While a given laser beam is polarized in a specified direction, the free electron laser experiments typically take place in an unpolarized environment. This is important, because the applied laser pulse will be most efficient in exciting a given bond if the direction of the electric field is along the direction of the bond. In this study, the native hydrogen bonds (H-bonds) (as listed in supplementary material Table S2) connecting the C–G base pairs (corresponding to bond indices 14 to 33 in supplementary material Table S1) are targeted, which are oriented almost perpendicular to the DNA axis.¹⁰¹ If the direction of the laser pulse (LPD) is parallel to the DNA axis, the laser pulse will have only a minimal effect on these bonds, even when the correct resonant frequency is used. Moreover, during the course of the simulation, the orientation of the DNA molecule may shift. To solve this issue, we used the following averaging procedure. For the simulations scanning the different frequencies, ten different LPDs were set and the results averaged. Specifically, the LPDs were set in parallel to the ten internal vectors $C1i′−C125−i′$⁠, with i representing the base index as shown in Fig. 1, and $C1i′$ the C1′ atom of the ith base of the starting configuration. For the laser melting simulation, averaging was performed as follows. First, the average vector (AV) of the ten $C1i′−C125−i′$ was calculated. Then, the angle between this AV vector and a given $C1i′−C125−i′$ was divided into four equal parts, thereby defining five different directions. Given the ten different $C1′−C125−i′$ vectors, this defines a total of fifty different LPDs. Simulations were performed for each of these directions and the results were then averaged.

For this project we carried out three kinds of simulations: equilibration, frequency scanning, and laser melting simulations. To equilibrate the solvated DNA (both with and without salt), we carried out NPT simulations at 310 K for 30 ns. From these simulations, 100 independent conformations were selected at random from the last 10 ns. These configurations provided us with initial configurations for the subsequent simulations. Following equilibration, we carried out laser scanning simulations in order to select the best frequency for the laser melting simulations. These simulations were carried out for both B- and Z-DNA in the absence of excess salt with laser parameters E₀, t₀ and σ set, respectively, to 1 V/nm, 10 ps, and 3 ps. The frequency ω was varied between 1600 to 1935 cm⁻¹, with frequency steps set to 5 cm⁻¹. The latter was lowered to 1 cm⁻¹ around the observed resonant peaks in order to provide for a better resolution. In total, 94 different frequencies were scanned. For each ω value, ten 30 ps simulations with different initial conditions selected from the 100 conformations in the equilibration step were then carried out, and the results averaged. From the mathematical expression of the laser pulse, we see that there are four parameters. Since the frequency ω is a property of the system, we actually have three parameters. For these simulations we chose two pulse widths, σ = 3 ps and σ = 6 ps and the two sets of simulations had the following parameters: (i) σ = 3 ps, t₀ = 10 ps, E₀ = (3.0, 4.0, 4.6, 5.0, 5.5, 6.0 and 6.5) V/nm, and total time for each run of 40 ps; (ii) σ = 6 ps, t₀ = 20 ps, E₀ = (3.0, 3.5, 4.0, 4.5 and 5.0) V/nm, and total time for each run of 100 ps. For each value of E₀, 50 simulations were carried out using different initial structures selected from the equilibration step. The values of the electric field used in our simulation range from about 6.0-20.0 J/cm² (depending on how one estimates the pulse shape), which is about 100 times higher than the current experimental values. Choosing such a high laser intensity is necessitated by the fact that experimentally a very large number of pulses are delivered over a relatively long period of time (microseconds). However, MD simulations are plagued by the so-called time scale problem, which precludes simulations over such long time scales. Here, we have chosen to add all the energy over a relatively short period of time.

The laser melting of the DNA structures was characterized in a number of standard ways. The DNA structure is held together by Watson-Crick H-bonds, so that their number represents a reasonable measure of the native structure. For both B- and Z-DNA, there are three H-bonds for each CG pair, giving a total of thirty-six H-bonds for the entire structure. Throughout the simulation, we tracked the percentage of the Watson-Crick H-bonds (WCHBs) as a function of time. As a further refinement, we also considered the base pair opening probability (BPOP), which represents the probability that a given base pair H-bond is broken.

Since the two major determinants of the stability of a helical duplex are the base-pair H-bonds and the base stacking energy, we have also defined a base-stacking index through the following procedure. A base is represented by a convex pentagon determined by five base atoms. In a C(G) base, these atoms are O2(O6), N4(N2), C5(N9), C6(C8), and N1(N7) (see Fig. 1). Given the stacking index (st) between two bases i and j (st_i,j), the base-stacking index (BSI) was calculated via the following equations:

In these equations, S_overlap is the overlapping area between successive bases, given by the projection of each pentagon plane on the next ones along the chain, r_i,j is the distance between the centers of mass of the two pentagons, and N is the total number of bases. Note that the choice of subindexes allows for the projection of each base pentagon onto both next ones.

Another useful measure of the duplex structure is given by its handedness (H).⁹⁵ Handedness is a natural choice here, since B- and Z-DNA are right-handed and left-handed helical structures. Our working definition of H is based on our previous work on polyproline peptides⁹⁵ and extended to transitions between B- and Z-DNA.⁸¹ It discriminates between left- (H < 0) and right-handed (H > 0) helical structures. For the DNA double helix, the position of the phosphorus (P) atoms of the backbone phosphate group was found to be a good choice for the definition of H.⁸¹ In brief, the definition of H for a portion of DNA between the base pairs n and m makes use of a sequence of P atoms: $Pn1,Pn2,Pn+11,Pn+12,…,Pm1,Pm2$⁠, where the upper index indicates the strand number (1 or 2) and the lower index indicates the base-pair number labeled in the 5′ → 3′ direction of strand 1. Note that this definition of H is independent of the labeling of the strands. Since the 5′ terminal is missing the backbone phosphate group, the first and last elements of the sequence are ignored in the definition of H. Therefore, the sequence used in the definition starts from $P12$ instead of $P11$ and ends with $PN1$ instead of $PN2$ (N being the total number of bases). The position of these P atoms then defines H via

in which each P_i is a point in the sequence discussed above, and

In this last equation, the points A, B, C, D define the vectors $AB⃗$ and $CD⃗$ and the midpoints of these vectors, called E and F, in turn form the vector $EF⃗$⁠.

The results of the laser frequency scans for DNA are shown in Fig. 2, which displays the percent fluctuation of all the different bond types for both B- and Z-DNA, as well as the change in energy ΔE, all as a function of the frequency. The fluctuations of the 34 different bond types (details of which are given in Table S2 in supplementary material)¹⁰¹ were computed as the average of the fluctuations of a given bond over a 6 ps time frame centered around the laser pulse maximum. This was then compared to the average fluctuation of the same bond in an equilibrium simulation in the absence of any laser pulse. The ΔE the system at a given laser frequency was defined to be the total energy difference (highest value) for a system with the laser minus the total energy of an equilibrium system. The ΔE spectrum displays a complex structure of approximately 7 peaks of varying height, which coincide with high levels of bond fluctuations. Both in terms of the bond fluctuations and ΔE, the spectrum for B- and Z-DNA track each other quite well, with only relatively small shifts. Table I lists the peak frequencies and the associated bond types for the DNA structures. Notably, the two highest peaks are at ω = 1659, 1870 cm⁻¹ and the associated fluctuating bond types are 22, 24, and 27 (former), and 22, 25, 26, 27, 28, and 29 for the latter. Clearly, the latter is associated with a larger number of bonds that are fluctuating. Moreover, bond type 28 is directly associated with a Watson-Crick H-bond. Given that ω = 1870 cm⁻¹ targets this important bond and that peak height for both B- and Z-DNA are almost identical, we chose to conduct our laser melting simulations at this frequency.

We now turn to the laser melting of the DNA helices. Both sets of simulations with pulse width σ = 3 ps and σ = 6 ps gave systematically consistent results, and here we report on the results obtained with the wider pulse. Selected results for σ = 3 ps are given in supplementary material.¹⁰¹ 

Thus, Fig. 3 shows the time evolution of the laser pulse with σ = 6 ps and t₀ = 20 ps, for two intensities of the field, E₀ = 3.0 V/nm and 4.5 V/nm. For both field values, the energy absorption by B-DNA and Z-DNA is exactly the same (except for noise), and the total energy shows a small lag (about 2 ps) of the duplex response with respect to the perturbation. Fig. 4 shows snapshots of both DNA structures following the application of the laser pulse. The duplexes stay almost intact immediately following the application of the laser pulse (which reaches its maximum at 20 ps). This is quickly followed by melting and considerable disruption of the H-bonds of the structure. After the pulse dies out at approximately 40 ps, the structures continue struggling with the unfolding provoked by the laser.

Fig. 5 shows two quantitative measures of the laser melting process: the percentage of Watson-Crick H-bonds (WCHBs, top panels) and base-pair RMSDs (bottom panels) taken with respect to the initial structure, for various electric field amplitudes. These quantities start to show the effect of the pulse at about 11-12 ps, with a subsequent sharp drop in the number of H-bonds and a sharp increase in the base-pair RSMDs that signal the disruption of the structure. The disruption becomes maximum at about 24-25 ps. After that, the duplex partially reforms, depending on the magnitude of the field. After about 40 ps, most of the rapid changes are over and there are only small changes. Naturally, the amount of unfolding depends on the applied field strength, which varies from E₀ = 3.0 V/nm, where the duplexes almost completely recover, to E₀ ≥ 5.0 V/nm, where the duplexes are denatured. Thus, for E₀ = 4.0 V/nm, about 43% of the H-bonds reform for B-DNA, and about 28% for Z-DNA; for E₀ = 4.5 V/nm, this drops to about 20% for B-DNA and 10% for Z-DNA. Note that in this regime, the percent WCHBs is always higher and the RMSD is smaller for B-DNA than for Z-DNA, reflecting the greater stability of the B-DNA structure (in zero or smaller fields B-DNA exhibits more end fraying than Z-DNA). These curves were averaged over 50 runs; we extended these simulations for up to 1 ns for only 10 of those 50 runs for each value of the electrical field (data not shown) and we found that essentially the values obtained at 100 ps remain constant or show a very slight increase. Thus, we can say that about 100 ps the system reaches a long-lived plateau. The WC hydrogen bond probabilities for each hydrogen bond as function of time for different magnitudes of the electric field at zero excess salt are shown in Figures S3 and S4 for B-DNA and Z-DNA, respectively, while the equilibrium distribution for various salts is shown in Figure S2.¹⁰¹ The terminal fraying of B-DNA both in equilibrium and at small fields is apparent. In both duplexes, the disruption of the bonds at around 20 ps for E₀ = 3.0,  3.5 V/nm is followed by the re-formation of the bonds, but this re-bonding quickly degrades as the magnitude of the field increases.

Fig. 6 shows the base-stacking index for B-DNA and Z-DNA at zero and 4M NaCl (in addition to the 22 neutralizing Na⁺). The initial base-stacking index is higher for B-DNA (0.25) than Z-DNA (0.15) reflecting the more stable base stacking of the right-handed duplex. This more favorable base stacking of B-DNA is maintained for all field strengths. Notice that the addition of 4M NaCl decreases base stacking for B-DNA while increasing it for Z-DNA. Finally, Fig. 7 shows the normalized handedness (nHDN) and BPOP of both DNA structures as a function of time for different field strengths at zero excess salt. The BPOP looks approximately as the mirror image of the WCHB images since, as the base pairs open, their WC H-bonds break. Lower values of BPOPs reflect higher stability of B-DNA with respect to Z-DNA. The nHDN is even more sensitive to the melting process, and is powerful to visually identify the more unwinded structure (the absolute value is shown here; handedness is positive for B-DNA and negative for Z-DNA). Fig. 8 gives a plot of the normalized handedness and base stacking index as a function of the magnitude of the electric field. B-DNA exhibits more favorable stacking throughout the entire range of electric field magnitudes, and B-DNA is considerably more resistant to unwinding by the laser pulse than Z-DNA.

We have also examined the behavior of all of these different structural measures as a function of excess salt concentration, as shown in Fig. 9. It is well known that increasing the salt concentration stabilizes the DNA duplexes, but stabilization is stronger in Z-DNA. Fig. 9 display this general trend, with the curves growing or decaying (depending on the quantity) towards more stability and exhibiting a greater slope for Z-DNA.

Finally, how was the running time chosen? Clearly, after the perturbation is gone the structures would tend to fold back if enough time (very long time if they need to completely refold) has elapsed. We find that for each magnitude of the field, the quantities characterizing the structures reach a plateau somewhere between 60 and 100 ps, depending on the quantity, as shown in Figs. 3, 5, 6, and 7. These figures were computed as an average over 50 runs for each field magnitude. We extended simulations for only 10 of each set of 50 runs up to 1 ns, and indeed the values at 100 ps are the same or extremely close to the values at 1 ns. Therefore, the run times were chosen in the first segment of the plateau, which is indicative of much slower dynamics. Thus, to build a comparative figure such as those in Fig. 8 one can choose the time where the pulse effectively goes to zero (at approximately 40 ps for σ = 6 ps) and add the “lag” in the molecular response (≃2 ps, such that quantities are compared at 42 ps). Alternatively, one could choose values somewhere in the plateau, say at 100 ps. In both cases, the same trends are observed, as shown in Fig. 8.

The classical laser pulse simulated in this work is typical of a free-electron based laser,^1–5 with specific oscillation characteristics of a picosecond pulse structure, tunable wavelengths within the infrared regime, and a high photon density. Using this laser pulse, it has been possible to target specific bonds within a given structure and thereby study the unfolding of the structure. Specifically, this method was used to investigate the dissociation of amyloid fibrils into soluble monomers,^3–6 the break-up of a peptide-based nanotube,⁸ and the dissociation of viruses.⁷ Comparing healing and, possibly, relative stability of different structures via laser induced melting is clearly not a universal technique. Key to the application of the technique is the fact that the perturbation must affect both structures in the same way (i.e., the damage to the structures must be the same). As an example, we have used this technique previously to compare the responses of polyasparagine and polyglutamine amyloid aggregates. In this case, the results of the laser frequency scanning simulations resulted in an optimum frequency of ω = 1671 cm⁻¹ that targeted the C = O main-chain bonds equally in both aggregates, thereby destabilizing the β sheets.⁹⁶ 

To find out the relative stability of two DNA (or RNA) structures with the same sequence but with different secondary structure can be quite expensive computationally. In this work, we set out to determine what can be inferred about the relative stability of two duplexes with exactly the same sequence, d(5′-CGCGCGCGCGCG-3′)₂ but with two very different structures, right-handed B-DNA and left-handed Z-DNA, that are subject to a laser pulse. A meaningful comparison can be carried out because it is possible to find a laser pulse frequency, mainly ω = 1870 cm⁻¹, that generates the same resonant peak affecting the same bonds, especially bond C6–O6 involving a Watson-Crick H-bond, for both structures (Fig. 2). This results in exactly the same energy absorption pattern for both B-DNA and Z-DNA (Fig. 3). The perturbation is tuned to vary from minor disruption at small field strength to complete melting at high field strength, which allows for an extensive comparison of the responses of the two duplexes. Since a laser pulse is by definition a nonequilibrium process, its results cannot be translated into an equilibrium melting curve with the ensuing free energy estimates. However, in all cases, Z-DNA lags behind B-DNA in the recovery process. Although this can be associated with the kinetics of refolding, we believe that because this occurs in the whole range of field strengths, from small perturbations to complete melting, the process actually reflects the lower stability of Z-DNA with respect to B-DNA under the solvation conditions employed in these simulations. In that sense, Fig. 8 can be interpreted as the corresponding non-equilibrium melting curves of the duplexes.

Finally, it should be noted that we explore the partial melting and subsequent healing of DNA strands by means of MD techniques. Given the time scale problem and the necessity of simulating relatively small system size, this requires the use of high electric field intensities and a temperature rescaling of the solvent. However, even if it were currently feasible to simulate much larger systems over very much longer periods of time, we believe that the system would still very much follow similar generic behavior of the kind that we have uncovered with our simulations.

Aside from uncovering the response of a biomolecular system to a fast, high intensity laser pulse, our simulations can viewed as representing a methodology in its own right, whose aim is to provide for a fast and cheap way of determining the relative stability of different DNA conformations. B-DNA is more stable than Z-DNA under physiological ionic strength and pH conditions. Experimentally, after a high nucleation barrier in a B- to Z-DNA transition, the propagation free energy of Z-DNA has been measured to be 0.66 kcal/mol per CG repeat at 0.1 NaCl,⁹⁷ and in general less than 2 kcal/mol per dinucleotide repeat.⁹⁸ Simulations of the base-extrusion zipper mechanism at 0.1M salt concentration in explicit water obtained 0.93 kcal/mol per dinucleotide,⁸⁰ while free energy maps with only neutralizing ions (no excess salt) resulted in 1.5 kcal/mol per dinucleotide.⁸¹ The computational measures required enhanced sampling techniques, with definition of reaction coordinates.

A melting experiment is conceptually more straightforward, but it faces some challenges. First, B-DNA is more flexible than Z-DNA, and the configurational entropy difference is important in stabilizing B-DNA.⁹⁹ However, due to its inherent flexibility, B-DNA is more prone to terminal fraying (Fig. S5)¹⁰¹ than more rigid Z-DNA, therefore B-DNA temperature melting is more susceptible to length effects than Z-DNA. While this does not affect long sequences, it could be a problem for relatively short oligomers. This is not an issue in vivo, where short DNA duplexes with four free ends are not encountered but could present a challenge for thermal stability measurements in vitro and in silico. Second, there is an evidence that Z-DNA is stabilized by higher temperature,^73,83–85 which complicates the interpretation of the traditional melting experiments. In addition, traditional melting simulations can be extremely taxing on the computational resources. We ran some tests to probe this. The melting temperature for d(5′-CGCGCGCGCGCG-3′)₂ in B-DNA form at 0.1M NaCl is estimated to be 331 K and 333 K by uMELT⁴⁹ and OligoCalc,¹⁰⁰ respectively. Preliminary tests show that well above 1 μs would be required to obtain the melting products for each temperature around the melting temperature, and because of noise, one needs tens of runs at each temperature interval to obtain a reliable value.

The nonequilibrium results presented here also lend themselves to the interpretation that B-DNA is more stable than Z-DNA. Some quantities are better than others to quantify this difference. Thus, for small fields the WC H-bonds, or the related variable, base-pair opening, might seem to slightly favor Z-DNA. This is because near equilibrium B-DNA frays more than Z-DNA, as shown in the equilibrium (zero field) fraying probability in Fig. S3.¹⁰¹ However, slightly larger values of the field definitely disrupt the H-bonds more in Z-DNA than in B-DNA. Since stacking interactions play such a significant role in polynucleotide stability, the base-stacking index defined in this paper is a better measurement of stability, as it captures the favorable stacking of B-DNA compared to Z-DNA, from zero applied field all the way to the maximum value of the field. Also, the normalized handedness is related to the normalized helicity that some experiments and codes measure. Fig. 7 shows that the application of the laser pulse unwinds Z-DNA (reflected in the loss of absolute handedness) considerably more than B-DNA.

In summary, we have explored a novel, infrared laser pulse melting technique that is suitable for melting polynucleotides is a very fast time. In addition to the intrinsic interest of the melting and healing process, this technique could be used to qualitatively compare relative stabilities of different polynucleotide structures, especially considering that the reduced number of building blocks of nucleotides can readily give common paths for the application of the pulse. The nonequilibrium process does not give such quantities as relative free energies. Rather, it provides one with a “more-or-less” (stable) comparison. This, however, could be extremely useful for stability rankings in competing candidate structures that cannot be easily predicted with other tools, including non-standard secondary structures with noncanonical base pairs and mismatches. In addition, the method applies to very realistic, fully solvated, fully atomistic systems that, in principle, can faithfully reproduce experimental setups.

We thank the NC State HPC Center for extensive computer support. This work has been supported by the National Science Foundation via Grant Nos. SI2-SSE-1148144 and 1534941, and also XSEDE Grant No. TG-MCB150114.