Unwinding Induced Melting of Double-Stranded DNA Studied by Free

Subscriber access provided by READING UNIV

Feature Article

Unwinding Induced Melting of Double-Stranded DNA Studied by Free Energy Simulations Korbinian Liebl, and Martin Zacharias J. Phys. Chem. B, Just Accepted Manuscript • DOI: 10.1021/acs.jpcb.7b07701 • Publication Date (Web): 24 Oct 2017 Downloaded from http://pubs.acs.org on November 3, 2017

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

The Journal of Physical Chemistry B is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

The Journal of Physical Chemistry

Unwinding Induced Melting of Double-Stranded DNA Studied by Free Energy Simulations Korbinian Liebl and Martin Zacharias∗ Physik-Department T38, Technische Universität M¨ unchen, James-Franck-Str. 1, 85748 Garching, Germany E-mail: [email protected] Abstract DNA unwinding plays a major role in many biological processes such as replication, transcription and repair. It can lead to local melting and strand separation and can serve as a key mechanism to promote access to the separate strands of a double stranded DNA. While DNA unwinding has been investigated extensively by DNA cyclisation and single-molecule studies on a length-scale of kilo base pairs, it is neither fully understood at the base pair level nor at the level of molecular interactions. By employing a torque acting on the termini of DNA oligonucleotides during Molecular Dynamics free energy simulations we locally unwind the central part of a DNA beyond an elastic (harmonic) regime. The simulations reproduce experimental results on the twist elasticity in the harmonic regime (characterized by a mostly quadratic free energy change with respect to changes in twist) and a deformation up to 7◦ was found as a limit of the harmonic response. Beyond this limit the free energy increase per twist change dropped dramatically coupled to local base pair disruptions and significant deformation of the nucleic acid backbone structure. Restriction of the DNA bending flexibility resulted in a stiffer harmonic response and an earlier onset of the anharmonic response. Whereas local melting with a complete disruption of base pairing

1 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

and flipping of nucleotides was observed in case of an AT rich central segment strong backbone changes and changes in the stacking arrangements were observed in case of a GC rich segment. Unrestrained MD simulations starting from locally melted DNA reformed regular B-DNA after 50-300 ns simulation time. The simulations may have important implications for understanding DNA recognition processes coupled with significant structural alterations.

1

Introduction

DNA deformability and in particular the twist deformability is of fundamental importance for vital genetic processes such as replication, transcription, DNA repair and protein-DNA interactions 1–5 . In vivo DNA is mostly found in circularly closed form including supercoiling or restricted by packing in complexes with proteins causing a mechanical stress (torque) on the DNA that can strongly affect recognition by proteins and influence transcription and replication fidelity 6,7 . The total twist of a double-stranded (ds)DNA and to some extent also its twist elasticity can be obtained by DNA ring closure ligation (cyclisation) experiments 8–11 . The technique has been used to study the twist flexibility of DNA fragments of various lengths and sequences in solution 12 . Time-resolved fluorescence polarization anisotropy (on DNA with intercalated ethidium bromide) can also be used to measure torsional DNA flexibility near the equilibrium state 13,14 . These studies indicate that the overall base composition has only a modest influence on the average twist flexibility of DNA. However, external stress on DNA such as DNA bending or stretching can influence the effective twist elasticity 8,15–18 . Large scale unrestrained Molecular Dynamics (MD) simulations have been performed to investigate the sequence dependent flexibility of DNA 19,20 . These studies and other longtimescale equilibrium simulations indicate that the local twist at individual base pair steps can follow a Gaussian distribution but frequently shows deviations from the Gaussian (harmonic) case or even bimodal behavior due to meta stable substates 21–25 . However, when considering the average (global) twist over a DNA segment these distributions average out to a 2 ACS Paragon Plus Environment

Page 2 of 40

Page 3 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


mostly harmonic deformability with respect to twist and other helical parameters (Gaussian distribution around the mean twist) 26–28 . Hence, in many cases a harmonic stiffness model seems to be sufficient to characterize global DNA flexibility at equilibrium 28–31 . However, many biological processes such as strong DNA supercoiling, transcription and replication, or DNA packing can lead to twist deformations beyond small equilibrium fluctuations. This can result in an anharmonic response and may also cause significant structural changes in DNA such as melting of the double strand 32,33 . For instance, Jeon et al. 32 found that supercoiling can induce melting bubbles in DNA and Kim et al. 2 suggested that torque induced unwinding underlies the effective basic mechanism driving promoter melting during transcription. Another striking example relates to the repair of damaged DNA through the NER (nucleotide excision repair) complex, where Velmurugu et al. have shown that lesions in the DNA are recognized by probing the DNA’s local unwinding deformability 3 . Simulation studies on deformed DNA also indicate that unwinding of DNA facilitates base pair opening and base flipping as part of the damage recognition and repair process. 34,35 Experimentally, the mechanical response to reducing the winding number between DNA’s both strands has been studied on a length scale of several kilo base pairs (kbp) 4,11,17,36–39 . In this context, not only global material characteristics of DNA have been determined, but also DNA supercoiling has been investigated in detail. A recent cryo-ET, biochemical and Molecular Dynamics (MD) simulation combined study has even detected localized denaturation in negatively supercoiled DNA 33 . However, the molecular mechanism of the release mechanism triggering localized denaturation of underwound DNA has still remained elusive.

Recently, it has also become possible to study the elastic response of single DNA molecules upon application of magnetic torque tweezers to un- or over-twist DNA 17,38–41 . Such experiments provide valuable insights into the physical and elastic properties of DNA. However, only relatively long DNA molecules (∼ 1 kbp) can be studied that can react on the external stress by local as well as global conformational changes (e.g. supercoiling that relaxes part



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

of the “twist stress”). Similar to the DNA ring closure experiments an overall twist elasticity can be obtained, however, with only modest insight into the conformational changes taking place at the molecular level. Molecular dynamics simulations on DNA have matured in recent years due to improved force field descriptions and improvements in the computer hardware resulting in DNA dynamics and deformability in close agreement with experiment for regular B-DNA 19,20,23,28 . This includes also the prediction of global flexibility of DNA for bending, twisting and stretching in agreement with experimental data 30,42 . The global flexibility variables can also undergo coupled changes, e.g. changes in twist can correlate with changes in bending or stretching. Recently, the molecular mechanism of opposite twist-stretch coupling in DNA 26 and RNA found experimentally in single molecule torque experiments 38 could be explained based on restrained and unrestrained MD simulations 43 . While equilibrium MD simulations and simple models based on treating DNA as a elastic material (e.g. twistable worm like chain model 16,44 ) and based entirely on empirical parameters seem to cover the harmonic regime of DNA deformability quite well, twist deformations beyond the harmonic response have so far not been investigated systematically using atomistic simulations. In the present MD study, we have extended an earlier approach of inducing changes in twist to DNA by applying an external torque 27,43 . Other previous efforts in this direction involved either very short MD simulations 45 or special setups that are based on fitting a DNA into a periodic box which may significantly restrict relaxation motions of the DNA 46 . In our approach a biasing potential causing a torque on the ends of DNA oligonucleotides is applied to induce unwinding in small steps combined with Hamilitonian replica exchange MD simulations (H-REMD) along the reaction coordinate. The combined Umbrella Sampling (US) and H-REMD approach allows extraction of free energy changes along the coordinate and rapid convergence of the calculations. The simulations are used to interpret the calculated free energy changes and associated conformational changes in the DNA and to follow the onset of local denaturation of the DNA. Comparative simulations were conducted on GC rich


