Kinetics of the Conformational Transformation between B- and A-Forms in the Drew-Dickerson Dodecamer.

Some DNA sequences in crystals and in complexes with proteins can exist in the forms intermediate between the B- and A-DNA. Based on this, it was implied that the B-to-A transition for any DNA molecule should go through these intermediate forms also in kinetics. More precisely, the helix parameter Slide has to change first, and the molecule should take the E-form. After that, the Roll parameter changes. In the present work, we simulated the kinetics of the B-A transition in the Drew-Dickerson dodecamer, a known B-philic DNA oligomer. We used the "sugar" coarse-grained model that reproduces ribose flexibility, preserves sequence specificity, employs implicit water and explicit ions, and offers the possibility to vary friction. As the control parameter of the transition, we chose the volume available for a counterion and considered the change from a large to a small volume. In the described system, the B-to-A conformational transformation proved to correspond to a first-order phase transition. The molecule behaves like a small cluster in the region of such a transition, jumping between the A- and B-forms in a wide range of available volumes. The viscosity of the solvent does not affect the midpoint of the transition but only the overall mobility of the system. All helix parameters change synchronously on average, we have not observed the sequence "Slide first, Roll later" in kinetics, and the E-DNA is not a necessary step for the transition between the B- and A-forms in the studied system. So, the existence of the intermediate DNA forms requires specific conditions, shifting the common balance of interactions: certain nucleotide sequence in specific solution or/and the interaction with some protein.

Introduction

The conformational flexibility of the DNA double helix plays an important role in its binding to proteins. Particularly, DNA can locally take the A-form (several riboses switch to C3′-endo, bases slide one relative to another and rotate around their long axes (roll), and the minor groove locally expands).[1] For example, DNA is in the A-form when bound with a cyclic AMP receptor protein,[2] TATA-box binding protein,[3] and some polymerases.[4−6] The transition to the A-form seems to be needed to open and bind the atoms of the minor groove to the non-polar surfaces of proteins.[7] The complex with the lactose repressor proves to be stable only if the DNA molecule is in the A-form.[8] However, the question how the molecule makes the transition from B to A in a cell is open. Therefore, the study of the A–B transition under different conditions and with different control parameters may shed light on the mechanisms of the DNA–protein interaction.

One of the early works on the A–B transition in DNA in vitro (by change of water content in crystalline DNA fibers)[9] reported a strong hysteresis. It means that the transition between the B and A DNA forms in the limit of large molecule length (“in the thermodynamic limit”) is of the first order. On the other hand, some DNA sequences can be crystallized or obtained in DNA–protein complexes in the forms in which some parameters of the DNA double helix take the values as in the A-form, and other parameters—as in the B-form.[10,11] Namely, one can obtain crystals from DNA molecules with the base pair sequence CATGGGCCCATG in which the parameter Slide and the conformations of the riboses are A-like, while the parameter Roll is B-like.[10] The molecules with the sequence (GGCGCC)2 with two brominated or methylated cytosine bases crystallize in a similar shape, called E-DNA (extended along its axis).[12] Moreover, a series of crystals in the intermediate between B- and A-forms was obtained by brominating or methylating the same dodecamer (GGCGCC)2.[13] In this series, the E-DNA is approximately in the middle. The crystalline E-DNA obtained within three weeks transformed to A-DNA after three months.[12] These facts were later referred to as the “Slide first, Roll later” path of the B-to-A transition.

All the DNA molecules detected in the intermediate forms had high content of G:C base pairs. A strong sequence dependence in a B-to-A transition with increasing tetrafluoroethylene (TFE) content was also reported.[14] In the Drew–Dickerson dodecamer (DDD) CGCGAATTCGCG, the G:C and A:T base pairs proved to make the transition at different TFE contents: to have different midpoints (in fact, only the trimeric model, taking into account not the kind of a base pair alone but also two of its neighbors, reproduces experimental data[15]). In addition, the riboses switched to C3′-endo at a lower TFE content than all the bases.

So, the picture of the B-to-A transition in DNA looks as follows. “A-DNA and B-DNA are not two separate, isolated conformations”.[16] Between these forms, there is practically a continuum of intermediate states found in methylated (GGCGCC)2 single crystals.[13] Sugars re-pucker first[14] together with changes in Slide.[16]Roll reacts later (for methylated single crystals, in three months[12]). The transition is strongly sequence dependent in a complicated way:[14,15] base pairs change their stacking at differing TFE contents. If the transition goes as described also in kinetics (as is commonly implied), then it is definitely a continuous, and not a first-order, phase transition in a very heterogeneous system. However, this picture is in direct contradiction with the observed hysteresis in the A–B transition.[9]

Sanyal[17] discussed this problem from the point of view of the Landau theory. The A and B helices belong to different symmetry groups, 111 and 101 (the helices have 11 and 10 base pairs per helix turn) correspondingly. None of these groups is a subgroup of the other, and this is the case when a continuous transition between phases usually takes place. Formally, one can imagine a phase with a higher symmetry (10 × 11)1 (110 base pairs per turn), i.e., almost chaotic, for which both the groups 111 and 101 are subgroups. Then, it follows from the analysis of the Landau free energy that the transitions between this phase and both the A and B phases should be continuous, and the transition between the A and B phases should be of the first order. However, Sanyal has not proposed physically meaningful order and control parameters. In addition, he erroneously interpreted the experimentally observed intermediate DNA forms[10,12] as the coexistence of the A- and B-forms.

