Sequence dependence of transient Hoogsteen base pairing in DNA.

Hoogsteen (HG) base pairing is characterized by a 180° rotation of the purine base with respect to the Watson-Crick-Franklin (WCF) motif. Recently, it has been found that both conformations coexist in a dynamical equilibrium and that several biological functions require HG pairs. This relevance has motivated experimental and computational investigations of the base-pairing transition. However, a systematic simulation of sequence variations has remained out of reach. Here, we employ advanced path-based methods to perform unprecedented free-energy calculations. Our methodology enables us to study the different mechanisms of purine rotation, either remaining inside or after flipping outside of the double helix. We study seven different sequences, which are neighbor variations of a well-studied A⋅T pair in A6-DNA. We observe the known effect of A⋅T steps favoring HG stability, and find evidence of triple-hydrogen-bonded neighbors hindering the inside transition. More importantly, we identify a dominant factor: the direction of the A rotation, with the 6-ring pointing either towards the longer or shorter segment of the chain, respectively relating to a lower or higher barrier. This highlights the role of DNA's relative flexibility as a modulator of the WCF/HG dynamic equilibrium. Additionally, we provide a robust methodology for future HG proclivity studies.

Introduction

In 1953, emblematic studies by Watson and Crick [1], and Franklin and Gosling [2] defined the structure of DNA for the first time. The discovery of a specific pairing of purine and pyrimidine bases via hydrogen bonds revealed not only DNA’s double-helical shape, but also the basis for its replication, which is essential for a genetic information carrier. In Fig 1A we depict a Watson-Crick-Franklin (WCF) A⋅T base pair, with the characteristic hydrogen-bonding pattern. Six years later, in 1959, Hoogsteen proposed an alternative base-pairing motif based on crystallographic data from A⋅T crystals, in which the purine base rolls 180° around the glycosidic bond [3], i.e. from anti to syn, with respect to the WCF geometry, as depicted in Fig 1B. In the Hoogsteen (HG) configuration, the pyrimidine forms hydrogen bonds with the 5-ring of the purine, rather than with its 6-ring, causing the opposite backbone C1’ atoms of the bases to be at a somewhat shorter distance, as well as some degree of twisting and bending of the double helix in the vicinity of the base pair [3, 4]. In the last decades, a number of studies have shown that the abundance of the HG conformation is non-negligible in canonical duplex DNA, and that its biological implications are very relevant. In 2011, Nikolova and coworkers reported the transient presence of HG base pairs in specific steps inside canonical duplex DNA using nuclear magnetic resonance (NMR) relaxation dispersion spectroscopy [5]. They reported populations of A⋅T and G⋅C HG base pairs of around 0.5%, with residence times of up to 1.5 ms. Recent studies have measured even larger HG populations, of 1.2%, in an A⋅T rich segment [6]. A few years later, Alvey and colleagues, also using NMR relaxation dispersion, demonstrated that HG base pairs appear in more diverse sequences than previously expected [7]. These studies shifted the perspective on HG base pairs, which were initially thought to appear mainly in distorted or damaged DNA, but are now considered to coexist with the WCF motif in a dynamic equilibrium [8]. A recent survey of structures of DNA-protein complexes in the Protein Data Bank (PDB) showed 140 HG base pairs out of a total of around 50,000 [4]. Some of these complexes are of particular biological relevance. For example, the human DNA polymerase-ι performs replication exclusively via HG base pairing [9, 10], a function that was previously thought to be unique of WCF base pairs. Other examples of the involvement of HG in DNA-protein complexes include the p53 tumor suppressor protein [11], the TATA-box binding protein involved in transcription [12], and the MATα2 homeodomain which regulates transcription in cells [13].

Fig 1

WCF-to-HG conformations and pathways.

(A) WCF conformation of the A16⋅T9 base pair in A6-DNA. (B) HG conformation of the same base pair. (C) Intermediate state of the WCF-to-HG transition via the inside pathway, i.e. with the rolling of A16 occurring within the double helix. (D) Intermediate state of the WCF-to-HG transition via the outside pathway. i.e. with A16 flipped out of the double helix via the major groove (E) Two cyclic paths in the space of the collective variables (CVs) χ′ and θ, which describe the rolling and the flipping of A16, respectively. Both paths start and end at the WCF state, describing a full revolution of A16 with intermediate states in which the 6-ring of A16 points toward the 5’ or the 3’ direction of DNA. The HG state is marked by an attractor, represented by the circle labeled aHG. The purple path samples the inside mechanism, while the green path explores the outside mechanism. Additional attractors, represented by the circles labeled a, are restrained between the WCF and the HG states along their respective paths and are also restrained along θ to prevent switching between the inside and the outside mechanisms. Restraints are represented by dashed lines.