Page 4 of 40

Page 5 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


vs AT rich DNA segments and with different levels of restricting the overall bending of the DNA which results in a stiffer response with respect to twist but also in an earlier onset of melting. Finally, unrestrained MD simulations starting from strongly unwound DNA (with disrupted and partially looped out nucleotides) indicate reversible structure formation of the intact double strand on the time scale of 50-250 ns. The current study gives new insights into how torsional stress affects the thermodynamics, local structure and global topology of DNA at molecular detail.

2 2.1

MATERIALS AND METHODS DNA Structures and Equilibration

All MD simulations were started from two regular B-DNA structures. Both duplexes comprise 15 base pairs and either a central core consisting of G:C base pairs, d(5’-CGCGCGCGCGCGCGC), or A:T base pairs, d(5’-CGCGCATATACGCGC). The DNA structures are termed centGC and centAT, respectively. Simulations were performed with the Amber14 Molecular Dynamics Package 47 . The DNA structures were solvated in explicit solvent (TIP3P water model 48 ) within a rectangular box and a minimum distance of 10 ˚ A between DNA and box boundaries. Potassium ions were added in order to neutralize the system 49 . All simulations were carried out with the pmemd.cuda module in combination with the parmbsc0 force field ( 50 ). The DNA structures are aligned along the z-axis of the box and the systems were first energy minimized in 5000 steps. The systems were heated up to 300 K in three stages including positional restraints on all heavy atoms of the DNA. Each stage entailed a temperature increase of 100 K and was simulated in the NVT ensemble for 100 ps. Positional restraints were subsequently reduced from 25 kcal/(mol·˚ A2 ) to 0.5 kcal/(mol·˚ A2 ) in five consecutive stages, whereby each stage was simulated in the NPT ensemble with a pressure of 1 bar and a temperature of 300 K. The equilibration was finalized by a 2 ns NPT simulation including restraints only on terminal base pairs (cylindrical restraints: see below). 5 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

2.2

Cylindrical restraints

In order to allow free mobility of the dsDNA along the helical z-axis but controlled bending motion, the terminal base pairs of the dsDNAs were cylindrically restrained with respect to the z-axis, i.e. only distances from the z-axis were restrained allowing full freedom in azimuthal directions (illustrated in Figure 1). In terms of cylindrical coordinates this means that just the distance variable (with respect to the z-axis) was restrained and not the angular variable. Hence, in Cartesian coordinates and using a quadratic penalty potential this corresponds to restraining the sum of the squares of x- and y-components. Three cases were investigated: In the case of strong bending suppression/restraining (SBR), the four terminal base pairs at both ends were cylindrically restrained using a force constant of 0.1 kcal/(mol·˚ A2 ). Importantly, the central 7 base pairs are unrestrained. Note, that even in this case most other helical variables of the entire DNA are kept flexible but the application of the cylindrical restraints can be used to control overall bending of the DNA fragment. From an experimental single molecule view this represents the case where supercoiling (i.e. local bending) of the helix is suppressed (typically by applying a small stretching force to keep the DNA straight in single molecule experiments). The second setup, weak bending suppression/restraining (WBR), includes restraints on just the two terminal base pairs at both ends. In this case still limited bending of the helical axis of the dsDNA is possible. Finally, in case of the fully mobile DNA (unrestraint bending: UB) the cylindrical restraints were applied just on the two terminal base pairs at one end of the DNA. This allows full bending mobility of the one end of the DNA with respect to the other end.

2.3

HREUS Simulations

In order to induce untwisting of the B-DNA molecules we used Umbrella Sampling (US) similar to previous studies on DNA twisting near the equilibrium state 27,43 . A torque was applied to the C1’ atoms of the 4th base pair and its symmetric counterpart, the 12th base pair (Figure 1). The torsion angle reaction coordinate was controlled using a quadratic 6 ACS Paragon Plus Environment

Page 6 of 40

Page 7 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


penalty potential,

V = k · ξ − ξiref

2

(1)

The force constant k was 150 kcal/(mol·rad2 ) = 0.0457 kcal/(mol·deg2 ). ξ denotes the reaction coordinate as defined in Figure 1. The reaction coordinate controls the winding of the central 8 base pairs steps and the reference value (ξiref ) was changed in 5◦ intervals from 70◦ to -35◦ during the Umbrella Sampling (US) simulations. This resulted in 22 US intervals. In order to enhance sampling along the reaction coordinate exchanges between neighboring parallel running US intervals were introduced (Hamiltonian Replica Exchange Umbrella Sampling: HREUS). Throughout the HREUS simulations, exchanges between neighboring replicas were attempted every 500 steps and overall each replica attempted at least 20000 exchanges (acceptance rate 0.2-0.5). With a time step of 2 fs this results in a simulation time of at least 20 ns for each replica. Due to convergence issues, specific regions were simulated for longer times in consecutive HREUS simulations of up to 120 ns per US interval. Simulation times for all setups and sequences are listed in the Supplementary Information Table S1. Resulting trajectories were analyzed using VMD 51 and CURVES+ 52 . Free energy profiles were calculated with the weighted histogram analysis method (WHAM) 53 and corresponding error bars were calculated based on thermodynamic integration and block averaging as derived by Zhu and Hummer 54 . Since unwinding induced bending of DNA is anisotropic we evaluated overall bending (curvature) of a segment according to the methodology of Strahs and Schlick 55 : curvature =

θT =

11 X j=4

τj cos(

p θT 2 + θR 2 , with

j X i=4

j X ti ) + ρj sin( ti ) i=4


(2)

(3)


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

θR =

11 X

j j X X −τj sin( ti ) + ρj cos( ti )

j=4

i=4

Page 8 of 40

(4)

i=4

Here, the variables τj and ρj denote tilt and roll of step j and ti is the twist of step i. Note that so defined global roll θR and global tilt θT allow for identifying the bending direction of the segment: θR < 0 means bending towards the minor, θR > 0 towards the major groove, and θT describes bending towards the backbone. In case of disrupted base pairing (local melting of the DNA) it is difficult to define consistent roll, tilt and twist angles and to calculate bending/curvature with the above methodology (from the local tilt, roll and twist at each base pair step). Alternatively, we employed the total bend angle computed with CURVES+ which assigns a smooth global helical axis for the whole DNA molecule and the total bending is largely determined by the relative placement of the two terminal segments of the duplex.

3 3.1

Results and Discussion Free energy change of DNA upon unwinding and comparison to harmonic models

In order to unwind dsDNA during Molecular Dynamics (MD) simulations and to calculate associated changes in free energy a torque was applied to the terminal segments of two dsDNA molecules during Hamiltonian Replica Exchange Umbrella Sampling (HREUS) simulations. The setup is analogous to experimental single molecule DNA torque experiments on DNA 38,39,41 except that experiments require much longer sequences than used in the present simulations. The two 15 bp dsDNA molecules contain either a GC-rich central segment (centGC) or a central AT-rich segment (centAT) (see Methods section). The systems are small enough to allow extensive MD sampling in explicit solvent and, in addition, application of the HREUS advanced sampling technique promotes exchanges of sampled DNA


Page 9 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