Many authors simulated the A–B transition in the framework of the AMBER or CHARMM all-atom force fields,[18−25] but no one addressed the question about the order of the transition. A few authors[26−28] touched upon the question of the sequence “sugars–Slide–Roll” in the kinetics of the B-to-A (or A-to-B) transition in all-atom molecular dynamics (MD) simulations. Orozco et al.[27] stated that in the studied A-to-B transition (in the framework of the B-philic[23] Amber force field), “the transition in RMSd only happens after major changes in sugar puckering”. Pastor[26] found that, in the B-to-A transition “pushed” by A-philic CHARMM force field, Roll finishes the transformation later than Slide, although it is hardly possible to say from Figure 3 of her work which parameter starts the transformation later. Kulkarni and Mukherjee[28] studied the dinucleotide step dependence of the “involvement” of the double helix parameters in the B-to-A transition. Slide, but not sugar pucker, proved to be involved for all the steps, and the general conclusion was that “the local B-to-A transition mechanism is not sequential or stepwise, rather a concerted one involving an orchestra of multiple factors”. Finally, despite all the recent efforts to improve both the AMBER and CHARMM phenomenological all-atom force fields, they still do not reproduce this transition properly (see, for example, Figure 5D by Minhas et al.[29] and Figure 1B by Zhang et al.[25]).

In the present article, we address the two questions, about the order of the A–B transition and about the sequence of the transformation of the DNA double helix. For the calculations, we use the “sugar” coarse-grained (CG) DNA model.[30] It is able to reproduce ribose flexibility with only six beads per nucleotide. The interactions between the bases are analogous to the all-atom ones (which provides correct sequence specificity), but the solvent is modeled through effective potentials for ion–ion and ion–DNA interactions, taking into account solvation effects. The exclusion of the explicit solvent and the small number of beads allowed us to reduce computation time so that we were able to fully study the A–B transition in water with change of available volume.

We choose the DDD as a model system because it is B-philic, and it does not contain G-tracts prone to A-like stacking.[31] The circular dichroism spectrum of the DDD[32] is close to the “reference” spectrum of long natural DNA in the B-form.[33] We study the transition of the DDD between the B- and A-forms induced by decreasing the volume (in water with a small amount of additional salt), as in the first step of the wet spinning method[34] for obtaining DNA fibers.

To determine the order of the B–A transition in the studied system, we build the needed histograms of the order parameter and the shape of the free energy. To investigate the flexibility and the mutual dependence of the helix parameters both in every DNA form and in the process of the transformation, we calculate the corresponding correlation functions. To determine the presence or absence of the E-DNA in kinetics in the process of the transition between the A- and B-forms, we construct a two-dimensional histogram of the distribution of the parameters Slide and Roll on a trajectory with a large number of transitions between the forms.

Results of MD Simulations

Control and Order Parameters of the Transition

As a control parameter, we chose the ratio of the volume available for one Na+ counterion to the volume of the ion. We calculated the parameter v by the formulawhere N is the number of counterions (or nucleotides, 24); V is the volume of the cell, VNa, of a sodium ion, and Vsalt, of the additional salt. We considered ions in the solution as spheres having radii 0.98 Å for Na+ and 1.81 Å for Cl–. These values are Goldschmidt ionic radii,[35] their sum coincides with the distance of the first solvation minimum (direct contact) in the adopted CG potential of interaction between sodium and chlorine ions. The volume of the CG DDD VDNA (5.9027 nm3, the mean between the A- and B-forms) is the volume unavailable for a sodium ion (with the radius 0.98 Å) having the kinetic energy 3k/2 (T is the temperature, 300K). When the cell volume changes from 266 to 26 nm3, the concentration of additional salt changes from 0.1 to 1 M, and the control parameter v, from 2950 to 260.

As an order parameter, we used the average over the molecule width of the major groove W (see Figure ). To calculate W, we took six consecutive base pairs and found the distance between the phosphorus belonging to the first base on the first strand and the phosphorus belonging to the last base on the second strand (counting in the C5′-C3′ direction of the first strand). There are seven such distances in a dodecamer. We regarded their average as the width of the major groove W for the molecule. In the A-form, the major groove narrows when, inside of the molecule, a cavity (filled with counterions) forms.

Figure 1

Histograms of some parameters of the DNA–ions system at three values of volume available for a counterion: before (v = 1470), near (v = 635), and after (v = 427) the transition between B- and A-DNA (NVE ensemble, inviscid water). The values of the order parameter W (the averaged over the molecule major groove width) are plotted against the frequency of the values. EP – Na is the interaction energy between the DNA phosphate beads and the sodium ions, and Na is the number of sodium ions at a distance less than 8 Å from the helix axis. We also traced (averaged over the molecule) values |C1′P| (indicating the sugar conformation), Roll, Slide, and Twist (the standard DNA geometric parameters). The total number of points in one histogram is 20,000. The time between counts is 0.1 ns.