The newly discovered relevance of HG base pairs, and the dynamic equilibrium in which they coexist with WCF base pairs, demands a more detailed understanding of the transition mechanisms between both conformations. Given the difficulty of observing the short-lived intermediate states experimentally, computational approaches based on molecular dynamics (MD) simulations, boosted with enhanced sampling, have provided valuable insights. In their original work from 2011, Nikolova and co-workers complemented their NMR results with conjugate peak refinement (CPR) simulations. They simulated the DNA sequence 5’-CGATTTTTTGGC-3’ (A6-DNA) in vacuo using the CHARMM27 force field [14], and studied the transition of the A16⋅T9 pair. Two collective variables (CVs) were used to steer the transition: (1) the χ glycosidic angle, which describes the rolling of A16 around the glycosidic bond and is defined by the atoms O4’-C1’-N9-C4; and (2) the θ base opening angle, which describes the flipping of A16 outside of the double helix towards the major groove and is defined in [15]. The CPR simulations revealed two types of pathways: one in which the rolling of the adenine, either clockwise or counterclockwise, occurs within the double helix, i.e. with a small base opening angle; and another in which the rolling occurs outside of the double helix, i.e. with a large base opening angle. We refer to these two types of pathways as inside and outside, respectively. Examples of intermediate states of the two types of pathways are shown in Fig 1C and 1D. The (χ, θ) pseudo-free-energy landscapes of Nikolova and coworkers showed an overwhelming preference for the inside mechanism, likely influenced by the lack of solvent to stabilize the flipped conformations. In 2015, Yang and colleagues used umbrella sampling (US) to obtain a (χ, θ) free-energy landscape of the same A16⋅T9 pair in A6-DNA [16]. They employed a version of the AMBER99-BSC0 force field [17] with modified parameters for the glycosidic torsion, and used explicit TIP3P water [18]. Their vast US calculations, with over 300 windows, revealed multiple pathways, including inside ones, and outside ones opening either towards the major or the minor groove. Yang and coworkers reported a free-energy difference from the WCF to the HG state of 4.4 kcal/mol, close to Nikolova and coworkers’ NMR result of 3.0 kcal/mol. The free-energy barrier was of 10–11 kcal/mol for the inside pathways, and up to 14 kcal/mol for the outside ones, again close to NMR results of 16 kcal/mol. A few years later, Yang and coworkers reported another (χ, θ) free-energy surface, this time for a base pair with an N1-methylated adenine [19]. Similarly, Ray and Andricioaei calculated a (χ, θ) free-energy surface for DNA and RNA [20]. While such 2D free-energy landscapes are rich in insight, they also demand significant computational expense to obtain several μs-long MD runs; hindering the systematic study of multiple DNA variations. In particular, investigations of the effect of the DNA sequence on HG base-pairing proclivity have remained beyond computational reach, even though there is experimental evidence that specific sequence patterns, e.g. A⋅T steps, can favor HG base pairing [21].

Path-based approaches offer an efficient alternative to study the WCF-to-HG transition. Techniques like path-metadynamics (PMD) [22]—i.e. metadynamics using as a unique CV the progress along an adaptive path connecting two known states in a multidimensional CV-space—focus sampling on the transition channels, rather than on the entire conformational landscape, e.g. the entire (χ, θ) plane. Recently, we performed a PMD mechanistic study [23], where we modeled the known A6-DNA with the state-of-the-art AMBER99-BSC1 force field [24] in explicit solvent. We focused on one direction of rotation, with the A16 6-ring pointing toward the A17 neighbor in the 3’ direction of DNA. We analyzed additional CVs, such as: the donor-acceptor distances of the characteristic hydrogen bonds of each type of base pairing, the distance between the backbone C1’ atoms, and the distance between the neighboring bases of the rolling A16 (C15 and A17). The latter CV displayed an increase at the mid-rotated intermediate state of inside pathways, indicating that the neighboring bases are displaced in order to accommodate the rotation of the adenine. In contrast, the distance between bases neighboring the rolling A16 remained unaffected in outside transitions, indicating that the inside pathways might be more subject to sequence dependence than the outside ones. We also found that, while the θ pseudo-dihedral as defined in [15] is effective at biasing base opening, the χ dihedral can induce rotations of the sugar ring, rather than of the adenine. This can be addressed by defining a less locally dependent pseudo-dihedral based on centers of mass, χ′, as done in [25]. Most importantly, and contrary to previous studies, we found a less pronounced preference for the inside pathway [23]. Together with other coauthors, we recently reported another pathway-focused study [25]. We used transition path sampling (TPS) [26] to generate unbiased trajectories, and transition interface sampling (TIS) [27] to calculate a WCF-to-HG free-energy difference of 3.2 kcal/mol, which compares well with Nikolova and coworker’s previous experiments [5]. Most importantly, TPS showed trajectories spontaneously changing from the inside to the outside mechanism, evidencing a preference for the latter one. This raises a debate with the previous results from free-energy surfaces, which favored inside pathways, although having either no explicit solvent [5] or restrained neighboring bases in their setup [16]. A predominance of outside pathways could explain why sequence variation does not seem to impact significantly the WCF-to-HG free-energy barrier, according to Alvey and coworkers’ NMR results [7]. Markov state modeling has also confirmed a dominant outside pathway [20].

In this work, we apply recent advances in PMD—which allow to treat multiple pathways in parallel—to study the WCF-to-HG inside and outside mechanisms in diverse DNA sequences with high efficiency. Unlike standard metadynamics [28], PMD does not suffer from exponential performance scaling with the number of CVs [29]. This is because the sampling is effectively done in only 1D, i.e. on the normalized progress component along the path, s, connecting the two known stable states in CV-space. This sampling along the path can be done with other well-established methods, like US, or using common algorithmic extensions, such as multiple walkers [23, 30]. The path curve is optimized based on the restrained cumulative sampling density, which is an estimator of the free-energy gradient, such that the method converges to the closest low free-energy path from the initial guess. However, this optimization can be challenging in systems with multiple pathways, especially when the transition can switch between mechanisms, as seen for the WCF-to-HG inside and outside pathways [23, 25]. To tackle this scenario, we recently developed multiple-walker multiple-path-metadynamics (MultiPMD) [31]. In this scheme, we initialize several paths, each with an associated group of walkers. Under normal conditions, all paths would converge to the same low free-energy channel. By introducing special walkers, i.e. repellers or attractors, the paths can be forced to diverge and find alternative mechanisms connecting the known states. Thus, with MultiPMD one can calculate the free-energy profiles along multiple pathways simultaneously, with sub-exponential performance scaling with the number of CVs, and exploiting the parallelism of current supercomputing resources.