conformations between neighboring US windows avoiding to become trapped in metastable states along the reaction coordinate 43 . The simulations allow us to calculate changes in the free energy along the torsion angle reaction coordinate ξ by applying the weighted histogram analysis method (WHAM). Application to a torque on both ends of a dsDNA molecule can affect twisting, bending and stretching of the central segment. In many biological processes DNA can experience a torque but at the same time due to other interactions the bending of DNA is restricted (e.g. due to protein binding or DNA packing). In order to cover such scenarios we performed not only simulations on dsDNA molecules without any restrictions on bending (UB case) but also considered free energy simulations including restraints to restrict DNA bending. This was achieved by applying either weak (WBR case) or strong ”cylindrical” restraints (SBR case) on the terminal segments of the DNA (i.e. keeping the distance of terminal base pairs relative to the helical axis within a preset range without affecting rotational motion around the axis, see Methods for details). In single-molecule experiments this is often mimicked by applying a small stretching force to the ends of a long DNA in order to prevent bending or to limit supercoiling of the DNA. Due to the coupling between twist and local DNA bending it is expected that this may increase the twist stiffness (see analysis below). In the next paragraph we investigate the free energy change vs. reaction coordinate ξ (illustrated in Figure 2) and how it is affected by twisting but also other global variables, the free energy change vs. twist alone (Figure 3) and beyond the harmonic elastic regime is analyzed in a subsequent section. In order to better relate the reaction coordinate to the average DNA twist per base pair step we consider in the following (if not indicated otherwise) mostly 1/8th of ξ as the relevant variable. Since the torque has been applied to a central element of 8 base pair steps this corresponds to a torsion angle per base pair step. For deformations by up to 7◦ per base pair step the calculated free energy change for dsDNA unwinding along the reaction coordinate ξ follows largely a quadratic function (harmonic response, see Figure 2) in case of the HREUS



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 10 of 40

simulations. As will be discussed below, the range is smaller in case of the WBR and SBR simulations. For larger deformations the free energy change flattens and even approaches a near plateau regime (with little or no further free energy change upon further untwisting, see below, and Figure 3). We will first investigate the harmonic regime which we term regime I. In principle, the applied torque can lead to changes in the global helical variables twist, bending and stretching along the central segment. It is of interest to compare the harmonic response of the dsDNA upon induced deformation of up to 7◦ per bp with the flexibility of the same dsDNA extracted from an equilibrium simulation (without torque restaints). Such comparison gives insight into the question how well twist deformability due to induced unwinding is described by the equilibrium properties of the dsDNA observed in the undeformed state. We investigated the dsDNA deformability by comparing calculated PMFs with different variations of an equilibrium harmonic model, relating the deformation free energy to twist, bending and stretching deformations. Such a model has been used by Lankas and coworkers to study local and global DNA deformability 28,30,42 .

F (∆w) =

1 2

∆wT · K∆w

(5)

with ∆w denoting deviations of the selected coordinates from equilibrium. In the harmonic model the deformability of DNA is fully described by the stiffness matrix K, which is determined by the covariance matrix C of the selected coordinates:

K = kB T C −1

(6)

It expresses the equivalence of a description with respect to harmonic force constants or coordinate variances. In the most simple approach we selected only the average twist (accumulated between 4th and 12th base pair) as relevant coordinate. Consequently, the stiffness matrix reduces to the inverse of this coordinate’s variance. 30 The second set of internal coordinates is composed of twist, bending and stretching, where K becomes a 4x4


Page 11 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


matrix since we distinguish between two orthogonal possible bending directions (termed global roll and global tilt of the central segment, see Methods). Furthermore, we applied the harmonic model also directly to the torsion angle reaction coordinate ξ as the only internal coordinate, so that we can assess directly the range of validity and accuracy of this model. Stiffness matrices were calculated from long MD simulations (500 ns) without any torque restraints and ∆w was measured for each replica window. The restraints for limiting DNA bending on the termini were included in the WBR and SBR cases (The elements are given in Supporting Information Table S2). Within the harmonic regime the relative free energies obtained from the HREUS sampling agree very well with the simple harmonic model for the reaction coordinate ξ extracted from equilibrium simulations (compare dotted black and pink curves in Figure 2). Remarkably, using just the simplest harmonic model with the inverse of the twist variance from an unrestrained simulation as the only stiffness force constant agrees also with the HREUS calculations (compare red curve and black dotted curves in Figure 2) in the harmonic regime. Note, that in this case in order to obtain the plot shown in Figure 2 the average twist in each US window was extracted and used to calculate the corresponding harmonic twisting free energy based on the inverse of the twist variance (at equilibrium) as stiffness force constant. Hence, this model includes indirectly also orthogonal degrees of freedom such as bending into account and relates to an effective twist stiffness, Cef f . Typically, the twist elasticity of DNA is given as effective twist persistence length (ratio of DNA length and twist variance). The calculated effective twist force constant in regime I yields a twist persistence length of Cef f =76 nm and 83 nm for centAT and centGC, respectively, (Supporting Information Table S3) in good agreement with experiment 17 . If one considers the second model with a 4x4 stiffness matrix that includes twist and bending deformabilities as diagonal elements and twist-bend (and other) coupling as non-diagonal elements it results also in very close agreement between harmonic model and HREUS-PMF calculations (blue curves in Figure 2) in regime I. However, just using the pure twist stiffness



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

diagonal element form the 4x4 stiffness matrix overestimates the free energy change upon torque deformation especially in the UB and WBR cases (see corresponding green curves in Figure 2)). This stiffness element represents pure twist deformability (without bending or twist-bend coupling) at equilibrium, and translates to an intrinsic twist persistence length, C, of between 102-120 nm (Supporting Information Table S3) in close agreement with experiment 17 and previous simulations 42 . The deviation with respect to the calculated free energy curve using the HREUS approach indicates the significance of the coupling between twisting and bending especially in the UB and to some degree in the WBR cases. Including bending flexibility and twist-bend coupling apparently lowers the free energy curve and results in quantitative agreement between the PMF-HREUS calculations and the harmonic model based just on single equilibrium MD simulations (Figure 2). Similar to the UB case, the result for the harmonic response regime can be almost quantitatively described by parameters extracted from equilibrium simulations on twist stiffness, bending stiffness and twist-bend coupling also in the WBR and SBR cases. Since bending is increasingly restricted in the WBR and SBR cases, the change in free energy along ξ can be described by just considering only the twisting stiffness diagonal element. In the SBR case the twist diagonal stiffness element alone gives already close agreement with the PMFHREUS simulations (compare green line and other lines in Figure 2). This result confirms the expectation that under conditions of strong bending suppression, the torque induced free energy change is entirely stored in twist deformation. Thus, we conclude that deformation energy caused by torsional stress is nearly fully absorbed in untwisting and local bending of the DNA. This result also emphasizes that DNA bending plays an essential role in relaxing torsional stress on DNA (for longer DNA exactly this torsion-induced bending is the cause of supercoiling). Another interesting result of this comparison concerns the validity range of a harmonic deformability model: Although the fluctuations in mean twist during equilibrium simulations in the UB case are considerably smaller than 7◦ (around 3 − 3.5◦ ) the extracted parameters


Page 12 of 40

Page 13 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