In addition, we traced several more parameters characterizing the molecule as a whole (see Figure ). The first is the (averaged over the molecule) length of the bond |C1′P| with imposed double-well potential reproducing ribose flexibility. The sugar ring conformations C2′-endo and C3′-endo correspond to the left and right wells, respectively, in the interaction potential between the C1′ and P beads (see the description of the sugar model in Section and Figure 3 in the article[30] about the model). The second feature parameter, which has different values in A- and B-DNA, is the interaction energy of P beads (phosphates) with sodium ions EP – Na. In A-DNA, it is lower because counterions spend more time in the cavity inside the major groove, next to phosphates and negative partial charges on the bases. Finally, we monitored the standard geometrical parameters of the helix averaged over the molecule Roll, Slide, and Twist.[36]

Variation of the Viscosity in the Model

In contrast to all atom models, the sugar model allows to change the friction in the system. For both the DNA beads and the ions, we use the Langevin thermostat, which naturally introduces the viscosity of water into the model. We can separately vary the friction coefficients for DNA beads and for ions. In particular, we can consider the idealized case of “inviscid” water with zero friction (NVE ensemble). In this case, the inertial properties of the DNA molecule determine its dynamics.[37]

We analyzed MD trajectories 2 μs long for 13 volumes v in a range from 260 to 2950. Besides the inviscid water, we considered the four following combinations of friction coefficients for DNA beads and ions: (γDNA = 5 ps–1, γions = 5 ps–1), (γDNA = 5 ps–1, γions = 70 ps–1), (γDNA = 50 ps–1, γions = 5 ps–1), and (γDNA = 50 ps–1, γions = 70 ps–1). The last combination corresponds to the friction in real water. For the last three combinations (with large friction), we carried out additional simulations near the transition point (v = 1060 – 720), totaling 20–30 μs.

First Order of the A–B Transition

Our simulations[38] showed that, in a wide range of variation of the control parameter, the DDD (forming one turn of the double helix) behaves like a small atomic cluster in a transformation, which, in the thermodynamic limit, corresponds to a first-order phase transition. At a first-order phase transition, large systems demonstrate a coexistence of two phases. Small atomic clusters in the transition region spend one part of their time in one phase and another part in the other phase.[39] Similarly, a single turn of a free DNA helix cannot be partly in A-form and partly in B-form. Thus, due to thermal fluctuations, the DNA molecule incessantly changes one form into another, which is well seen in the time dependence of the (average for the molecule) major groove width (Figure ).

Figure 2

Time dependence of the (average for the molecule) major groove width W at different values of volume available for a counterion (v = 1470 before the transition from B- to A-DNA, v = 720 near it, and v = 427 after the transition) and for different values of friction coefficients for the DNA molecule γDNA and for the ions γions (γDNA = 50 ps–1 and γions = 70 ps–1 correspond to the real water). The value of W for the B-DNA is about 25 Å, and for A-DNA, it is 10 Å. In the bottom row of the graphs, we show the autocorrelation functions for the trajectories at the volume 720.

We have carefully checked (see Figure ) that there were no intermediate forms, like E-DNA (with a wide major groove (small Roll) like in the B-form and large negative Slide like in the A-form), at the used conditions including the transition area. Therefore, one order parameter (reaction coordinate) W is sufficient for the studied transition.

Figure 3

Histogram of the (average for the molecule) Slide and Roll parameters on a trajectory of 100 ns long (step 1 ps) calculated in the NVE ensemble (“inviscid” water) for the available volume v = 635 (near the transition). We show two portions of the trajectory with transitions from A to B (length 25 ps, black points) and from B to A (length 28 ps, magenta points). The E-DNA region (Roll is close to zero, with large negative Slide, lower left corner) is not populated.

The histogram for W appeared to be bimodal (see Figure ). We define the transition point as the volume at which the areas under the two peaks in the histogram are equal. This point, v= 635, corresponds to the transition also for the parameter EP – Na which also characterizes the state of the molecule as a whole. On the histograms for the parameters Slide, Twist, and Na (see Figure ), the area between the peaks at the transition point is much more populated just because these parameters have large dispersion in both the A- and B-forms. The fraction of riboses in unfavorable conformation (tracked by parameter |C1′P|) at a small (v = 425) as well as at large (v = 1470) volume is rather large, which well agrees with experiment.[16]

Besides the good sampling in all regions of the phase space within available to our CG model simulation times of a few tens of microseconds, we also carried out a more accurate calculation (using the metadynamics method[40,41]) of the free energy profile depending on the reaction coordinate W for different values of the control parameter (see Figure ) to find the quantitative characteristics of the transition. One can see that the observed change in the free energy shape is classical for a first-order phase transition.

Figure 4

Analysis of the free energy profile for the Drew–Dickerson dodecamer in the framework of the sugar CG model (metadynamics method). The reaction coordinate is the width of the major groove W (average for the molecule). The control parameter is the (normalized) volume available for a counterion v. Left: the free energy as a function of the reaction coordinate. Right: the difference between the local minima for the B- and A-forms and the barriers for the transitions from A- to B- and from B- to A-form when the control parameter changes. For every v, the length of the trajectories was 31.2 μs (the total for eight trajectories). Friction coefficients for the DNA and the ions: γDNA = 5 ps–1 and γions = 5 ps–1.