Here, we apply MultiPMD to find inside and outside pathways, and free-energy profiles, for the WCF-to-HG transition of the A16⋅T9 pair in seven sequence variations based on the well-studied A6-DNA sequence. In order to verify the effect of local variations, we test four different nucleobases (A,T,G,C) as direct neighbors on the 5’ side of the rolling A16. We repeat the procedure for the neighboring base pair on the 3’ side. Our MultiPMD methods provide, for the first time, sufficient computational agility for a systematic evaluation of HG base pairing in multiple DNA chains. This enables us to extract trends about the structure and free energy of the base-pairing transition across sequences.

Materials and methods

Sequence variations

From [5] we take the original A6-DNA sequence 5’-CGATTTTTTGGC-3’, with its complementary strand 3’-GCTAAAAAACCG-5’, as shown in Fig A in S1 Text. We study the WCF-to-HG transition of the A16⋅T9 base pair, which is produced by the 180° rotation of A16 around its glycosidic bond. The original direct neighbors of A16 are C15 and A17, in the 5’ and in the 3’ direction respectively, which gives the local environment CAA. To test the local effect of the direct neighbor of A16 in the 3’ direction, we replace A17 with T, G and C. This yields the new local environments CAT, CAG and CAC. We apply the same criteria for the direct neighbor of A16 in the 5’ direction, and replace C15 with A, T and G. This gives the new local environments AAA, TAA and GAA. All sequences are generated as ideal B-DNA duplex structures using the make-na tools [32], and can be regenerated with the current W3DNA server [33]. Since make-na generates WCF base pairs, we manually rotate A16 to obtain the HG state for each sequence.

Molecular dynamics

System preparation is done with GROMACS 5.4.1 [34]. We employ the AMBER99-BSC1 force field [24] to model the DNA’s interatomic interactions. Each sequence is placed in a periodic dodecaedron box, with a distance of 1 nm between the DNA sequence and the edge of the box. Each system is then solvated in water, which is modeled by the TIP3P [18] force field. We add 25 mM of NaCl—modeled with the AMBER99 force field [35]—to mimic physiological conditions and neutralize the charge of each system. The systems are energy minimized at each stage of the preparation process. To run MD, the canonical sampling through velocity rescaling (CSVR) thermostat [36] is set at at 300 K, and the Parrinello-Rahman barostat [37] at 1 bar. Each sequence, both at the WCF and at the HG state, is equilibrated for 100 ns with a time step of 2 fs. From the equilibrations, we analyze key structural features and their variations from sequence to sequence.

Multiple-path-metadynamics

MultiPMD production runs are performed with GROMACS 5.4.1 [34], patched with PLUMED 2.3.1 [38] and with the added PMD code, available at: https://github.com/Ensing-Laboratory/PathCV. We use two CVs: the base opening pseudo-dihedral angle θ, as defined in [15], which has been repeatedly proven effective to bias base flipping [5, 16, 23]; and the base rolling pseudo-dihedral angle χ′, as defined in [23], which is a correction to the glycosidic torsion χ to prevent sugar rotation. The fitness of these two CVs has been discussed in [39]. A value of χ′ = 0 rad implies a mid-rotated A16 with its 6-ring pointing in the 3’ direction, while χ′ = ±π implies that the A16 6-ring points toward the 5’ direction. Negative values of θ imply base flipping toward the major groove. We do not consider base flipping toward the minor groove, since our previous TPS study showed spontaneous switching only from the inside to the major groove pathway [25]. To handle both directions of rotation of A16, the χ′ angle is treated with its sine and cosine. Then, the paths are curves in the CV-space spanned by [cos(χ′), sin(χ′), θ]. These paths are cyclic, starting and ending at the WCF configuration, whose coordinates in the [cos(χ′), sin(χ′), θ]-space are determined by averaging the CVs from the corresponding equilibration run. The initial guess for all paths describes a full revolution of A16, transitioning from WCF to HG to WCF, with no flipping. All outside paths are discretized as strings with 39 nodes, with the initial and the final node fixed at the same point in the [cos(χ′), sin(χ′), θ]-space, which marks the WCF state. The inside paths are discretized as strings of 11 nodes, since they require less flexibility. The progress component along the path, s, which usually grows from s = 0 to s = 1 from the initial to the final node, is now set to grow from s = −1 to s = +1, and to be periodic in the same range, in order to match the cyclic nature of the path. This implies that the WCF state corresponds to s = ±1. on the other hand, the HG state corresponds to s ≈ 0, but the exact value cannot be known a priori. Assigning the HG state to a specific fixed node in the middle of the path would require to make an assumption about the length of each section of the path. Instead, the HG state is marked by an attractor, i.e. a walker that does not participate in the free-energy calculation. The HG attractor is harmonically restrained, with a force constant of 50 kcal/mol, at the average value of χ′ and θ during the corresponding equilibration run. Values of −1 < s < 0 imply a rotation with the A16 6-ring in the 5’ direction, while values of 0 < s < +1 signify a rotation in the 3’ direction.

For each sequence, we initialize two paths, one for the inside, and one for the outside mechanism. We only consider the outside mechanism with opening toward the major groove, as previous work already demonstrated that opening towards the minor groove is unlikely [16, 25]. As we reported in [23], switching between the inside and the outside mechanisms can occur during a PMD calculation of base rolling. Since path-switching prevents convergence, we use a new strategy to keep the paths apart. The separation is induced by special walkers, called attractors, that do not take part in the free-energy calculation and guide the paths through specific intermediate states. We add two attractors on each path, which are steered to s = 0.5 and s = −0.5 by a moving harmonic restraint with a force constant of 5000 kcal/mol per squared path unit during the first 20 ps of the simulation, such that they are located at intermediate states in the WCF-to-HG and the HG-to-WCF sections of the cyclic path. These intermediate states between the two kinds of base pairing have the A16 mid-rotated, perpendicular to the other bases. To keep the inside path from flipping, its two attractors are restrained by a harmonic potential with a force constant of 5000 kcal/(mol rad2) at θ = 0.0, which prevents them from leaving the confines of the double helix. In turn, the repellers of the outside path are restrained by a harmonic potential with a force constant of 5000 kcal/(mol rad2) at θ = −π/2, which ensures that the A16 stays flipped toward the major groove at the mid-rotation. The large values for all the force constants are chosen because the attractors are located near the top of the free-energy barriers, and we require them to remain in these high-energy, mid-rotated, conformations. Then, each path has one attractor at the HG state, and two attractors at intermediate states; one at s = 0.5 and one at s = −0.5. In Fig 1E, we show a scheme of the inside and outside paths, the fixed nodes, the HG attractor and the intermediate-state attractors.