are apparently sufficient to accurately describe larger deformations (see Figure 2). A third remarkable result concerns the breakdown of the harmonic response which depends on the restriction of bending deformations: Similar to the UB case we observe a harmonic response regime of the PMF along ξ for both the WBR and SBR case for small deformations but it is restricted to a smaller deformation range and with a stiffer response (Figure 2). It indicates that as explained above restricted bending also causes an effective increase in twist stiffness. Note, that this is also the case if one restrains bending of DNA not to something close to a straight DNA but also to other defined bending magnitudes, e.g. in case of forming circularily closed DNA or upon DNA packing or in single molecule experiments. Indeed, for such cases an increase of twist stiffness has been experimentally observed 8,17,18 .

3.2

Free energy change upon DNA unwinding beyond the harmonic regime

Comparison of the calculated free energy from the HREUS simulations as a function of average twist of the central dsDNA segment indicates a slightly larger stiffness of the GC rich centGC segment compared to the centAT case. Indeed, a small sequence dependence of the torsional stiffness of DNA has also been observed experimentally with GC rich sequences being slightly stiffer than AT rich regions 12 . This may have implications for the specific recognition by minor-groove binding proteins (e.g. TATA box binding) since local untwisting leads to opening of the minor groove. In the UB case the breakdown of the harmonic response sets in around an average twist of 26◦ (Figure 3). Note, that the calculated equilibrium twist (for both dsDNAs) indicates undertwisting of DNA relative to experimental B-DNA (Figure 3). However, this is a known deficiency of the parmbsc0 force field 23,50 . The recent updates of the force field (parmbsc1 25 or the ol15 forcefields 56 ) largely correct this deficiency without appreciably changing the deformability of B-DNA 23 . The calculated free energy per base pair step at the limit of the harmonic response 13 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

amounts to approximately 2kB T (=1.2 kcal/mol) in the centGC and is slightly smaller for the centAT case (by 0.2 kcal/mol). Note, that this corresponds to a collective onset of melting of a segment of base pairs. The disruption of a single base pair in an otherwise intact DNA may require larger free energy barriers. As already mentioned upon further induced untwisting below an average of 25◦ − 26◦ the harmonic description breaks down with a significant reduction of the free energy increase (Figure 3). Associated structural changes correspond to a disruption of base pairs, hence, onset of local DNA melting (analyzed in the next paragraph). In case of the WBR-HREUS simulations, the breakdown of the harmonic free energy response is predicted to occur at a significantly higher average twist of 27◦ and in the SBR case the DNA melting sets in at even higher values of the reaction coordinate corresponding to an average twist of 28◦ . The free energy per base pair step required to cause a breakdown of the harmonic response is generally slightly lower for centAT compared to centGC but of similar magnitude for UB, WBR and SBR setups (Figure 3). Note, that the calculated PMFs do not indicate any local minima in the melting regime but a rather smooth free energy landscape. It should be emphasized that this result is still compatible with several possible locally stable states that densely fill the space along the reaction coordinate and are of similar stability. Although the absence of distinct visible free energy minima can still be due to insufficient sampling of relevant states in this regime, it is important to note and discussed in the last paragraph that the sampling of states appears to be largely reversible. The transition of a locally strongly disrupted DNA with a ”melted” central segment was revealed by starting unrestrained MD simulations from severely untwisted DNA duplex structures which yield eventually intact B-DNA in less than ∼ 50 ns for the centAT case and within 250 ns in case of the centGC case (see last paragraph of the Results and Discussion section).


Page 14 of 40

Page 15 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.3

Torque induced helical conformational changes

Application of the torque on the DNA results in unwinding of the DNA. As indicated in Figure 4 and expected for a harmonic response we find a linear relation between average twist and external torsion reaction coordinate within the harmonic regime which is indicated as regime I in the plots of Figure 4. For all sequences and bending restraint setups we find that the average twist (over the central 8 bp steps) is simultaneously decreased with ξ and one can distinguish two regimes in addition to the harmonic regime I (Figure 4). If the external stress causes only untwisting a linear relation with a slope close to 1 would be expected in all regimes. As discussed in the previous section, in addition to untwisting, the external torque also results in bending of the DNA (illustrated in Figure 4A,B). This bending is locally due to roll-twist and twist-tilt coupling. For the harmonic regime I, we already compared the twistbend coupling between unrestrained equilibrium simulations and HREUS simulations and found excellent agreement for these two approaches. Exact correlations between torsional angle ξ and average twist were determined by linear interpolation in each of the three regimes. The fitted slopes indicate a correlation of about 88% for UB case and 92% for the WBR setup (Table 1). Similar to our conclusions on the energetics of DNA untwisting (see previous paragraph) these results reveal that the external torsional stress is not fully accommodated by untwisting but also by bending deformation. In line with this, the correlation in regime I amounts to almost exactly 100 % for the SBR case (Table 1), which largely suppresses bending. The torque induced bending has a specific direction and we also compared the relation between bending direction and induced unwinding. The direction of bending of the central segment is given by the above defined global roll θR and global tilt θT of the segment. Remarkably, the correlation between average twist and global roll as well as global tilt extracted from the unrestraint equilibrium simulations is in close agreement (Supporting Information Figure S1). It means that in regime I the torque induced bending follows closely the bending/twisting direction that is sampled in unrestrained equilibrium simulations. In regime I 15 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

the HREUS simulations sample conformations that are significantly deformed compared to B-DNA but show still largely an intact base pairing and stacking. The various correlations between twist changes and other local base pair helical variables in this regime have been extensively investigated in previous equilibrium MD simulations 19,23,28,30 and in our previous studies that analyzed the effect of a torque on dsDNA but limited to small twist deformations 27,43 . Thus, we limit the discussion to helical variables that show dramatic changes upon transition between the regimes along the reaction coordinate. While the total bend angle increases steadily upon unwinding in regime I, it begins to drop in regime II (Figure 4) resulting in a loss of the strict bend-twist coupling observed in regime I and discussed above. Furthermore, the total bend angle is no longer composed of continuous bending (due to uniform small contributions of each base pair step of the central DNA segments) but is determined by local strong kinking at bp steps due to local structural failure (this is illustrated in Figure 5). Due to this rather abrupt loss of bend-twist coupling (as a mechanism to absorb torsional stress) in regime II the external torsional stress must be compensated by increased untwisting in the central segment. Exactly this is observed in regime II resulting in an excessive coupling between ξ and average twist (Table 1). In regime III, however, again a linear relation between external torque and twist is observed with a slope closer to 1 (Table 1, there is still rather irregular bending of the DNA, however, without a strict dependence on twist, see Figure 6). The response of the DNA to the external torque changes rather abruptly upon reaching the limits of the harmonic regime which corresponds to a qualitative change in the calculated PMF for twist deformation from a quadratic to a nearly flat linear regime (see last paragraph). In this transition phase (regime II) DNA base pairings start to disrupt in the central region. Especially in the centAT case various flipping events of bases and unstacking events are observed. For the centGC case, in contrast, several of the base pairs remain intact but the nucleic acid backbone structure starts to adopt altered conformations and the bases form altered stacking patterns (illustrated in Figure 5 and 7). Whereas mainly stacking dis-


Page 16 of 40

Page 17 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