A–B Transition Parameters Depending on Viscosity

The dodecamer jumps between the forms not only in the NVE ensemble, without friction, but also at any considered combination of the friction coefficients for the DNA beads and the ions. If we define the major groove width value WAB separating the A- and B-forms as 14.5 Å (the position of the minimum between the two peaks on the histogram in Figure ), then we can plot the dependence of the A-form fraction (the area below the left peak) on the volume. These curves were obtained for all the friction coefficients considered (Figure ). It turned out that the curves coincide within the measurement accuracy. The only result of the water viscosity is a substantial slowdown of the dynamics, the time required for the helix transformation from one form to another and the time between the jumps increasing in equal proportion (see Table ). Namely, as compared to inertial dynamics of the molecule in inviscid water, these time intervals for the dynamics in the real water increase by two orders of magnitude, the time of transformation up to 5 ns, and the time between the jumps up to 80 ns.

Figure 5

Fraction of the A-DNA (the share of time the DNA spends in the shape with the average width of the major groove less than 14.5 Å) depending on the (normalized) volume available for a counterion v at different friction coefficients for the DNA and the ions. In the NVE ensemble, water is inviscid. The label 50/70 means friction coefficients γDNA = 50 ps–1 and γions = 70 s–1, and this combination corresponds to the real water. The dashed curves mark the measurement error.

Table 1

Number N of Jumps of the DNA Molecule between the Forms and the Time of the Transformation τ for Different Pairs of Friction Coefficients for the DNA and the Ionsa

Friction	N	τA→B (ns)	τB→A (ns)
γDNA = 0, γions = 0	2590	0.06 ± 0.03	0.07 ± 0.04
γDNA = 5 ps–1, γions = 5 ps–1	211	0.9 ± 0.5	1.0 ± 0.5
γDNA = 5 ps–1, γions = 70 ps–1	79	2.6 ± 0.9	2.3 ± 0.9
γDNA = 50 ps–1, γions = 5 ps–1	28	4 ± 2	5 ± 2
γDNA = 50 ps–1, γions = 70 ps–1	24	5 ± 2	5 ± 2

The calculations were performed for the volume v = 720 (near the transition) on trajectories 2 μs long. The time τ was calculated from the change of the width of the major groove (average for the molecule).

Motion in Configuration Space Depending on Viscosity

The comparison of the MD trajectories with different friction coefficients (Figure ) allows to deduce that it is the DNA molecule, and not the counterions, that initiates the jumps between the forms. Indeed, for the same viscosity for the ions, at the small viscosity for the DNA beads, there are more jumps (per time unit) between the forms than at the large viscosity for the DNA beads. The reason for that is, evidently, the presence of the two locally stable configurations for the A- and B-forms.

The larger ion viscosity leads not directly to a decreased number of jumps but rather to a delay in the current form, even if it is unfavorable for the given volume v. Fast ions, instead, quickly return the DNA molecule into the form with lower free energy at the given volume. Interestingly, on average, only four counterions more (15 vs 11 out of 24) reside near the DNA (at a distance less than 8 Å from the axis of the molecule) in the A-form than in the B-form. In this connection, we should note that, in real water, ions are the slowest to equilibrate among the parts of the DNA–ion conglomerate (apart from the water). Therefore, the DNA molecule may stay for a long time in a non-equilibrium state, for example, when binding to proteins.

DNA Flexibility in A- and B-Forms and in the Process of Transformation

For such a small molecule in water at 300 K, even the average helix parameters (not to mention local characteristics) have very large dispersion in both the A- and B-forms, as well as in the process of transformation. Therefore, it is practically impossible to precisely define the time at which a transformation starts, passes the landmark between the forms, and is completed. The error is large for all the parameters, even for the best one, W (see Figure ), and unacceptable (as compared with the transformation time) for Slide, Twist, and Na (see Figure ). For some individual transformations (not for all), one can say that one of the parameters starts to change (or passes the landmark or completes the transition) first, then goes the second parameter, and then the third.[38] However, this sequence changes from transformation to transformation, and there is no permanent “leader”, at least for the dodecamer as a whole. As a result, it is impossible to reliably say which parameter makes the transition first and which does later.

However, one can extract the information about the character of flexibility of the DNA double helix from autocorrelation and correlation functions of the helix parameters. We obtained two kinds of these functions. The “full” functions were calculated on the full trajectories with many jumps between the forms. The “short” functions characterize only the mobility of the A- and B-forms separately. The short functions were calculated on the pieces of the trajectories between the jumps. The results for the friction corresponding to the real water are presented in Figure . The characteristic times for the functions are listed in Tables and 3. For the other friction coefficients, only the time scale proved to be shorter (for example, 80 times shorter for the inviscid water). The shapes of the correlation functions are very similar.

Figure 6

Analysis of the DNA flexibility at volume v = 720 (near the transition point) for the case of the real water (γDNA = 50 ps–1, γions = 70 ps–1). The full length of the trajectory used in the calculations is given above every plot. For the A- and B-forms, it consisted of 20 pieces chosen between the jumps from one form to another. The time between the used points on the trajectories was 0.1 ns.