Nine standard walkers—for a total of 12 walkers per path—perform metadynamics on the s component of the path. The metadynamics Gaussian potentials have a width of 0.1 path progress units and a height of 0.05 kcal/mol, and are deposited every 1 ps. The paths are updated every 1 ps. The value of the half-life parameter is infinite for inside paths and 20 ps for outside paths. This parameter determines the flexibility of the path by setting the amount of simulation time that it takes for previous samples to weight only 50% of their original value for path updates. We use a tube potential—i.e. an upper harmonic wall at a distance of 0.0 on the component perpendicular to the path, z—with a force constant of 50 kcal/mol per path unit to maintain all walkers near the path. For all the sequences, we analyze the adapted paths after 7 ns of sampling. We verify that the dynamics have reached the “free-diffusion” regime along the biasing coordinate, indicating that the free-energy profile has been reasonably approximated. For all sequences, this occurs before 2 ns of biasing. We obtain metadynamics-based free-energy profiles—i.e., the negative of the sum of Gaussian potentials—every 0.1 ns, starting at 2 ns and finishing at 7 ns. We obtain our final free-energy profiles and error bars from the average and standard deviation of these metadynamics-based estimations. There are a few exceptions to the general protocol, which we detail in S1 Text.

All the data and PLUMED input files required to reproduce the results reported in this paper are available on PLUMED-NEST (www.plumed-nest.org), the public repository of the PLUMED consortium [40], as plumID:21.033.

Results and discussion

Stable-state structures

From the 100 ns equilibration runs at the WCF and HG stable states, we analyze structural variations between the different sequences. Table A in S1 Text presents the definitions, averages and standard deviations of several CVs during the equilibrations. See [23] and Fig C in S1 Text for more details about the CVs. First, we note that all equilibrations are stable at the respective base-pairing configurations, as evidenced by the characteristic hydrogen-bond distance, dWCF or dHG, and the conserved dHB with average values of ∼ 3 Å; and by the base rolling angle χ′ close to either +1.5 rad in WCF, or to -1.5 rad in HG. However, the standard deviations of dHG for the HG states, which range from 0.2 to 0.7 Å, already indicate that the stability of HG base pairing is not equal for all sequences. In contrast, the standard deviation of dWCF for WCF base pairs is of ∼ 0.1 Å, indicating more rigid hydrogen bonds. Another structural signature of the two kinds of base pairing is the distance between the C1’ atoms of A16 and T9, dCC, which shows an expected constriction from ∼ 10.6 in WCF, to ∼ 9.1 Å in HG conformations [4]. A key feature that varies significantly from sequence to sequence is the distance between the neighboring bases, dNB, with a range from 7.4 to 8.0 Å. In the following sections, we show how this CV relates to the free-energy barriers from the WCF-to-HG transition.

Free-energy differences and barriers

Inside and outside pathways, together with their corresponding free-energy profiles, are shown in Fig B in S1 Text. For all free-energy profiles, we report error bars. Additionally, we change the sampling time by 2 ns and show that the free energies do not change beyond the error bars (See dotted lines in Fig B in S1 Text). The attractors successfully keep the paths separated, with the inside ones near θ = 0 and the outside ones flipping toward the major groove. The larger error bars for the outside paths of sequences GAA and CAT indicate that, unlike the other sequences, their optimal outside transition channel might not intersect with the attractor and have a different base opening angle. Note that our attractors mark only two landmarks in a continuous spectrum of inside-to-outside pathways with different barriers, which overlap at the stable states. We do not relocate the attractors in order to maintain a direct comparison with the rest of the sequences. The free-energy profile for the CAA sequence compares well to our results from [23] and [39], which consider only the 3’ direction of rotation. To simplify the analysis of the numerous free-energy profiles, in Fig 2 we show only the barriers, in both the 3’ and the 5’ directions of rotation, as well the free-energy difference from WCF to HG. The free-energy differences calculated for the same sequence along inside and outside paths are consistent within ∼ 2 kcal/mol and have overlapping error bars; providing a sensible check for our profiles (see Fig 2B). The small inconsistencies are due to the implicit error of the metadynamics estimation and to the tube potentials, which restrain the sampling of the free-energy valleys in the direction perpendicular to the paths, thereby modifying the distributions with respect to unbiased sampling of the minima. As expected, WCF base pairing is preferred over HG in all cases. The free-energy difference between both states ranges from ∼ 0.5 to ∼ 6 kcal/mol across all sequences. This range agrees with free-energy differences obtained by NMR relaxation dispersion for other sequences [7]. According to our calculations, the sequence TAA presents the most favored HG state, with a free-energy difference of 0.6 to 2 kcal/mol with respect to WCF; and with the lowest barrier overall (9.7 kcal/mol) via the outside 3’ rotation. This result agrees with experimental reports of HG stability being increased by A⋅T steps [21]. The CAT sequence shows the most disfavored HG state, with a difference of 4.9 to 5.9 kcal/mol and the highest barrier overall (20 kcal/mol) via the 5’ inside rotation.

Fig 2

WCF-to-HG free-energy differences and barriers for the seven sequences.