turbances, left handed segments and base flippings were sampled for the centGC sequence, the centAT sequence was locally melted with increasing frequency for lower twist values. This different behavior is reflected also in higher opening angles for the AT sequence (Figure 6). Hence, the simulations indicate that the accessibility of the melted DNA segment may also depend on the base pair composition in the relevant segment. In regime I, the helical rise of the DNA decreased continuously upon unwinding down to values characteristic for A-DNA (Figure 6). This decrease in helical rise upon unwinding is a reflection of negative twist-stretch coupling of B-DNA 38 . However, in regime II the helical rise increases abruptly, relaxing the tension by rapid extension (increase in helical rise) of the DNA in the helical axis direction (due to the structural failure of the central segment and local melting, illustrated in Figure 5) and loss of the compact bend A-type like central structure and a reduced diameter of the double helix . We also looked at the distribution of the mean twist (average over the 8 central base pair steps) in individual Umbrella Sampling (US) windows. Interestingly, single narrow approximately Gaussian distributions are observed in regime I (Figure 8) but upon reaching regime II the histograms represent rather bimodal distributions as a consequence of a greater variety of sampled partially melted subsegments mixed with still regular helical substructures in the central DNA segments (Figure 8). In regime III the twist distributions appear more narrowly distributed around single Gaussian distributions again. It should be emphasized that a Gaussian distribution of the mean twist (over the 8 central base pair steps) does not mean that the local twist distribution at every base per step also follows a Gaussian distribution. Indeed, it has been observed that the local twist at every step is sequence dependent and can show even a bimodal distribution at certain steps, e.g. 19–25 . We have found in a previous study 27 that in the elastic regime the response to an external torque is also significantly sequence dependent (e.g. pyrimidine-purine steps react much more strongly to an external torque strain than purine-purine or purine-pyrimidine steps). In future studies we will systematically study the change in substate distribution and local twist distribution



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

upon external torques.

3.4

Torque induced backbone conformational changes

As discussed in the previous paragraph the unwinding of the dsDNA in the regime I leads to local bending (towards the major groove) and opening of the minor groove as well as reduction in helical rise towards a local A-DNA-like conformation. This is also reflected in an increase of the proportion of the O1’-endo (A-form type) sugar pucker (Figure 9) although the typical C3’-endo (A-form) sugar pucker was only marginally sampled. For both the centAT and the centGC cases the proportion of O1’-endo increased with decreasing twist at the cost of the C2’-endo pucker characteristic for B-DNA. At the breakdown of the harmonic response regime I this tendency abruptly changes and upon further unwinding (and local breakdown of the double helical structure) the proportion of C2’-endo increases whereas the contribution of O1’-endo decreases towards a ratio similar to B-DNA. The proportion of the C1’-exo sugar pucker remains for both cases approximately constant (Figure 9). In addition, the induced twist deformation also affects the distribution of dihedral conformational substates in DNA (Figure 9). Transitions from canonical α/γ (g-/g+) to non-canonical states as well as significant reductions of the BI/BII population point out that untwisting has a stronger influence on the backbone conformation for the centGC sequence than for the centAT. Here, the BI substate corresponds to the regular canonical combination of the coupled /ζ backbone dihedral angle pair and the BII is a frequent non-canonical substate (Figure 9). The result is in accordance with equilibrium MD simulations of B-DNA that indicate low twist of DNA coupled to decreased proportions of BII states in DNA 19 . In the regime II and III due to the disruption of base pairing and partial unstacking a variety of conformational substates become accessible which results in decreased population of the canonical BI and also increased substates involving non-canonical α/γ states.


Page 18 of 40

Page 19 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


3.5

Conformational transition of locally melted DNA towards intact B-DNA

An interesting question is if local DNA melting results in conformational states that can relax back rapidly to regular B-DNA upon removal of the external torque or are kinetically trapped in long-lived meta-stable subminima. Starting from several distorted DNA conformations obtained from the regime III (at the final stage of the HREUS simulations) unrestrained MD simulations were performed. As illustrated in Figure 10 the initial centrally disrupted duplex structure underwent transitions back to fully intact B-DNA after ∼ 50 ns in the centAT case with several flipped out melted nucleotides in the central segment in the start structure. Surprisingly, however, the average twist approaches its equilibrium value already after only ∼ 10 ns, thus indicating that the global helical structure relaxes on a shorter timescale than the total local relaxation process. Inspection of the first 10 ns of simulation indicate that the twist relaxation is triggered by a characteristic preceding alteration in the helical structure (illustrated as a sequence of snapshots in Figure 10). The bending increases transiently up to ∼ 60◦ and the helical rise of the central segment drops to ∼ 2.5 ˚ A. Remarkably, this configuration resembles globally the transition state from phase I to II that was observed in the HREUS unwinding simulations. It is characterized by a minimal helical rise and an excessive bending of the central segment. After the twist and subsequently other helical parameters have relaxed the whole DNA approaches an equilibrium B-DNA configuration and subsequently fluctuates closely around an equilibrium ensemble. For the centAT case, this includes opening of disrupted base pairs (Figure 10) but also the relaxation of completely flipped nucleotides (in a looped out solvent exposed conformation) back into a fully based-paired B-DNA structure. In contrast to the centAT case, the centGC structure adopted mostly irregular but frequently still base paired conformations of the central segment in the strongly unwound regime III (Figure 7). Starting unrestrained simulations from such conformations resulted in spontaneous abrupt transitions to regular B-DNA on longer times scales of > 200 ns (Figure 19 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

10).

4

Conclusions

Local unwinding of DNA is essential in many biological processes such as supercoiling, packing or damage repair 1,2,5 . It is even proposed to be the underlying mechanistic basis for the melting and disruption of base pairing to trigger subsequent enzymatic processes during transcription and replication 2,32 . The twist deformability near an equilibrium B-DNA state and its coupling to other helical parameters can be described quite well with a harmonic stiffness model 28,30 or based on empirical parameters in the framework of treating DNA as an elastic deformable material 12,16 . It should be emphasized that here the mean twist over a stretch of base pair steps is considered and not a local twist at each base pair step that can follow non-Gaussian or even bimodal distributions. For the harmonic regime of mean twist changes the present MD study provides realistic predictions on the effective and intrinsic twist stiffness and the magnitude and direction of coupled bending and stretching. Beyond previous simulation studies, it indicates the limits of the validity of a harmonic model both with respect to the twist coordinate and in addition with respect to the required free energy to onset local melting. It was found that the unwinding process can be best described by distinguishing three regimes with regime I representing the harmonic response on the external torque. It is characterized by a quadratic free energy change for twist deformations and an associated strict coupling of twist to changes in stretching/bending. At regime II, the onset of DNA melting, the coupling of twist to bending and stretching breaks down. It results in an extension of the DNA and stochastic bending of the DNA due to the occurrence of local kinks as a result of deformed base pair steps and flipped out nucleotides. In this regime, the free energy change associated with further unwinding starts to drop to reach a relatively flat form in regime III with much reduced resistance of the system to further unwinding. Interestingly, torsional stress disrupted the base pairing between adenine and thymine bases


Page 20 of 40

Page 21 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