Table 2

Time Intervals Characterizing the Flexibility of the DNA Double Helix (“Sugar” CG Model) for the Cases of Inviscid (NVE) and Real (γDNA = 50 ps–1, γions = 70 ps–1) Watera

τ	inviscid	real
first zero of full autocorrelation function	2 ns	∼70 ns
interval between jumps	0.77 ns	83 ns
duration of transformation	0.06 ns	5 ns
decay of (auto)correlation in the forms between jumps	0.025–0.03 ns	2–3 ns

The listed times are close for all the helix parameters.

Table 3

Maximal Absolute Values of Full (for Transitions) and Short (for the A- and B-Forms Separately) Correlation Functions for Some Helix Parameters (in the “Sugar” CG Model) for the Case of Real (γDNA = 50 ps–1, γions = 70 ps–1) Watera

correlation	transitions	A	B
C1′P–Roll	0.82	0.4	0.2
C1′P–Slide	0.82 → 0.55	0.7	0.68
C1′P–Twist	0.7 → 0.45	0.48	0.68
Slide–Twist	0.75 → 0.4	0.45	0.65
Slide–Roll	0.6	0.25	0.2

The second value after the arrow is the level of full (long) correlation after the time of short correlation. We listed these values only for the pairs of parameters highly correlated in the forms.

First of all, the kinetics of the forms and of the transformations is similar for the inviscid and real water, with the time scale of the former approximately 80–100 times shorter than of the latter. The time of decay of the full autocorrelation functions is of the order of time between the jumps (more exactly, it is three times more for the inviscid water and the same for the real water, but the last estimate is not reliable because of insufficient statistics (a single trajectory 2 μs long)). The autocorrelation functions for the A- and B-forms drop much sooner, the time of their decay is about two times less than the time of the A–B transformation and two orders less than the loss of the full autocorrelation (Table ).

The short autocorrelation functions for C1′P (sugars) and Slide fade out perfectly synchronously and at the same rate as their short correlation function, which starts from the value 0.7. This means that there is really high correlation of these parameters, evidently through valent and torsion angles connecting sugars and bases. The mechanical explanation of this known experimental fact has been given by Dickerson and Ng.[16] On the full trajectory, C1′P–Slide correlation drops from 0.82 down to 0.55 exactly after the time of their correlation in the forms. The similar situation is with sugars and Twist. Twist has also proved to be tightly connected with the sugar conformation. The parameters Slide and Twist seem less than W and Roll connected with the jumps between the forms: the full autocorrelation functions of Slide and Twist drop sooner than the ones of W and Roll. However, the reason of this is that their values in the forms have large dispersion and are much worse separated (see Figure ).

The situation with sugars and Roll is opposite to the situation with sugars and Slide. In both the forms, Roll (and W) is correlated with C1′P very weakly and rather not at all in B-DNA. These helix parameters are mechanically independent. However, the correlation on the full trajectory is really high and long, indicating that, in the transformations, they change coherently.

Not unexpectedly, the pair Slide–Roll is not correlated in the forms while the pair Slide–Twist is. On the full trajectory with jumps, these two correlation functions are close, except for the initial part of the curves, which reflects the presence or the absence of the correlation in the forms. The correlation function Slide–Roll has no peak at nonzero time, and, correspondingly, there is no the “Slide first, Roll later” regularity for the molecule as a whole. The kinetics of the transformations is not stepwise, the helix parameters change synchronously on average (see Figure ).

Figure 7

Kinetics of changes in the shape of the double helix in the transition region between B- and A-forms in the Drew–Dickerson dodecamer. The geometry of the helix is such that the parameters Slide and Roll are almost uncorrelated in either the A- or B-forms. Therefore, in principle, an E-form is possible. However, on the full trajectory in the transition region with a large number of jumps between the forms, these parameters do noticeably correlate (Figure ), with a maximum at zero delay. This means that, during the transformation of the helix, they change, on average, synchronously. The E-DNA is absent on the trajectory in the transition area (see Figure ).

Discussion: Time of the B- to A-Transition in the Experiment on Long Chains

Although it is aside from the main problem with which the article is concerned, it seems inevitable to discuss the time of transition from the B- to the A-form experimentally obtained by Jose and Porschke.[42−45] In the framework of our CG model, in the area of the transition, the time of the transformation is 5 ± 2 ns, the same for both the A-to-B and the B-to-A transformations. This time interval is close to the estimates made for the transitions in water in the framework of both the CHARMM[26,46] (B to A) and the AMBER[47] all-atom force fields. We think that, for short DNA oligomers, this result corresponds to reality, including the same time for both the transformations, because this time depends mainly on the (mechanical) properties of the molecule itself and, to a lesser extent, on the properties of the solvent.