The sequences are shown in Fig A in S1 Text. The free-energy differences and barriers are extracted from the average free-energy profiles and error bars in Fig B in S1 Text. Inside paths are shown in purple and outside paths are shown in green. (A) Free-energy barriers via the rotation with the 6-ring of A16 pointing in the 5’ direction. (B) Free-energy differences. (C) Free-energy barriers via the rotation with the 6-ring of A16 pointing in the 3’ direction.

For all sequences, the lowest free-energy barrier is that of the outside pathway with a 3’ rotation (see Fig 2C). Such a prevalent outside mechanism agrees with our previous results for the CAA sequence using TPS [25], as well as Markov state modeling by Ray and Andricioaei [20]. Most of the outside 3’ rotations also show a transition state closer to HG (s ≈ 0) than to WCF, which agrees with experiments [7]. These results imply that the conformational penalty of A16 rolling within the double helix is higher than that of it flipping out and back in.

Sequence dependence of the free-energy barriers

In [23], we showed how the inside rotation induces an increase in the distance between the neighbors of A16, which suggests that inside pathways have a stronger neighbor-dependence. Specifically, one may hypothesize that less flexible—i.e. double-ringed (A,G) or triple-hydrogen-bonded (G,C)—neighboring bases can hinder the inside rotation of A16. We observe such a trend in the barriers for the 3’ inside rotation (see Fig 2C). The 5’ neighbor variations with respect to CAA yield the following ranking of sequences, from the highest to the lowest barrier: GAA, CAA, TAA, AAA. While the 3’ neighbor variations yield the following order, again for the highest to the lowest barrier: CAG, CAC, CAT, CAA. These rankings would indicate that the barriers are increased mostly by G neighbors, followed by C, T and A. However, this trend is not repeated for the inside rotations in the 5’ direction (see Fig 2A), or for the outside rotations. Instead, we notice that the most dominant influence on the free-energy barrier is the direction in which the 6-ring of A16 points during the rotation. Remarkably, the rotation in the 3’ direction has a consistently and significantly lower barrier than in the 5’ direction (see Fig 2A and 2C). This is due to the asymmetrical length of the DNA sequence in each direction (see Fig A in S1 Text), i.e. the position of the rolling base along the sequence. Starting from the rotating base, A16, there are three base pairs in the 5’ direction and eight in the 3’ direction, which imply a difference in flexibility. We observe that, when the 6-ring rotates in the 5’ direction, the 5-ring protrudes in the opposite direction toward the longer segment of DNA, which is less flexible, causing a higher barrier. In contrast, when the 6-ring rotates in the 3’ direction, the 5-ring pushes toward the shorter, more flexible, segment of DNA, causing a lower barrier.

We analyze trends between the free-energy barriers and a few significant CVs. The average values of the CVs are taken from the attractors restrained at the intermediate states s = −0.5 and s = 0.5, which are close to the peaks of the free-energy barriers. In Fig 3A, we show the distance between the two neighbors of the rolling A16, dNB, in relation to the free-energy barrier. Several trends can be identified. First, dNB remains at low values (< 8 Å) for all outside paths, close to those of the stable states (see Table A in S1 Text); confirming that the transition with base flipping induces much less deformation of the neighbors. On the other hand, all the inside paths show larger values of dNB, ranging from ∼ 9 to ∼ 14 Å. Among the inside paths, two distinctive trends can be observed. Inside rotations in the 3’ direction reach much larger values of dNB, with a relatively low increase in the free-energy barrier (0.3 kcal/mol per Å). In contrast, for the inside rotations in the 5’ direction, the free-energy barrier quickly rises (1.7 kcal/mol per Å). This again highlights the role of the DNA segment’s relative length in each direction of the rolling base, as well as the favored 3’ direction of rotation.

Fig 3

Trends between the WCF-to-HG free-energy barriers and the average value of a few selected CVs.

The free-energy barriers are reported in Fig 2. The average values of the selected CVs are taken from the attractors restrained at s = −0.5 and s = 0.5. Inside paths are shown in purple and outside paths are shown in green. The barriers for the rotation with the A16 6-ring pointing in the 5’ direction are represented in solid triangles pointing up, while the ones in the 3’ direction are represented in outlined triangles pointing down. (A) Free-energy barrier vs. the distance between the neighboring bases of A16, dNB. The fit of the inside 5’ direction barriers is shown with a thick purple line (free-energy barrier = 1.7 kcal/(mol Å)dNB + 1.5 kcal/mol). The linear fit of the inside 3’ direction barriers is shown with a narrow purple line (free-energy barrier = 0.3 kcal/(mol Å)dNB + 12.6 kcal/mol). (B) Free-energy barrier vs. the number of water oxygens within 6 Å of the N6 atom of A16, Nwater. The linear fit of the inside 5’ and 3’ direction barriers is shown with a thick purple line (free-energy barrier = 0.2 kcal/mol dNB + 13.9 kcal/mol). (C) Free-energy barrier vs. the distance between the C1’ atoms of A16 and T9, dCC.

We also analyze the influence of the water solvation on the free-energy barriers (see Fig 3B). We measure the solvation using the parameter Nwater [25], i.e. the number of water oxygens within 6 Å of the N6 atom of A16, which is the atom involved in the conserved hydrogen bond of both WCF and HG base pairs. The intermediate states of the inside pathways show a wide range of Nwater values, from ∼ 8 to ∼ 21. For both the 3’ and the 5’ direction, an increase of Nwater in the intermediate state relates with an increase in the free-energy barrier (0.2 kcal/mol per water molecule). Rather than a direct solvation of A16, the increase of Nwater in the inside pathways is due to the separation of the neighboring nucleotides, as measured before by dNB, which generates an opening accessible to water. On the other hand, outside pathway intermediates present larger and mostly constant values of Nwater ≈ 24, with the flipped conformations of A16 being stabilized by the water solvent. Our reported ranges of Nwater for inside and outside mechanisms agree with those in [25].