and hence also caused local strand separation, whereas for G-C base pairs the backbone structure and stacking pattern was predominantly deformed in regime III. Remarkably, although the structural alterations in the central segments upon unwinding were significant, it was possible to recover an intact B-DNA structure starting from such deformed structures during unrestrained MD simulations in the range of 50-250 ns. This result indicates that a local melting event in an otherwise dsDNA may have only a relatively short lifetime (few tenths or hundreds of ns) after removal of the external torque stress. This agrees with NMR experiments that predict nanosecond life times of flipped out bases in DNA. Although the present study gives also some insight into the sequence dependence of the effect of external unwinding of DNA, future studies on other sequences are important to characterize the response of mixed GC and AT containing segments. In a previous study 27 the local response at the the base pair step level to an external torque was already investigated. It will be of interest to also consider the distribution of local conformational substates of individual base pair steps under an external torque and how it differs from the distribution without a torque. In addition, it is also of interest to investigate the response of mismatched or damaged DNA and of the surrounding ion atmosphere to external torque. Furthermore, the overwinding of DNA that we already investigated in a previous study within the elastic regime will be tackled in future work although it is biologically of lesser importance than DNA unwinding. Finally, it is also interesting to compare the unwinding response of dsRNA with dsDNA.

5

Supporting Information Available

Tabulated Umbrella Sampling simulation details, stiffness parameters, calculated twist presistence lengths and global bending direction coupled to twist in unrestrained and Umbrella Sampling simulations



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

6

ACKNOWLEDGEMENTS

The authors thank Jan Lipfert and Filip Lankas for insightful discussions. The project was funded by the Deutsche Forschungsgemeinschaft (DFG) grant number SFB749 project B06. Computational resources were provided by the Munich Leibniz Supercomputing Center (LRZ) grant number pr48po.

References (1) Lavelle, C. DNA Torsional Stress Propagates Through Chromatin Fiber and Participates in Transcriptional Regulation. Nat. Struct. Mol. Biol. 2008, 15, 123–125. (2) Kim, T.-K.; Ebright, R. H.; Reinberg, D. Mechanism of ATP-Dependent Promoter Melting by Transcription Factor IIH. Science 2000, 288, 1418–1421. (3) Velmurugu, Y.; Chen, X.; Slogoff Sevilla, P.; Min, J.-H.; Ansari, A. Twist-Open Mechanism of DNA Damage Recognition by the Rad4/XPC Nucleotide Excision Repair Complex. Proc. Natl. Acad. Sci. U.S.A. 2016, 113, 2296–2305. (4) Strick, T. R.; Croquette, V.; Bensimon, D. Single-Molecule Analysis of DNA Uncoiling by a Type II Topoisomerase. Nature 2000, 404, 901–904. (5) Collier, C.; Machón, C.; Briggs, G. S.; Smits, W. K.; Soultanas, P. Untwisting of the DNA Helix Stimulates the Endonuclease Activity of Bacillus Subtilis Nth at AP Sites. Nucleic Acids Res. 2011, 40, 739–750. (6) Collins, I.; Weber, A.; Levens, D. Transcriptional Consequences of Topoisomerase Inhibition. Mol. Cell Biol. 2001, 21, 8437–8451. (7) Rohs, R.; Jin, X.; West, S. M.; Joshi, R.; Honig, B.; Mann, R. S. Origins of Specificity in Protein-DNA Recognition. Annu. Rev. Biochem. 2010, 79, 233–269.


Page 22 of 40

Page 23 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(8) Heath, P. J.; Clendenning, J. B.; Fujimoto, B. S.; Schurr, M. J. Effect of Bending Strain on the Torsion Elastic Constant of DNA. J. Mol. Biol. 1996, 260, 718–730. (9) Crothers, D. M.; Drak, J.; Kahn, J. D.; Levene, S. D. [1] DNA Bending, Flexibility, and Helical Repeat by Cyclization Kinetics. Methods Enzymol. 1992, 212, 3–29. (10) Taylor, W. H.; Hagerman, P. J. Application of the Method of Phage T4 DNA LigaseCatalyzed Ring-Closure to the Study of DNA Structure: II. NaCl-Dependence of DNA Flexibility and Helical Repeat. J. Mol. Biol. 1990, 212, 363–376. (11) Bryant, Z.; Stone, M. D.; Gore, J.; Smith, S. B.; Cozzarelli, N. R.; Bustamante, C. Structural Transitions and Elasticity from Torque Measurements on DNA. Nature 2003, 424, 338–341. (12) Fujimoto, B. S.; Schurr, J. M. Dependence of the Torsional Rigidity of DNA on Base Composition. Nature 1990, 344, 175. (13) Barkley, M. D.; Zimm, B. H. Theory of Twisting and Bending of Chain Macromolecules; Analysis of the Fluorescence Depolarization of DNA. J. Chem. Phys. 1979, 70, 2991– 3007. (14) Millar, D.; Robbins, R.; Zewail, A. Torsion and Bending of Nucleic Acids Studied by Subnanosecond Time-Resolved Fluorescence Depolarization of Intercalated Dyes. J. Chem. Phys. 1982, 76, 2080–2094. (15) Strick, T.; Allemand, J.; Bensimon, D.; Bensimon, A.; Croquette, V. The Elasticity of a Single Supercoiled DNA Molecule. Science 1996, 271, 1835. (16) Moroz, J. D.; Nelson, P. Torsional Directed Walks, Entropic Elasticity, and DNA Twist Stiffness. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 14418–14422. (17) Nomidis, S. K.; Kriegel, F.; Vanderlinden, W.; Lipfert, J.; Carlon, E. Twist-Bend



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Coupling and the Torsional Response of Double-Stranded DNA. Phys. Rev. Lett. 2017, 118, 217801. (18) Schurr, J. M. Possible Origin of the Increased Torsion Elastic Constant of Small Circular DNAs: Bending-Induced Axial Tension. J. Phys. Chem. B 2017, 121, 5709–5717. (19) Pasi, M.; Maddocks, J. H.; Beveridge, D.; Bishop, T. C.; Case, D. A.; Cheatham, T., III; Dans, P. D.; Jayaram, B.; Lankas, F.; Laughton, C.; et al., µABC: A Systematic Microsecond Molecular Dynamics Study of Tetranucleotide Sequence Effects in B-DNA. Nucleic Acids Res. 2014, 42, 12272–12283. (20) Beveridge, D. L.; Cheatham, T. E.; Mezei, M. The ABCs of Molecular Dynamics Simulations on B-DNA, circa 2012. J. Biosci. 2012, 37, 379–397. (21) Dans, P. D.; Pérez, A.; Faustino, I.; Lavery, R.; Orozco, M. Exploring Polymorphisms in B-DNA Helical Conformations. Nucleic Acids Res. 2012, 40, 10668–10678. (22) Dans, P. D.; Faustino, I.; Battistini, F.; Zakrzewska, K.; Lavery, R.; Orozco, M. Unraveling the Sequence-Dependent Polymorphic Behavior of d (CpG) Steps in B-DNA. Nucleic Acids Res. 2014, 42, 11304–11320. (23) Dans, P. D.; Danilane, L.; Ivani, I.; Drsata, T.; Lankas, F.; Walther, J.; Pujagut, R. I.; Battistini, F.; Gelp´ı, J. L.; Lavery, R. Long-Timescale Dynamics of the Drew–Dickerson Dodecamer. Nucleic Acids Res. 2016, 44, 4052–4066. (24) Balaceanu, A.; Pasi, M.; Dans, P. D.; Hospital, A.; Lavery, R.; Orozco, M. The Role of Unconventional Hydrogen Bonds in Determining BII Propensities in B-DNA. J. Phys. Chem. Lett. 2016, 8, 21–28. (25) Ivani, I.; Dans, P. D.; Noy, A.; Pérez, A.; Faustino, I.; Hospital, A.; Walther, J.; Andrio, P.; Go˜ ni, R.; Balaceanu, A. Parmbsc1: A Refined Force Field for DNA Simulations. Nat. Methods 2016, 13, 55–58. 24 ACS Paragon Plus Environment