However, Porschke experimentally estimated (by an electric field jump technique) the time of the B-to-A transition for poly(AT)[42] in an ethanol–water mixture as 2.5–25 μs (depending on the percent of ethanol; the same time for the molecules consisting of 70 and of 1600 base pairs within experimental accuracy ±10%). For natural DNA,[43] the main part of the magic amplitude (75%) in the center of the transition was associated with time constants from 1 to 10 μs (depending on the GC content) for the lengths between 2500 and 8000 base pairs. The second large component of the magic amplitude appeared with time constants of 50–100 μs. We see that even the shortest estimates are three orders of magnitude longer than the time of transformation obtained in MD simulations. Pörchke also tried[48] to obtain the time of stacking in a model compound with two adenine residues connected by a −(CH2)3– bridge. The time could not be resolved, and so it was less than 10 ns. On the other hand, the stacking rearrangement in an RNA hairpin[49] is on the same time scale of microseconds.

In 2007, Orozco et al[27] explained the difference between the MD results and the Porschke’s experiments by a huge barrier (11.4 kcal/mol) in the potential of mean force on the way from the B-form to the A-form appearing in an ethanol–water mixture. However, Orozco et al. used the B-philic AMBER force field, which, as is now known, does not reproduce the A-minimum in ethanol correctly.[25] However, we do not think that the currently used all-atom force fields, although not quite balanced, do not take into account something that can change the time of the B–A transition by three orders after changing water to the ethanol–water mixture. If we also discard the hypothesis that the time of the B–A transition for long molecules in an ethanol–water mixture is less than 10 ns, and so cannot be resolved, then we should examine the experiments[42,43] more closely.

A short pulse of a constant electric field was applied to a DNA solution. Under the action of this field, DNA molecules (poly(AT)[42]) polarized and (partially) oriented along the direction of the field. After pulse termination, the anisotropy of the solution relaxed, which was recorded by optical methods. The experiment was carried out in a mixture of water with ethanol, in the range of concentrations including the B–A transition. According to the circular dichroism spectrum, the B-to-A transition took place between 65.87 and 72.97% of ethanol, the midpoint being at approximately 69.13% (Figure 2 in ref (42)). In this experiment, if one observes the relaxation process in the light polarized at the magic angle (54.7°) to the direction of the pulsed electric field, then the rotational relaxation is not visible. If any relaxation is observed at this angle, then it should correspond to another process (“reaction”) that simultaneously occurs in the system. At the end of the B–A transition, in the area of the predominant A-form, between 68 and 74% of ethanol and with the midpoint at 71%, another reaction was observed with differing rise and relaxation times (these time intervals were of the same order of magnitude as the rotational ones, they might be longer or shorter depending on the chain length[43]). Porschke ruled out all the possible reactions, such as bending of molecules, their aggregation, and denaturation. Therefore, he interpreted the observed reaction as the transition of DNA from the A- to the B-form under the action of the field and the subsequent transition back to the A-form after turning off the field. Correspondingly, Porschke identified the magic rise time with the time of the A-to-B transition and the magic relaxation time with the time of the B-to-A transition. The length of the chains in the described experiments varied from 70 to several thousands of base pairs. The magic time constants weakly depended on the chain length (Figure 4 in ref (43)), contrary to rotational diffusion time constants.[44] The co-operativity length of the B–A transition was estimated a long time ago, it lies in the range from 10 to 30 base pairs.[33] It coincides with a couple of helix turns. So, one can conclude that the transition on the quite short fragment belonging to a long (70–8000 base pairs) DNA molecule takes several microseconds. The question is, why so long?

The process of structural transformations in polymer chains, quasi-one-dimensional systems, differs from a phase transition in low-molecular compounds like water. We saw that a free short duplex could not be partially in the A-form and partially in the B-form, it jumped between the states in the area of transition. The reason is that the boundary between the forms should be highly unfavorable, energetically. Since the 1980s, it was common to interpret the boundary between the A- and B-forms as a “soliton” or “kink”, a topological localized defect moving from one end of the duplex to another and changing the state of the molecule in the process (see a short review in chapter 9.1 of the book[50] by Yakushevich).

However, due to the high viscosity of water and due to the mobility of the DNA molecule, such defects cannot move along the chain at a high rate. Since the transition time is practically independent on the chain length, large displacements of the boundaries between the forms can be excluded, as well as the propagation of a new form only from the ends of the chain. Therefore, in the process of the transition, regions of a new form (with two boundary “solitons”) should appear independently in different parts of the chain and propagate over only short distances, a couple of helix turns. The formation of such a nucleus on a long chain requires overcoming a large energy barrier (to create two boundaries between the forms) and, therefore, should take much longer than the transition on a short oligomer with free ends. This hypothesis seems to be the most plausible explanation for the three orders of magnitude difference in the transition time between simulations on short chains and the experiment on long chains, although the detailed understanding of the process requires further theoretical and experimental study.

Conclusions

In the present work, we studied the process of conformational transformation between the B- and A-forms in a B-philic oligomer, the Dickerson–Drew dodecamer. We chose the (normalized) volume available to a counterion as a control parameter and decreased the volume of the cell from 266 to 26 nm3 with a small amount of additional salt. We used the sugar CG model[30] on time intervals of tens of microseconds. For large volumes, the DDD is in the B-form. At small volumes, the molecule takes the A-form. Under the described conditions, the transition that we observed corresponds to a first-order phase transition in the limit of large molecule lengths.