Additionally, we study the relation of the free-energy barrier with the distance between the C1’ atoms of A16 and T9, dCC (see Fig 3C). Inside pathway intermediates show values of dCC mostly from ∼ 8 to ∼ 11 Å. This spans a larger range than that of the stable states (∼ 9.1 to ∼ 10.6 Å), indicating another possible conformational penalty for the inside pathways, but there is no clear correlation with the free-energy barrier. The intermediates of the outside pathway show larger dCC values, from ∼ 12 to ∼ 16 Å, which are expected given that the A16⋅T9 base pair is completely broken. There is no clear correlation of the free-energy barriers via the outside pathways with dCC, or with the other analyzed parameters.

Conclusion

Our efficient MultiPMD protocol allows, for the first time, a systematic study of HG base-pairing proclivity in diverse DNA chains. We investigate the rotation of the A16-T9 base pair in seven sequences, based on the previously studied A6-DNA [5, 23, 25]. The seven sequences are variations of the direct neighbor of the rolling A16 in the 3’ and the 5’ direction of DNA (see Fig A in S1 Text). For all sequences, we study the inside and outside pathways, i.e. without and with base flipping, as depicted in Fig 1. We only consider the outside pathway with flipping toward the major groove, because our previous TPS study showed spontaneous opening only in that direction [25]. However, opening toward the minor groove could be included in a future study by adding a third path and set of attractors. We obtain WCF-to-HG free-energy profiles for the inside and the outside pathway of each of the seven sequences (see Fig B in S1 Text). Our profiles are validated by: 1) the previous result for the CAA sequence reported in [39]; 2) the consistent WCF-to-HG free-energy difference between our inside and outside calculations; and 3) the relatively favored HG base pairing for the TAA sequence, which agrees with experimental evidence about the effect of A⋅T steps [21]. We run ∼ 7 ns with twelve walkers to calculate each free-energy profile; placing our total runtime to obtain the inside and outside profiles of one sequence at ∼ 168 ns. This requirement is well below the μs-long runs required for calculating previous free-energy surfaces [16, 19]. Our 1D free-energy profiles along the path progress, s, could be compared with 2D free-energy surfaces by a posteriori optimizing a string on the (χ′, θ) plane, as described in [41]. To this end, we have made a postprocessing tool for standard metadynamics free-energy surfaces available at: https://www.compchem.nl/software_package/trace-irc/.

From our free-energy calculations, we observe that all sequences have a preferred outside pathway, in which the A16 rotates with its 6-ring pointing in the 3’ direction (see Fig 2). This dominant outside pathway agrees with our previous results using TPS [25], as well as with published reports using Markov state models [20]. Our result is likely to settle the debate arising for previous simulations that showed favored inside pathways, but either with no explicit solvent [5], or with a modified force field and restrained neighboring bases upon initialization [16].

We analyze the possible influence of the varied direct neighbors of A16 on the free-energy barrier for its rotation. Based on the mechanistic analysis done in [23], we expect inside paths to be more sensitive to neighbor-dependence, since they require to increase the distance between the neighbors of A16, dNB, in order to accommodate the rotation. In Fig 2C, we observe that barriers for an inside 3’ rotation are increased mostly by G neighbors, followed by C, T and A. While this could point to triple-hydrogen-bonded neighbors hindering the transition, the trend is not reproduced as prominently for outside pathways. Instead, we observe that the most impactful factor for the free-energy barrier is the direction, either 3’ or 5’, in which the 6-ring of A16 points during the rotation. The barriers for the rotation in the 3’ direction are significantly lower than their counterparts. This difference is due to the asymmetric length of the DNA chain in each direction of the transitioning base pair. The rotation of the 6-ring in the 5’ direction causes the 5-ring to protrude in the 3’ direction, toward the longer and more rigid side of the DNA chain, causing a higher free-energy barrier. In contrast, the 3’ direction of rotation causes the 5-ring to push against the shorter and more flexible side of the DNA chain, which comes with a lower free-energy barrier. This trend is confirmed in Fig 3A. We observe that increasing the distance between the neighbors of A16, dNB, quickly raises the free-energy barrier for inside rotations in the 5’ direction, while the inside rotation in the 3’ direction can reach larger neighbor separations with a much lower free-energy penalty. We also observe that dNB is almost undisturbed with respect to stable-state values during the outside transitions. In Fig 3B, we analyze the number of water molecules surrounding A16, Nwater, which is large and mostly constant for all outside intermediates. The value of Nwater for inside intermediates is expectedly lower, and also evidences the energetically costly opening of the neighbors already described by dNB. Additionally, we analyze the distance between the C1’ atoms of the A16⋅T9 base pair, which separate significantly during the outside transitions. Nonetheless, we do not find a CV that correlates with the free-energy barrier of the outside pathways. In summary, our WCF-to-HG free-energy differences and barriers reveal that HG proclivity is weakly dependent on the direct neighboring nucleotides and strongly controlled by local flexibility. This provides a valuable filter for studies searching for HG base pairs—which might be mislabeled as WCF—in experimentally determined structures [4, 42]. Finding these base pairs, and their transitions pathways, is key to understand how, while WCF base pairs store static information, HG base pairs encode dynamic functions in living organisms [43].