Page 24 of 40

Page 25 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


ˇ (26) Lankaˇs, F.; Sponer, J.; Langowski, J.; Cheatham, T. E. DNA Basepair Step Deformability Inferred from Molecular Dynamics Simulations. Biophys. J. 2003, 85, 2872–2883. (27) Kannan, S.; Kohlhoff, K.; Zacharias, M. B-DNA Under Stress: Over- and Untwisting of DNA during Molecular Dynamics Simulations. Biophys. J. 2006, 91, 2956 – 2965. (28) Drsata, T.; Pérez, A.; Orozco, M.; Morozov, A. V.; Sponer, J.; Lankas, F. Structure, Stiffness and Substates of the Dickerson-Drew Dodecamer. J. Chem. Theory Comput. 2012, 9, 707–721. (29) Lankaˇs, F.; Gonzalez, O.; Heffler, L.; Stoll, G.; Moakher, M.; Maddocks, J. H. On the Parameterization of Rigid Base and Basepair Models of DNA from Molecular Dynamics Simulations. Phys. Chem. Chem. Phys. 2009, 11, 10565–10588. (30) Drsata, T.; Spavkova, N.; Jurecka, P.; Zgarbova, M.; Sponer, J.; Lankas, F. Mechanical Properties of Symmetric and Asymmetric DNA A-tracts: Implications for Looping and Nucleosome Posi tioning. Nucleic Acids Res. 2014, 42, 7383–7394. (31) Bignon, E.; Drsata, T.; Morell, C.; Lankas, F.; Dumont, E. Interstrand Cross-Linking Implies Contrasting Structural Consequences for DNA: Insights from Molecular Dynamics. Nucleic Acids Res. 2017, 45, 2188–2195. (32) Jeon, J.-H.; Adamcik, J.; Dietler, G.; Metzler, R. Supercoiling Induces Denaturation Bubbles in Circular DNA. Phys. Rev. Lett. 2010, 105, 208101. (33) Irobalieva, R. N.; Fogg, J. M.; Catanese Jr, D. J.; Sutthibutpong, T.; Chen, M.; Barker, A. K.; Ludtke, S. J.; Harris, S. A.; Schmid, M. F.; Chiu, W.; Zechiedrich, L. Structural Diversity of Supercoiled DNA. Nat. Commun. 2015, 6 . (34) La Rosa, G.; Zacharias, M. Global Deformation Facilitates Flipping of Damaged 8-OxoGuanine and Guanine in DNA. Nucleic Acids Res. 2016, 44, 9591–9599.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(35) Knips, A.; Zacharias, M. Both DNA Global Deformation and Repair Enzyme Contacts Mediate Flipping of Thymine Dimer Damage. Sci. Rep. 2016, 7, 41324. (36) Strick, T. R.; Allemand, J. F.; Bensimon, D.; Croquette, V. Behavior of Supercoiled DNA. Biophys. J. 1998, 74, 2016–2028. (37) Strick, T.; Allemand, J.-F.; Bensimon, D.; Croquette, V. Stress-Induced Structural Transitions in DNA and Proteins. Annu. Rev. Biophys. 2000, 29, 523–543. (38) Lipfert, J.; Skinner, G. M.; Keegstra, J. M.; Hensgens, T.; Jager, T.; Dulin, D.; Köber, M.; Yu, Z.; Donkers, S. P.; Chou, F.-C. Double-Stranded RNA under Force and Torque: Similarities to and Striking Differences from Double-Stranded DNA. Proc. Natl. Acad. Sci. U.S.A. 2014, 111, 15408–15413. (39) Kriegel, F.; Ermann, N.; Lipfert, J. Probing the Mechanical Properties, Conformational Changes, and Interactions of Nucleic Acids with Magnetic Tweezers. J. Struct. Biol. 2017, 197, 26 – 36. (40) Gore, J.; Bryant, Z.; Stone, M. D.; Nöllmann, M.; Cozzarelli, N. R.; Bustamante, C. Mechanochemical Analysis of DNA Gyrase Using Rotor Bead Tracking. Nature 2006, 439, 100–104. (41) Lipfert, J.; Kerssemakers, J. W.; Jager, T.; Dekker, N. H. Magnetic Torque Tweezers: Measuring Torsional Stiffness in DNA and RecA-DNA Filaments. Nat. Methods 2010, 7, 977–980. (42) Lankas, F.; Sponer, J.; Hobza, P.; Langowski, J. Sequence-Dependent Elastic Properties of DNA. J. Mol. Biol. 2000, 299, 695–709. (43) Liebl, K.; Drsata, T.; Lankas, F.; Lipfert, J.; Zacharias, M. Explaining the Striking Difference in Twist-Stretch Coupling between DNA and RNA: A Comparative Molecular Dynamics Analysis. Nucleic Acids Res. 2015, 43, 10143–10156. 26 ACS Paragon Plus Environment

Page 26 of 40

Page 27 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


(44) Marko, J. F.; Siggia, E. D. Fluctuations and Supercoiling of DNA. Science 1994, 506– 506. (45) Wereszczynski, J.; Andricioaei, I. On Structural Transitions, Thermodynamic Equilibrium, and the Phase Diagram of DNA and RNA Duplexes under Torque and Tension. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 16200–16205. (46) Randall, G. L.; Zechiedrich, L.; Pettitt, B. M. In the Absence of Writhe, DNA Relieves Torsional Stress with Localized, Sequence-Dependent Structural Failure to Preserve B-form. Nucleic Acids Res. 2009, 37, 5568–5577. (47) Case, D.; Babin, V.; Berryman, J.; Betz, R.; Cai, Q.; Cerutti, D.; Cheatham III, T.; Darden, T.; Duke, R.; Gohlke, H. The FF14SB Force Field. Amber 2014, 14, 29–31. (48) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926–935. (49) Joung, I. S.; Cheatham III, T. E. Molecular Dynamics Simulations of the Dynamic and Energetic Properties of Alkali and Halide Ions Using Water-Model-Specific Ion Parameters. J. Phys. Chem. B 2009, 113, 13279–13290. (50) Pérez, A.; Marchán, I.; Svozil, D.; Sponer, J.; Cheatham, T. E.; Laughton, C. A.; Orozco, M. Refinement of the AMBER Force Field for Nucleic Acids: Improving the Description of α/γ Conformers. Biophys. J. 2007, 92, 3817–3829. (51) Humphrey, W.; Dalke, A.; Schulten, K. VMD: Visual Molecular Dynamics. J. Molec. Graphics 1996, 14, 33–38. (52) Lavery, R.; Moakher, M.; Maddocks, J.; Petkeviciute, D.; Zakrzewska, K. Conformational Analysis of Nucleic Acids Revisited: Curves+. Nucleic Acids Res. 2009, 37, 5917–5929. 27 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