Namely, the DNA dodecamer behaves like a small cluster, jumping from one form to another in a wide range of volumes. The control parameter determines the residence time of the molecule in one or another form. Using the metadynamics method, we obtained the change in the shape of the free energy with a change in the control parameter. It proved to be classical for a first-order transition. Finally, for a long trajectory with a large number of jumps between the A- and B-forms, we have built a two-dimensional histogram of the distribution of the Slide and Roll parameters. The area of the “intermediate” E-DNA is not populated.

We studied in detail the flexibility of the DNA molecule in both the A- and B-forms and in the process of the jumps. We calculated the autocorrelation and correlation functions of the (average for the molecule) helix parameters and sugar conformation (length of the double-well CG bond C1′P). We analyzed both “short” (on the pieces of the trajectories between the jumps) and full correlation functions.

The possibility of varying friction in the used model separately for the DNA and for the ions allowed us to determine that the “driving force” of the transformation is the DNA molecule with its two local minima of energy. The ions determine the character of the system behavior in the configuration space. With a large friction of the ions, the DNA molecule tends to linger in the current, even unfavorable, conformation. The “fast” ions quickly return the molecule to the global minimum of free energy. The introduction of the friction corresponding to real water slows down the dynamics of the system by about two orders of magnitude as compared to “inviscid” water. There is no influence on the parameters of the transition, such as its midpoint. In our model with real friction (γDNA = 50ps–1, γions = 70 ps–1), the characteristic decay time of the (auto)correlation functions in the A- and B-forms is about 2 ns (the B-DNA is softer, and the correlations decay longer in it), the time of the transformation between the forms is 5 ns, the time between the jumps in the transition region (v = 720) is 83 ns, and the time after which the full autocorrelation function reaches the first zero is of the same order.

It turned out that, indeed, as follows from mechanical considerations,[46] for each of the forms, especially for the B-DNA, the correlation of C1′P and Slide is very large. The correlation of C1′P and Twist is slightly weaker (only in the A-form), which also follows from the analysis of the geometry of the double helix.[51] On the other hand, the Roll parameter does not correlate with either C1′P or Slide in the B-form. For the stiffer A-DNA, a slight Roll–C1′P correlation is present. So, the Slide and Roll parameters in the DNA double helix are really almost not interdependent. It is this feature of the geometry of this helix that makes possible an intermediate E-form between the A- and B-DNA: the form with a large magnitude of Slide, like in A, and almost zero value of Roll, like in B.

However, in the B-philic DDD, in water, and with a small amount of additional salt, the parameters of the helix change synchronously in average in the process of transformation between the forms. Therefore, the full correlation function of Slide and Roll is quite large (starts at 0.6), although less than, for example, the correlation between C1′P and Roll (starts at 0.82). The difference between these two functions results from a large overlap of the possible values of the Slide parameter in the A- and B-forms. We did not observe the sequence “Slide first, Roll later” in kinetics. Also, none of the helix parameters passed the landmark between the forms earlier than the others on average (we did not study the local properties depending on the type of base pairs).

Thus, the transition conditions used by us do not cause the appearance of E-like forms between A and B. Here, we should note that we have observed the E-DNA in the framework of our model under differing conditions. So, the E-DNA or other forms of the DNA molecule intermediate between A and B, known in DNA crystals[10,12,13] or in complexes with proteins,[11] are not necessary steps that the molecule goes through when it passes from B to A or vice versa.

For the existence of the intermediate between the A- and B-forms of DNA, special conditions shifting the balance of interactions are necessary. In particular, they can appear in a cell when a DNA molecule interacts with a protein in a solution of a certain composition. Additionally, these conditions may be satisfied only selectively for some sequences of base pairs. DNA crystals in intermediate forms[10,12,13] were obtained for oligomers with a large portion of G:C base pairs at a low temperature (4 °C)[10] or after the methylation or bromination of several bases.[12,13] In addition, the DNA molecules crystallized in the solution having very complex composition with a large fraction of different molecules (spermine, Ca acetate, 2-methyl-2,4-pentanediol, cobalt hexamine, and MgCl2 to name a few). All these factors (especially base modification in works by Ho et al.[12,13]) were considered only as the means to “catch” DNA molecules in the intermediate states. However, they rather played a decisive role in the appearance of these unusual helices.

Methods

“Sugar” Coarse-Grained DNA Model

The conformation of a DNA molecule crucially depends on external conditions, including the volume accessible to the molecule, the solvent composition, surrounding molecules, both small, like ions,[34] and large, like DNA itself,[24] or other polymers[52] including proteins.[11]

An isolated DNA with sodium counterions (NaDNA) in a large volume of physiological saline takes the B-form. The most well-known factor causing the B-to-A transition is water reduction. It can be achieved by simple drying of living cells[53,54] or by adding ethanol (or TFE) to DNA in water.[33] It is less known that an increase in the concentration of the additional NaCl salt at a low concentration of DNA in solution also leads to the transition of some DNA oligomers to the A-form. The other oligomers may keep the initial B-form or change to the Z-form (see the review of experimental results in the work[55] by Peticolas et al.). Likewise, DNA oligomers can crystallize in one of the A-, B-, C-, or Z-forms depending on the nucleotide sequence.[55] However, the fibers of long NaDNA molecules with a random nucleotide sequence are in the A-form under normal relative humidity.[56] An all-atom MD simulation of a DNA molecule in a small water drop[47] showed that the B-to-A transition is accompanied by collecting a number of counterions (together with water) in the major groove of the DNA. These ions constrict the major groove due to electrostatic interactions with negative charges on DNA.