We believe that this work provides a robust and efficient methodology for future investigations of HG base pairing in various sequences. This includes studies with equal number of base pairs on both sides of the transitioning base, in order to elucidate the role of the neighbors exclusively, without the effect of an asymmetric length. One could also study transitions of more than one base pair, similarly to Chakraborty and Wales discrete path sampling work [44]. More importantly, our observation about the dominant influence of the relative chain length, and associated flexibility, towards the 5’ and 3’ directions of the rolling base, has major mechanistic implications. Recent studies have observed HG base pairing in stressed, protein-bounded DNA [42]. Protein-DNA complexes that function via HG base pairs might not only recognize, but even induce the transition by modulating the rigidity of a DNA sequence. Our simulation protocol can also enable a fast investigation of protein-DNA complexes. Finally, the efficiency of the MultiPMD method also enables the use of higher levels of theory, such as hybrid quantum mechanics/molecular mechanics (QM/MM) studies, in which the HG transition of G⋅C base pairs could be simulated with flexible protonation states [45].

The supplementary text contains further details about the simulation protocol, collective variables and all supplementary figures and tables.

(PDF)

Click here for additional data file.

20 Dec 2021

Dear Dr. Ensing,

Thank you very much for submitting your manuscript "Sequence dependence of transient Hoogsteen base pairing in DNA" for consideration at PLOS Computational Biology.

As with all papers reviewed by the journal, your manuscript was reviewed by members of the editorial board and by several independent reviewers. In light of the reviews (below this email), we would like to invite the resubmission of a significantly-revised version that takes into account the reviewers' comments.

Reviewers raised serious concerns regarding the data reliability of the simulations. In the revised manuscript, the issue must be addressed in a satisfactory way. Reviewer 3 wonders about the implications for experimental biology. We realize that it is challenging to address this point in a study that addresses a fundamental question. However, we ask you to show biological relevance to concrete cases as much as possible.

We cannot make any decision about publication until we have seen the revised manuscript and your response to the reviewers' comments. Your revised manuscript is also likely to be sent to reviewers for further evaluation.

When you are ready to resubmit, please upload the following:

[1] A letter containing a detailed list of your responses to the review comments and a description of the changes you have made in the manuscript. Please note while forming your response, if your article is accepted, you may have the opportunity to make the peer review history publicly available. The record will include editor decision letters (with reviews) and your responses to reviewer comments. If eligible, we will contact you to opt in or out.

[2] Two versions of the revised manuscript: one with either highlights or tracked changes denoting where the text has been changed; the other a clean version (uploaded as the manuscript file).

Important additional instructions are given below your reviewer comments.

Please prepare and submit your revised manuscript within 60 days. If you anticipate any delay, please let us know the expected resubmission date by replying to this email. Please note that revised manuscripts received after the 60-day due date may require evaluation and peer review similar to newly submitted manuscripts.

Thank you again for your submission. We hope that our editorial process has been constructive so far, and we welcome your feedback at any time. Please don't hesitate to contact us if you have any questions or comments.

Sincerely,

Shi-Jie Chen

Associate Editor

PLOS Computational Biology

Nir Ben-Tal

Deputy Editor

PLOS Computational Biology

***********************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: Ensing and co-workers report on the pathways and free energies of Watson-Crick to Hoogsteen base pairing in A6 DNA with multiple sequence contexts in MultiPMD (multiple-path-metadynamics) simulations exploring two collective variables that enable the conversion. Overall the manuscript is well written and well summarizes the findings from the simulations. Of some concern is the short sampling time (< 10 ns per walker), results from single (sets of) simulations , and potential bias of attracters, repellers, and force constants on the results. For example, if simulations are run twice as long per walker or with twice as many walkers are consistent results obtained? In the conclusion, it is stated you need much less simulation time that the microsecond timescales seen previously – to better justify, perhaps present results that show convergence – for example, do independent runs on same sequence (with different arrangement of ions or initial velocities) give equivalent results? Figure S2 suggests poor convergence (or high variance) in the GAA and CAT sequences.

Minor points

- What ion parameters?

- (see Fig. S1, i.e. – needs closing “)”

Reviewer #2: Perez et al. reported the in- and out-HG pairing pathways using the adaptive path sampling method. Also, they investigated sequence dependent effects on those pathways.

In this computational study, they concluded that the out-path is more favorable than in-path. This is interesting. Overall, this manuscript was well written and the results are well presented.

I recommend this manuscript be published in the PLOS comp. Bio. after some minor revisions:

1. a.The authors need to describe in detail how the free energy profile along the paths were obtained.

b.I am wondering if these free energy profiles on their optimized paths be exactly reproduced by more

conventional US-WHAM or other non-equilibrium schemes along those paths? If they are, please give a proper

reference on this.

2. In Figure 2B, the free energy difference between WC and HG in various sequences were shown. In principle, this

free difference should be the same for the inside and outside paths. But non-negligible differences (up to several

kcal/mol) were found in some sequences (CAC, CAG, ….). This may raise a concern on the reliability of this data.

Also, in GAA, the free energy difference from the outside path (in green) has too much error. The authors should

provide some reasonable explanations for these two issues.

3. For the outside paths, only the major groove opening path (negative bf angle) was considered. (In this paper, the

bf angle (theta) is defined as negative for the major groove opening, while it was defined as positive in the original

paper). Although they claim that the minor groove path is unlikely, previous papers on full 2D free energy maps

(chi, theta) showed that the minor-groove path is quite possible. By matching the two

CVs (chi', theta) used in this study, re-projection of the previously reported (chi, theta) to (chi', theta) map would

not change the popularity of the minor groove path. The authors are invited to address this difference in

discussion.

4. In Figure S1, in the 3’ variation panel, there is a typo: CAG CAG �  CAT CAG?.

Reviewer #3: In this manuscript, the authors describe an enhanced sampling method, as applied to a base-flipping motion associated with the transition between canonical and Hoogsteen orientations. The primary conclusion is very intuitive: there is a smaller energetic barrier when the base first flips "out", rather than rearranging within the stacked bases of dsDNA. One certainly expects this result, since the sterics associated with base rotation in a confined environment will be extremely strong.