(53) Kumar, S.; Rosenberg, J. M.; Bouzida, D.; Swendsen, R. H.; Kollman, P. A. The Weighted Histogram Analysis Method for Free-Energy Calculations on Biomolecules. I. The Method. J. Comput. Chem. 1992, 13, 1011–1021. (54) Zhu, F.; Hummer, G. Convergence and Error Estimation in Free Energy Calculations Using the Weighted Histogram Analysis Method. J. Comput. Chem. 2012, 33, 453–465. (55) Strahs, D.; Schlick, T. A-Tract Bending: Insights into Experimental Structures by Computational Models. J. Mol. Biol. 2000, 301, 643–663. (56) Zgarbova, M.; Sponer, J.; Otyepka, M.; Cheatham III, T. E.; Galindo-Murillo, R.; Jurecka, P. Refinement of the Sugar–Phosphate Backbone Torsion Beta for Amber Force Fields Improves the Description of Z-and B-DNA. J. Chem. Theory Comput. 2015, 11, 5723–5736.


Page 28 of 40

Page 29 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Table 1: Slope for the average twist of the central 8 base pair steps vs. torsion reaction coordinate ξ AT GC regime slope UB slope WBR slope SBR slope UB slope WBR slope SBR I 0.886 0.920 0.998 0.861 0.905 1.01 II 1.55 2.86 1.33 1.93 2.04 1.37 III 1.27 1.21 1.03 1.24 1.07 1.01 The indicated slopes corrrespond to the slopes of the average twist vs. ξ as illustrated in Figure 4



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 1: Illustration of the Umbrella Sampling free energy simulation setup. The dsDNA is indicated as stick model (unrestrained part: atom-color coded). In the weak or strong bending restraining (WBR or SBR) simulations the terminal 2 or 4 base pairs, respectively, were cylindrically restrained to the helical (z-)axis (illustrated as cylinders surrounding the terminal base pairs). These restraints allow complete rotational freedom but restrict bending deformability. The torque to induce unwinding (black arrow) acts on four centers (blue spheres) labelled A-D.


Page 30 of 40

Page 31 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 2: Comparison of the calculated free energy vs. reaction coordinate ξ (black dotted curves) obtained from HREUS simulations with different harmonic models based on the stiffness parameters from equilibrium simulations without torque restraints (in the following termed equilibrium simulations). The different systems and restraints are indicated in the insets of each panel. In case of the dotted pink curve a stiffness parameter for ξ was directly calculated from an equilibrium simulation (UB case, upper panels), including weak bending restraints (WBR, middle panels) or strong bending restraints (SBR, lower panels). For the red curve a stiffness parameter for twist was calculated directly from the twist fluctuations during equilibrium simulations. To translate this into a free energy vs. ξ, the average twist in each HREUS window was extracted and used to calculate a free energy based on a quadratic model using the twist stiffness as force constant. For the blue curve the same procedure was used but employing a stiffness matrix including a twist, stretch and bending stiffness as diagonal components and coupling terms each extracted from fluctuations extracted from equilibrium simulations. The green curve represents a harmonic model based on just the diagonal twist stiffness parameter of the 4x4 stiffness matrix. 31 ACS Paragon Plus Environment


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3: Calculated free energy vs. average twist of the central segment for all setups and both sequences extracted from the HREUS simulations. The calculation of uncertancies in the free energy profiles is described in the methods section.


Page 32 of 40

Page 33 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4: Relation between torsional angle ξ/8 and average twist over the central 8 bp steps. Only data of the last 20 ns of each replica were used, ξ was measured every ps and snapshots were taken every 10 ps. For each replica, an average was calculated for ξ and for the average twist over the central 8 bp steps. Average twist values were determined with CURVES+. Thereby, only the 4th, 8th and 12th base pair were taken into account, the twist values computed for the steps 4 → 8 and 8 → 12 were added and divided by 8. This procedure was found to perform robustly in regime II and III where the significant distortion (opening and flipping of nucleotides) of the DNA often prevents a reliable extraction of the twist variable for each individual base pair step. Linear interpolation was performed for each regime.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 5: Illustration of global conformational changes during HREUS unwinding simulations of the centAT (upper panels) and centGC cases. (A) Snapshots taken from simulation without external torques (helical axis calculated by CURVES+ is indicated). (B) Snapshots from the unwinding regime I prior to unwinding induced melting. (C) Snapshots taken from regime II/III indicating less bending of the DNA but structural distortions in the central segments

Figure 6: Changes in global bending (left panels), helical rise (middle panels) and base pair opening angle (right panels) induced during HREUS simulations along ξ for the UB case. The plot was generated similarly as described in the caption of Figure 4.


Page 34 of 40

Page 35 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 7: Structural snapshots illustrated as stick models (four base pairs) taken from Umbrella sampling windows in regime III, A:T base pairs are in blue (panel A, B), G:C base pairs are shown in green (panel C, D)



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 8: Probability distribution of the mean twist (average over the 8 central base pair steps) in individual Umbrella Sampling windows. In regime I narrow Gaussian type distributions were sampled in each US window (A:centAT, D:centGC) whereas broader and more complicated distributions were obtained in the regime II (B:centAT, E:centGC) followed again by more narrow distributions in regime III (C:centAT, F:centGC). The broader distribution in regime II (panel B, E) is due to a greater variety of states of the central segment sampled in the regime II and III compared to regime I which includes also states with different twist and bending variables.


Page 36 of 40

Page 37 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 9: Relation between reaction coordinate ξ and sampled nucleic acid backbone conformations characterized by the ratio of backbone dihedrals /ζ (BII/BI ratio) and α/γ in the central segment (A: centAT (UB) case, C: centGC (UB) case). The fraction of sampled sugar pucker states in the central segment is indicated in (B: centAT case) and (D: centGC case).



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 10: Relaxation of unwound DNA towards regular B-DNA in unrestraint MD simulations. (a) Time evolution of global parameters twist, helical rise and bending during 100 ns unrestraint MD simulation starting from a strongly deformed unwind centAT structure. The global helical parameters are given relative to the values in B-DNA. Snapshots of the starting structure and characteristic intermediate structures are illustrated in the inset of panel (a). The corresponding simulation times associated with each snapshot are labelled in the plot (snapshots A-D). (b) Time evolution of twist, helical rise and bending in case of starting an unrestraint MD simulation from a strongly unwind centGC conformation.


Page 38 of 40

Page 39 of 40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 11: Toc-Figure

Figure 12: Korbinian Liebl (1991) obtained his M.Sc. in Physics at the Technical University of Munich in 2016. Since then, he is engaged as PhD student at the chair for Theoretical Biophysics. His major research interests concern the mechanics of DNA and RNA as well as the repair of DNA damages.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 13: Martin Zacharias is Professor of Theoretical Biophysics and Head of the Biomolecular Dynamics group at Technical University Munich. He obtained his PhD from Max Planck Institute of Molecular Genetics and the Free University of Berlin in 1991. After postdoctoral stays at University of Houston (with Prof. J. Andrew McCammon) and University of Colorado (with Prof. Paul Hagerman) he headed a research group at the Institute of Molecular Biotechnology in Jena, Germany (1999-2003). He became affiliated with Jacobs University Bremen as Professor of Computational Biophysics in 2003 and in 2009 joined Technical University Munich to take his current position. His research interests focus on the structure and dynamics of proteins and nucleic acids using molecular simulation and docking techniques. It also includes the study of structure formation processes and the understanding of nucleic acid flexibility and how it is coupled to recognition and binding of proteins.


Page 40 of 40

Unwinding Induced Melting of Double-Stranded DNA Studied by Free

Recommend Documents