All the described features of the interaction of DNA with its environment shifting the balance to the B- or A-form are taken into account in the “sugar” CG model of DNA.[30,38] The ions are modeled explicitly. The solvent is excluded, but the interactions between ions, phosphate beads, and charges on bases are described by CG potentials, taking into account the effects of solvation. This makes it possible to simulate the change in the properties of the solvent, like adding alcohol or drying. To keep the sequence specificity, we model the interactions between the bases by the all-atom AMBER potentials (parm99,[57] the parameters of interactions between the bases have not been changed since then). In the model, both the A- and B-geometric forms of the helix are local minima of potential energy with a barrier between them. This is achieved (besides the all-atom minima in the interactions between the base pairs) by introducing simplified ribose conformational flexibility (see Figure ). More exactly, we use one double-well potential for the (short) CG bond between C1′ and P(2) and one three-particle potential providing correlation between the distances |P(1)P(2)| and |C1′P(2)|. The model uses the Langevin thermostat, which naturally introduces water viscosity.

Figure 8

Introduction of the conformational flexibility of the ribose ring into the CG sugar model.

The sugar model was built to study the A–B transition. However, the structure of the model is universal, and minor modifications make it possible to use it for modeling DNA dynamics at almost any conditions. Namely, in the current version, only the basic conformational mobility of the DNA backbone is realized: the flexibility of sugars. BI–BII (ε/ζ) and α/γ mobility can be taken into account in a similar manner. Also, in this version, sequence dependence is not fully introduced (it is absent in the base rotation potential), and the effective potentials of the interaction between ions and charges on DNA need tuning. In addition, the backbone is a little biased to the A-DNA.[30]

Details of MD Simulations

We performed an MD simulation of the Drew–Dickerson dodecamer (12 base pairs, and all the 24 phosphate beads were present) with counterions (24 Na+ ions) and 32 additional salt ions: 16 Na+ and 16 Cl–. In our model, on the time intervals used for calculations, there was no base pair fraying or other “end effects” overly present in most all-atom MD simulations.[58] We suppose that the reason is that we do not have explicit water molecules facilitating opening events, and we have not included the corresponding part into our CG potentials. Therefore, in the calculations of the model parameters (W, averaged Slide, Roll, and others), there was no need to omit terminal base pairs.

To reproduce the B-to-A phase transition, we reduced the volume of the computational cell from 266 to 26 nm3 at a constant amount of additional salt (actually repeating the first step of the wet spinning method[34] for obtaining DNA fibers). The cell was a cylinder with a height of 40–64 Å and a diameter of 29–72 Å.

To obtain the results of the present article, we implemented the sugar CG DNA model[30] into the LAMMPS[59] MD package[60] with slight differences from the program used to obtain the results of the previous articles.[30,37,38] First of all, the bases were replaced with rigid bodies (“fix rigid” command) having inertia tensors equal to the ones of the all-atom representations of the bases. We always used the reflecting boundary conditions, which is closer to the situation for a DNA in a cell than the periodic boundary conditions. Due to the impossibility of using hard reflective walls for rigid bodies, we employed, for both the DNA grains and the ions, soft walls with the Lennard-Jones (6–12) interaction potential. Its coefficients were ε = 0.2 kcal/mol, Rmin = 5 Å. At r = Rmin, the potential was cut. In calculating the cell volume, we accepted that at temperature 300 K, the mass centers of the particles can approach the wall not closer than 4.32 Å. In addition, the maximum value of the angle C1′C3′P(2) was, for technical reasons, limited to 180°. To the very weak torsion potential for the angle C3′C1′NC, we added a soft term preventing free rotation of the bases by 90°. Optimization in the LAMMPS package allowed, for the DNA dodecamer at a salt concentration of 0.5 M, calculating the trajectory at a rate of 14.5 ns per 1 h on one core of the Skylake unit at the Joint Supercomputer Center of the Russian Academy of Sciences. The LAMMPS’ parallelization in space reduces computation time for high salt concentrations by 25% when using two cores, and further segmentation for such a small system is ineffective. Invoking LAMMPS in multi-partition mode allows obtaining trajectories for all the cases (for the system at different viscosities and salt concentrations) in parallel.

For the metadynamics[40] calculations, we used the Colvars module[41] in LAMMPS. We carried out the calculation for the system γDNA = 5 ps–1, γions = 5 ps–1. To provide the fulfilment of the condition τζ≪100 ns (the longest correlation time τζ is about 10 ns at the specified friction; see Figure ), we chose the following values for the parameters. The hill height was V = 5.962·10–5 kcal/mol (kBT/10,000), width (2·σ) = 0.5 Å, δt = 10 ps. To keep the reaction coordinate W within the grid, we used Colvars’ harmonic walls with lower boundary = 5.5 Å, upper boundary = 28.5 Å (the area much wider than the range observed in MD simulations), and force constant = 10 kcal/mol. The calculations were carried out until stabilization of the resulting mean force potential in the needed area was reached.