The manuscript is easy to read, and the work appears to be sound. However, the study is not suitable for PLOS Comp Bio. First, there is no clear connection to biology. That is, do the different barrier heights for different sequences explain a known biological phenomenon, or suggest new biological physics? As written, the results are primarily a display of the methodology. The second issue is that the general topic of enhanced sampling methods is not particular new. There are many methods available in the literature, and there is no comparison of the results with other methods. If this is a methods paper, then rigorous comparison should be provided. Comparisons would have to describe differences in accuracy and performance. This work would be more suitable for a simulation-methods-focused journal, such at JCTC.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: Yes

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

Figure Files:

While revising your submission, please upload your figure files to the Preflight Analysis and Conversion Engine (PACE) digital diagnostic tool, . PACE helps ensure that figures meet PLOS requirements. To use PACE, you must first register as a user. Then, login and navigate to the UPLOAD tab, where you will find detailed instructions on how to use the tool. If you encounter any issues or have any questions when using PACE, please email us at .

Data Requirements:

Please note that, as a condition of publication, PLOS' data policy requires that you make available all data used to draw the conclusions outlined in your manuscript. Data must be deposited in an appropriate repository, included within the body of the manuscript, or uploaded as supporting information. This includes all numerical values that were used to generate graphs, histograms etc.. For an example in PLOS Biology see here: http://www.plosbiology.org/article/info%3Adoi%2F10.1371%2Fjournal.pbio.1001908#s5.

Reproducibility:

To enhance the reproducibility of your results, we recommend that you deposit your laboratory protocols in protocols.io, where a protocol can be assigned its own identifier (DOI) such that it can be cited independently in the future. Additionally, PLOS ONE offers an option to publish peer-reviewed clinical study protocols. Read more information on sharing protocols at https://plos.org/protocols?utm_medium=editorial-email&utm_source=authorletters&utm_campaign=protocols

22 Mar 2022

Submitted filename: plos_reply.pdf

Click here for additional data file.

19 Apr 2022

Dear Dr. Ensing,

We are pleased to inform you that your manuscript 'Sequence dependence of transient Hoogsteen base pairing in DNA' has been provisionally accepted for publication in PLOS Computational Biology.

Before your manuscript can be formally accepted you will need to complete some formatting changes, which you will receive in a follow up email. A member of our team will be in touch with a set of requests.

Please note that your manuscript will not be scheduled for publication until you have made the required changes, so a swift response is appreciated.

IMPORTANT: The editorial review process is now complete. PLOS will only permit corrections to spelling, formatting or significant scientific errors from this point onwards. Requests for major changes, or any which affect the scientific understanding of your work, will cause delays to the publication date of your manuscript.

Should you, your institution's press office or the journal office choose to press release your paper, you will automatically be opted out of early publication. We ask that you notify us now if you or your institution is planning to press release the article. All press must be co-ordinated with PLOS.

Thank you again for supporting Open Access publishing; we are looking forward to publishing your work in PLOS Computational Biology.

Best regards,

Shi-Jie Chen

Associate Editor

PLOS Computational Biology

Nir Ben-Tal

Deputy Editor

PLOS Computational Biology

***********************************************************

Reviewer's Responses to Questions

Comments to the Authors:

Please note here if the review is uploaded as an attachment.

Reviewer #1: The authors largely addressed the concerns raised in previous review.

Reviewer #2: The authors have made minor revisions accordingly. This revised version appears to be sound and clear. I find this version of the manuscript to be suitable for the publication in the Plos Comp. Biol.

Reviewer #3: With the revised introduction and discussion, the context for the study is now clear. While the technical side is interesting, the revised focus on biological implications makes the study appropriate for PLOS CB.

**********

Have the authors made all data and (if applicable) computational code underlying the findings in their manuscript fully available?

The PLOS Data policy requires authors to make all data and code underlying the findings described in their manuscript fully available without restriction, with rare exception (please refer to the Data Availability Statement in the manuscript PDF file). The data and code should be provided as part of the manuscript or its supporting information, or deposited to a public repository. For example, in addition to summary statistics, the data points behind means, medians and variance measures should be available. If there are restrictions on publicly sharing data or code —e.g. participant privacy or use of data from a third party—those must be specified.

Reviewer #1: Yes

Reviewer #2: Yes

Reviewer #3: None

**********

PLOS authors have the option to publish the peer review history of their article (what does this mean?). If published, this will include your full peer review and any attached files.

If you choose “no”, your identity will remain anonymous but your review may still be made public.

Do you want your identity to be public for this peer review? For information about this choice, including consent withdrawal, please see our Privacy Policy.

Reviewer #1: No

Reviewer #2: No

Reviewer #3: No

20 May 2022

PCOMPBIOL-D-21-02065R1

Sequence dependence of transient Hoogsteen base pairing in DNA

Dear Dr Ensing,

I am pleased to inform you that your manuscript has been formally accepted for publication in PLOS Computational Biology. Your manuscript is now with our production department and you will be notified of the publication date in due course.

The corresponding author will soon be receiving a typeset proof for review, to ensure errors have not been introduced during production. Please review the PDF proof of your manuscript carefully, as this is the last chance to correct any errors. Please note that major changes, or those which affect the scientific understanding of the work, will likely cause delays to the publication date of your manuscript.

Soon after your final files are uploaded, unless you have opted out, the early version of your manuscript will be published online. The date of the early version will be your article's publication date. The final article will be published to the same URL, and all versions of the paper will be accessible to readers.

Thank you again for supporting PLOS Computational Biology and open-access publishing. We are looking forward to publishing your work!

With kind regards,

Anita Estes

PLOS Computational Biology | Carlyle House, Carlyle Road, Cambridge CB4 3DN | United Kingdom ploscompbiol@plos.org | Phone +44 (0) 1223-442824 | ploscompbiol.org | @PLOSCompBiol