Investigation of dynamical flexibility of D5SIC-DNAM inside DNA duplex in aqueous solution: a systematic classical MD approach

Abstract

Incorporation of artificial 3rd base pairs (unnatural base pairs, UBPs) has emerged as a fundamental technique in pursuit of expanding the genetic alphabet. 2,6-Dimethyl-2H-isoquiniline-1-thione: D5SIC (DS) and 2-methoxy-3-methylnaphthalene: DNAM (DN), a potential unnatural base pair (UBP) developed by Romesberg and colleagues, has been shown to have remarkable capability for replication within DNA. Crystal structures of a Taq polymerase/double-stranded DNA (ds-DNA) complex containing a DS–DN pair in the 3′ terminus showed a parallelly stacked geometry for the pre-insertion, and an intercalated geometry for the post-insertion structure. Unconventional orientations of DS–DN inside a DNA duplex have inspired scientists to investigate the conformational orientations and structural properties of UBP-incorporated DNA. In recent years, computational simulations have been used to investigate the geometry of DS–DN within the DNA duplex; nevertheless, unresolved questions persist owing to inconclusive findings. In this work, we investigate the structural and dynamical properties of DS and DN inside a ds-DNA strand in aqueous solution considering both short and long DNA templates using polarizable, and non-polarizable classical MD simulations. Flexible conformational change of UBP with major populations of Watson–Crick–Franklin (WCF) and three distinct non-Watson–Crick–Franklin (nWCFP1, nWCFP2, nWCFO) conformations through intra and inter-strand flipping have been observed. Our results suggest that a dynamical conformational change leads to the production of diffierent conformational distribution for the systems. Simulations with a short ds-DNA duplex suggest nWCF (P1 and O) as the predominant structures, whereas long ds-DNA duplex simulations indicate almost equal populations of WCF, nWCFP1, nWCFO. DS–DN in the terminal position is found to be more flexible with occasional mispairing and fraying. Overall, these results suggest flexibility and dynamical conformational change of the UBP as well as indicate varied conformational distribution irrespective of starting orientation of the UBP and length og DNA strand.

---

Introduction

In every DNA-based organism, genetic information is represented by a four-letter genetic alphabet composed of deoxyadenosine (dA), deoxyguanosine (dG), deoxycytidine (dC), and deoxythymidine (dT).¹ The storage and retrieval of this information depend on the formation of two base pairs, (d)A–dT/U and (d)G–(d)C. Synthetic biology,² which emerged over a century ago, aims to create new biological forms with potential applicability towards biomedical and bio-engineering fields. One promising approach to achieving this goal is to expand the amount of information that can be stored and retrieved in a cell.^3–5 As a result, scientists have dedicated considerable effort over the last decade to discovering fifth and sixth nucleotides that can form a third, unnatural base pair (UBP) with increased functionality that can be orthogonally replicated in DNA.^4,6–37 This would also expand the usefulness of nucleic acids for biological and biotechnological applications.

While several unnatural base pair candidates have been identified, only few have been shown to be able to be efficient for central dogma process,^24–37 among them, d5SICS (DS)–dNaM (DN) reported by Romesberg and co-workers, have been PCR amplified without sequence bias and efficiently transcribed in both directions.³² What makes the DS–DN UBP particularly interesting is that it forms an intercalated structure in duplex DNA and does not rely on complementary hydrogen bond formation for inter-strand pairing. The underlying cause of the unconventional orientations of DS–DN remains uncertain. It is yet to be verified whether the polymerase stabilizes the structure of the UBP-incorporated DNA, or if the UBP itself can be stabilized inside the DNA duplex.

Betz et al. reported several crystal structures of the large fragment of T. aquaticus (Taq) DNA polymerase with bound UBP-incorporated DNA duplex. Among these, for the structure in the pre-insertion phase, the UBP are observed to be arranged in a parallel stack (WCF orientation), whereas in the post-insertion structure, the DS–DN pair adopts an intercalated structure (nWCF) inside the DNA duplex.^32,38 Several research groups have used computational tools to investigate the structural properties of the UBP-incorporated DNA.^39–45 Datta et al. investigated the structure of DS–DN incorporated DNA through both QM and MD simulations and showed that the DS–DN distance is found to be consistent with a WCF pairing pattern during MD simulations.³⁹ Wetmore et al. considered UBP incorporated three nucleotides-long double strands and observed that for DS and DN, QM calculations suggest DS and DN adopt an intercalated nWCF structure, whereas a planar WCF-ike configuration has been predicted through MD simulations.⁴³ Barroso-Flores et al. have reported several conformers of DS and DN inside DNA duplex through their extensive MD simulations.^44,45 They concluded that an equilibrated structure of DS–DN incorporated DNA duplex may not have been achieved due to sampling time and/or forcefield incompatibility. In summary, unlike natural base pairs, UBPs appear to adopt different conformations in aqueous solution and during different steps of the replication processes as observed from X-Ray crystal structures and computational simulations.

Under these circumstances, we have divided our work into two parts. In the first part we investigate the UBP-incorporated DNA duplex in aqueous solution, aiming to investigate the inherent structural attributes of DS–DN within DNA duplex strands through classical MD simulations. Recognizing the significance of non-covalent interactions in UBP stabilization, we have carried out simulations with the multipolar/polarizable AMOEBA forcefield, alongside the fixed-charge non-polarizable AMBER force field. The organization of the paper is as follows: in the next section, we describe the development of forcefield parameters for UBP and simulations details. Next, we discuss the results of the investigation of dynamical structural properties of UBP-incorporated DNA by placing DS–DN in the middle of a dsDNA strand with different orientations with both short and long fragments of DNA. Subsequently, we describe the investigation of the same with UBP placed in a terminal position within the DNA, followed by concluding remarks.

Computational methods

DFT calculations

All gas phase geometry optimizations for the UBs (DS and DN) and UBP (DS–DN) have been performed using Gaussian 16 A.03⁴⁶ at the ωB97X-D⁴⁷/6-311++G(d,p)^48,49 level. Symmetry Adapted Perturbation Theory-DFT, SAPT(DFT)⁵⁰ analysis has been done using the PSI4 1.2 software package.⁵¹ NCIPLOT⁵² has been employed to investigate the topology of non-covalent interactions between DS and DN. To predict the stability of the UBP in the gas phase, we have calculated the interaction energy (IE_UBP) of the complexes by employing the following equation.

IE_UBP = E^opt_UBP − E^Frag_DS − E^Frag_DN

where E^opt_UBP = energy obtained from optimized UBP and E^Frag_UB = energy of the fragments (DS and DN) in the optimized geometry.

Molecular dynamics (MD) simulations

System setup. The DNA templates considered here are represented in Scheme 1. The UBP is placed inside the DNA duplex in two possible ways. In one case DS–DN is incorporated into the middle of a 9-mer DNA duplex designated as MUD (5′-GCGCDSGCGC-3′, Scheme 1). In the MUD structure, a DS–DN pair has been placed with different orientations. The non-Watson–Crick DNA models have been created through the placement of the UBP as intercalated forms, denoted MUD_SYN and MUD_ANTI. The parallel model, corresponding to a canonical Watson–Crick DNA duplex, is denoted MUD_PAR. Further, we have considered long DNA (MUDL) with 21 base pairs (Scheme 1) having UBP positioned in the middle of the DNA to investigate the size effect on the stability of the UBP incorporated DNA. We have also investigated the UBP incorporated DNA by placing the DS–DN in the 3′ terminal position, designated as UUD. It should be mentioned that here glycosidic orientation of the UBP has been considered based on the crystal structure (Fig. S1, ESI†).

Scheme 1 Schematic representation of the sequence of DS–DN incorporated DNA duplex of (A) MUD, (B) MUDL and (C) UUD. (D) depicts the 2D representations of DS and DN. E. Distance and angles associated with DNA.

General MD setup

AMOEBA. All simulations with the polarizable AMOEBA^53–57 (Atomic Multipole Optimized Energetics for Biomolecular) force field were performed with the TINKER HP software.⁵⁸ For the DS and DN AMOEBA parameters, we have used the parametrization tools available in the TINKER software in tandem with GDMA 2.3 for the atomic multipoles⁵⁹ for multipole generation (all parameters are provided in ESI†). The systems were built using the packmole⁶⁰ software. Initially, the UBP-incorporated DNA duplex complex was minimized using the BFGS nonlinear optimization algorithm with a convergence criterion (RMS gradient) of 0.1 Å. Subsequently, relaxation via MD in vacuum followed by implicit water with the GBSA model for 2 ns to obtain the starting system was performed. After that, the structure was solvated in explicit water in the center of a box with a volume 50 × 50 × 50 ų containing 24 000 water molecules and neutralized by adding Na⁺ using packmole. The system was heated to 300 K in 4 simulation steps (2 ns each) with an NVT ensemble removing all positional restraints (100.0–0.0 kcal Å⁻¹). After the equilibration step, MD simulations were carried out for 125 ns in an NPT ensemble (1 atm and 298 K) for 3 replicates each (total simulation time 375 ns). The Monte Carlo barostat and Bussi thermostat were used to maintain the pressure and temperature fixed respectively. The duration of the time step was 2 fs using RESPA integrator. The smooth particle mesh Ewald (PME) method⁶¹ was used in the calculation of charge, atomic multipole, and polarization interactions. A cutoff of 10 Å was used for van der Waals potential energy interactions and the real-space distance cutoff in the Ewald summation.⁶² For the analysis we have picked 10 000 frames from the entire 125 ns trajectory with equal intervals.

AMBER. The DS and DN parameters have been calculated with the PYRED program⁶³ to generate AMBER parameters.^64–66 The LEaP module⁶⁷ in AMBER20⁶⁸ was used to set up the simulation box with UBP-incorporated DNA duplex in water. Neutralization of the system with the required number of counterions (Na⁺), and solvation of the system in a cubic box filled with TIP3P water,⁶⁶ extending at least 12 Å from the DNA duplex was done with the LEAP module in AMBER. All MD simulations were performed with the AMBER20 pmemd.cuda program using the Ol15 AMBER force field.⁶⁵ Seven minimization steps were done with decreasing restraint (10.0–0.0 kcal mol⁻¹ Å⁻²) on the solute's heavy atoms. In each stage, the system was minimized with 5000 cycles of minimization of steepest descent, followed by 5000 cycles of conjugate gradient minimization. Subsequently, each system was heated to 300 K using Langevin dynamics^69,70 with a collision frequency of 2 ps⁻¹ followed by 7 ns of NVT equilibration with decreasing restraints (10.0–0.0 kcal mol⁻¹ Å⁻²) on the system's heavy atoms every ns. Production calculations for each system were performed for 1 μs in the NPT ensemble without restraints in triplicate–a total of 3 μs for each system. Total simulated time is 9 μs for all MUD structures, 3 μs for the MUDL structure, and 3 μs for the UUD structure. Long-range Coulomb interactions⁷¹ were handled with the smooth particle mesh Ewald method^61,62 using a 10 Å cutoff for real-space non-bonded interactions. For the analysis, we have picked 25 000 frames from the entire 1 μs trajectory with equal intervals.

MD analysis

Selected geometrical parameters associated with the UBP have been monitored including the UBP distance (d_DS–DN) (Scheme 1). In Scheme 1, we have shown how the distances and angles are calculated. We have also calculated (UB–NB) (NB = natural base), for DS the calculated distances are DC4-DS5 (d_4-5) and DG6-DS5 (d_6-5) whereas for DN, DC13-DN14 (d_13-14) and DG15-DN14 (d_15-14) are the calculated distances (Scheme 1). The distance between DC4-DG6 (d_4-6) is also calculated. The distance between the sulphur of DS and oxygen (–OMe) of DN is designated as d_O–S. Parameters related to angles have been calculated (<NB–UB–UB) to predict the conformational change of the UBP inside the DNA duplex. The measured angles include <DC4DS5DN14(a_4-5-14), <DG6DS5DN14(a_6-5-14) for DS and <DC13DN14DS5(a_13-14-5) and <DC15DN14DS5(a_15-14-5) for DN. Sampling, distribution plot and population analysis (Fig. 5) of the conformers (WCF and nWCF) have been carried based on DC4-DG6. Interstrand flipping has been detected based on UB–NB distances. d_O–S gives a general idea about the SYN and ANTI conformers; high O–S indicates ANTI orientation whereas SYN orientations are detected by low O–S. For every sampling visualization technique is also taken under consideration. The CPPTRAJ module⁷² in AMBER18 was used to analyze production dynamics, i.e., RMSD, RMSF and geometrical parameters. In addition, Python libraries NumPy,⁷² Matplotlib,⁷³ Pandas,⁷⁴ were also employed for further data processing and graphing.⁷⁵ Energy decomposition analysis (EDA) has been employed to investigate the intermolecular interactions between the UBP and residues of the rest of the systems. An in-house Fortan90-based EDA code was employed to calculate the nonbonded intermolecular interaction energies.⁷⁶ For EDA calculations we have picked 8000 frames from different trajectories for each conformer and then perform EDA calculations.

Results and discussion

DFT analysis

QM calculations have been carried out to investigate the possible geometries of a UBP comprised by DS and DN. DFT calculations indicate two possible intercalating conformers designated as SYN and ANTI through which DS and DN can interact with one another (Fig. 1). In the SYN conformer, the sulfur of DS and the methoxy group of DN are on the same side.

Fig. 1 Optimized structures of (A) SYN and (B) ANTI conformers of DS–DN at ωB97x-D/6-311++g(d,p) level. NCIPLOT of (C) SYN and (D) ANTI represented the non-covalent interactions between DS and DN in optimized structures.

The interaction energy calculations indicate that the SYN conformer (IE_UBP = −10.8 kcal mol⁻¹) is slightly more stable than the ANTI conformer (IE_UBP = −9.1 kcal mol⁻¹) (Table S1, ESI†).

Our DFT results are consistent with the UBP orientation observed in the post-insertion structure in Taq,^32,38 which show DS and DN form an intercalated structure inside the DNA duplex with the sulfur and methoxy group on the same side. SAPT analysis suggests that the dispersion component is the major contributing factor in total energy to stabilize the UBP as an intercalated structure (Table S1, ESI†). The NCI index analysis also shows non-covalent interaction between DS and DN in the intercalated structures (Fig. 1).

AMOEBA simulations

MD simulations with AMOEBA were performed for both the SYN and ANTI orientations of DS and DN for three replicates, spanning 125 ns each. It has been observed that UBP predominantly forms nWCF structures with occasional flanking and distortion. This distortion leads to the generation of WCF structures on a few occasions during the simulations.

It has been noticed from our simulation that nWCF orientations of the UBP are not static, rather the system explores different geometries. Smaller DS–DN distance generally indicates nWCF geometries of UBP whereas distorted and WCF structures are characterized through higher DS–DN distances (Fig. S2, ESI† and Table 1). Conformational changes of UBP are reflected in d_4-6 distances (Fig. 4); for WCF structures, d_4-6 shows reduced distances whereas nWCF and distorted structures display a range of d_4-6 values, which indicates that the flexibility of UBP has an impact on the adjacent BPs.

Table 1 Average DS–DN values with standard deviation for each system for each replicate. Distance values are in Å

System	DS–DN (AVG/STDV)
	Rep1	Rep2	Rep3
AMOEBA
MUD	10.74/1.23	10.31/1.90	11.43/1.69
AMBER
MUDSYN	8.68/1.20	8.46/0.89	9.48/1.04
MUDANTI	8.91/1.08	9.51/1.07	9.12/1.04
MUDPAR	9.90/1.00	9.88/1.02	9.84/1.02
MUDL	9.45/1.18	9.27/1.13	9.86/1.02

---

For nWCF orientations, a larger d_4-6 represents in-phase placement of UBP whereas a decreased d_4-6 denotes outer-phase orientations of UBP inside the DNA duplex. Two distinct structural transformational processes have been observed during the simulations that lead to the generation of several conformers (Fig. 2) (i) intra-strand flipping which transforms the geometry of the UBP from SYN to ANTI and vice versa, (ii) inter-strand flipping leading to reorient the DS and DN upside down (Fig. 3). SYN orientation can be recognized by high O–S value whereas low O–S distance generally depicts ANTI orientations. It should be mentioned that apart from SYN and ANTI, due to the dynamical movement of the DS and DN several other intermediates are also produced during the simulations. Inter-strand flipping can be recognized by the pattern shift of UB–NB distance, transitioning between high and low values (Fig. 4). Flexibility of the UBP is also synchronized with the RMSD values (Fig. S2, ESI† and Table 2); distorted UBP-incorporated DNA shows higher RMSD whereas sudden change of RMSD values implies structural transformation. Three sets of replicates have yielded varying conformer distributions, suggesting a stochastic arrangement of the conformers. Overall, AMOEBA force field simulations suggest a dynamical nature of the UBP intercalated ds-DNA, with multiple conformational orientations inside the DNA duplex, where nWCF-DNA structures exhibit the largest occurrence.

Fig. 2 Snapshot of different conformers i.e. (A) nWCFP1, (B) nWCFP2, (C) nWCFO, (D) WCF of UBP inside DNA duplex.

Fig. 3 Conformational change through intra- and inter-strand flipping observed during the MD simulations.

Fig. 4 (A) d_4-6 value: high value depicts nWCP conformers, low values depict WC and nWCO conformers, (B) DN–NB value: circles point out inter-strand flipping, (C) d_O–S values: low value generally indicates DS and DN are in same phase (SYN), high-value generally indicates DS and DN are in opposite phase (ANTI) for three replicates obtained from AMOEBA forcefield.

Table 2 Average RMSD values with standard deviation for each system for each replicate. RMSD values are in Å

System	RMSD (AVG/STDV)
	Rep1	Rep2	Rep3
AMOEBA
MUD	3.73/1.03	2.49/1.04	3.21/1.4
AMBER
MUDSYN	2.50/0.63	2.71/0.58	2.68/0.65
MUDANTI	2.42/0.62	2.53/0.66	2.44/0.48
MUDPAR	2.78/0.99	2.46/0.92	2.64/1.03
MUDL	4.49/1.13	4.49/1.19	4.56/1.09

---

AMOEBA provides an improved description of the non-bonded interactions due to the inclusion of permanent atomic multipoles and explicit polarization, albeit at a higher computational cost. Therefore, in order to enable longer sampling times, we have also investigated the same systems using AMBER simulations.

AMBER simulations

AMBER simulations have been carried out to further investigate the dynamical properties of MUD structure. Here we have considered MUD_SYN, MUD_ANTI and MUD_PAR conformations as the initial intercalated structures. For both cases simulations have been done for 3 replicates with simulation time of 1 μs each.

Structural analysis of MUD

Starting from the nWCF conformation with SYN orientation (MUD_SYN), the distribution plot and population analysis for each replicate indicate a prevalent population of nWCF conformers, with occasional occurrence of WCF structures (Fig. 6). Delving into the nWCF structures, it becomes evident that their orientations are not static; instead, the system explores various conformers with frequent transitions from one conformer to another specifically for replicate 2 and 3 (Table 3). Interestingly, dynamical conformational characteristics of UBP predicted by the AMBER forcefield align closely with the results derived from the AMOEBA simulations discussed above. With these frameworks, the majority of nWCF structures falls into three distinct categories: in-phase-intercalation1 (nWCFP1), in-phase-intercalation2 (nWCFP2), and outer-phase-intercalation (nWCFO) (Fig. 5). Calculated d_4-6 distances in nWCFP1 structures are found to be notably high (Fig. S4, ESI†), facilitating the accommodation of the intercalated UBP inside the DNA duplex. For nWCFP2, reduced d_4-6 has been observed and it becomes lowest for nWCFO insisting the UBP to settle at the outer phase of the DNA. Interestingly in these two cases, occasional distortion has been noticed leading to form mis-paired and flanked structures. During the simulation, SYN-to-ANTI transformation or vice versa through intra-strand flipping is witnessed whereas DS–DN are found to be upside down in their position through inter-strand flipping. In the distribution plot, flipping has been pointed out. SYN conformers can be recognized by shorter d_O–S whereas higher d_O–S represent ANTI and distorted orientations of UBP (Fig. S4, ESI†).

Table 3 Transition count of the conformational change obtained from AMBER simulations

Systems	Conformational transition count
	Rep1	Rep2	Rep3
MUDSYN	2	1150	979
MUDANTI	589	662	407
MUDPAR	757	757	797
MUDL	550	435	316

---

Fig. 5 RMSF values for all the systems obtained from AMBER simulation.

Generation of several conformers with different orientations suggests flexibility of the UBP inside the DNA duplex. Higher RMSF for DS and DN further confirms the flexible nature of the UBP (Fig. 5). RMSD values of the entire DNA are found to be synchronized with conformational orientations of the UBP; analogous to what is observed with the AMOEBA force field. Here also elevated RMSD values correlate with distorted structures whereas nWCFP1 and WCF structures exhibit comparatively lower RMSD values (Fig. S3, ESI†). It is noticed that the conformational change of the UBP is reflected on the associated geometrical parameters related to the UBP depicted in Fig. S4 (ESI†).

The nWCFP1 structures, characterized by intercalation, exhibit the smallest DS–DN distances (d_DS–DN). By contrast, the WCF structure, akin to the natural BP orientations, displays an evident increase in distance. Notably, the nWCFP2 and nWCFPO structures reveal fluctuating d_DS–DN values, suggesting the formation of distorted intercalated arrangements. Overall, AMBER simulations are in agreement with the AMOEBA simulation which suggests that flexibility and dynamical conformational change of the DS–DN are not an artifact of the force fields; rather it is a features of the UBP-incorporated DNA.

To study the impact of the starting conformation on the geometry of UBP we explored MUD_ANTI and MUD_PAR structures as starting points for the simulations. Similar to the MUD_SYN scenario, commencing with MUD_ANTI also revealed frequent conformational change from one conformer to another (Table 3). Here also distribution plot and population diagram (Fig. 6) suggest that WCF and nWCF structures are generated during the simulations where nWCF structures are found to be predominant. Noticeably, in this case outer-phase nWCFO structures are not observed. Conformational change between SYN and ANTI further verifies the occurrence of intra-strand flipping, a phenomenon evident through the corresponding O–S distances (Fig. S5, ESI†). Conformational change through Inter-strand flipping is also discernible from the UB–NB distance curves (Fig. S5, ESI†), which is also pointed out in the distribution curve. Here the pattern of the RMSF (Fig. 5) looks similar to the one obtained for MUD_SYN, underscoring the flexibility of the UBP, which remains dynamic and doesn’t depend on the initial structure. The stability of WCF and nWCFP1 structures are confirmed by low RMSD values whereas high RMSD value of distorted structures indicates that they are comparatively less stable (Fig. S3, ESI†).

Fig. 6 Distribution and population of the conformers for all the replications for (A) MUD_SYN, (B) MUD_ANTI, (C) MUD_PAR and (D) MUDL. WC, nWCP1, nWCP2, nWCO are designated by Green, Purple, Brown and Pink respectively. Black pointer indicates inter-strand flipping with the change of SYN to ANTI and vice versa. Red pointer indicates inter-strand flipping with no conformational change. Blue pointer indicates intra-strand flipping.

In the context of MUD_PAR, the distribution plot suggests the frequent generation of both WCF and nWCF structures across all the replicates as observed from the distribution plot (Fig. 6). Notably, the WCF structure exhibits a population exceeding 40%, signifying a higher prevalence than the other scenarios as observed in population analysis (Fig. 6). In this case alongside in-phase nWCF (nWCFP1) structures, outer-phase (nWCFO) are also generated which include occasionally distorted structures. Interestingly in this case no distinct intra and inter-flipping processes have been noticed during the simulation. Flexibility of the UBP is evident from the RMSF plot (Fig. 5) whereas dynamical conformational change has been reflected in the distribution plot (Fig. 6) and conformational transition count value (Table 2). High RMSD values are observed in the region of nWCFO structures indicating the generation of UBP-distortion mainly in the outer-phase region. Here change of pattern of UB–NB distance and <UB–NB–NB angle curve depict the transition between nWCF(P/O) and WCF conformers which is also reflected in O–S distance (Fig. S6, ESI†).

Collectively, simulations conducted using both the AMBER and AMOEBA force fields consistently highlight the flexible and dynamical conformational changes of UBP inside DNA duplex which leads to generating both WCF and nWCF structures. It is also evident from the simulation that UBP has an inclined tendency to stay as a nWCF forms throughout the simulations.

We have used energy decomposition analysis (EDA) to investigate the interactions between the base pairs using the WCF, nWCFP1, nWCFP2, and nWCFO structures (Table 4). We have examined the interactions between the DS–DN as well as adjacent complimentary DC–DG base pairs (cBPs) (Table 1). It is observed that for nWCFP1, DS–DN is stabilized through vdW interactions, where the vdW energy is −11.3 kcal mol⁻¹. cBPs are stabilized through coulomb interactions with associated E_Coul = ∼−9.0 kcal mol⁻¹, indicating the stability of the base pairs. nWCFP2 structures also show similar interactions between DS–DN as well as cBPs as obtained from EDA analysis. During the calculation of energy decomposition analysis (EDA) for WCF structures, we have identified two distinct interaction regions between cBPs. In one instance (WCF1), Coulombic interactions are approximately around ∼−7.3 kcal mol⁻¹, while in the alternative scenario (WCF2), they have escalated to around ∼−8.0 kcal mol⁻¹. The calculated van der Waals interaction energy, E_vdW, for WCF structures is reduced to ∼−1.5 kcal mol⁻¹ between DS and DN. We have also calculated the interaction energies between DS and DN as well as adjacent DG–DC base pairs for nWCFO structures.

Table 4 Non-covalent interactions value along with the standard deviations between BPs and UBP for short DNA (MUD). E_col (Coulamb energy), E_vdW (van daar Wal energy) are in kcal mol⁻¹

cBPs	E Coul/stdv
	nWCFP1	nWCFP2	nWCFO	WCF1	WCF2
DC3-DG16	−9.2/2.7	−9.7/2.5	−7.0/2.6	−7.2/	−7.9
DC4-DG15	−9.0/2.6	−8.9/2.6	−6.3/3.1	−7.3	−7.5
DG6-DC13	−9.2/2.7	−8.0/2.5	−7.2/2.6	−7.5	−8.9
DC7-DG12	−9.2/2.4	−9.4/2.7	−7.2/2.8	−7.5	−9.1

---

UBP	E vdW/stdv
	nWCFP1	nWCFP2	nWCFO	WCF1	WCF2
DS–DN	−11.3/1.5	−10.9/2.7	−6.4/2.9	−1.5	−1.7

---

Here the Coulomb interactions between cBPs are further decreased to E_Coul ∼−7.0 kcal mol⁻¹ along with a significant reduction of E_vdW (−6.4 kcal mol⁻¹) for DS–DN. Overall, the population of the different conformers is directly synchronized with the UBP and neighboring cBPs interactions where nWCF which correspond to the highest populations (39.4%) shows the largest interactions as obtained from the EDA analysis.

MUDL. A system comprising ds-DNA with 21 base pairs with the placement of the UBP at the middle of the DNA duplex has also been considered to investigate the impact of a larger strand on UBP conformational stability (Scheme 1). Our results suggest the frequent generation of both WCF and nWCF structures during the simulations to different extents for different replicates as observed from the distribution plot (Fig. 6) and conformational transition count (Table 3). Both in-phase and outer-phase (nWCFP1 and nWCFO) nWCF structures are generated where the distribution and population of the conformers appear different for the replicates. The d_O–S curve (Fig. S7, ESI†) predicts a majority of the time the system is in the SYN form during the simulation, whereas occasional distorted structure has been generated with larger d_O–S. It is found that RMSD values are comparatively higher than that of short DNA (Fig. S3, ESI†), which indicates that flexibility of the UBP transfers to the entire DNA duplex systems, which is further reflected in high RMSF values of all the nucleotides (Fig. 5). RMSF values of the DS and DN are also found to be higher for MUDL than other MUD systems which implies that larger DNA duplex is able to give more freedom to the UBP to maintain and amplify its flexible nature. Interestingly like short DNA, here also the RMSF shows similar patterns indicating despite high flexibility, the characteristic of the UBP- incorporated DNA remains similar. High fluctuations of the UBP are also reflected in the high range of DS–DN distance as observed in Fig. S3 (ESI†).

EDA analysis has also been employed to predict the interactions between DS and DN as well as complimentary natural base pairs for both WCF and nWCF (nWCFP1 and nWCFO) conformers. vdW energy between DS–DN is found to be the highest for nWCFP1 followed by nWCFO and lowest for WCF. Interestingly, here the interactions between complementary base pairs are similar irrespective of their orientations (Table 5). This suggests that as the length of the DNA increases, flexibility and different conformational orientations do not influence the interactions between the adjacent cBPs. Consequently, unlike MUD, the average population of the conformers for MUDL is almost equal. It further suggests that instead of UBP stability, interaction with adjacent base pairs is a major factor for the conformational distribution and population of the conformations.

Table 5 Non-covalent interactions along with the standard deviations between BPs and UBP between BPs and UBP for long DNA (MUDL). E_col (Coulamb energy), E_vdW (van daar Wal energy) are in kcal mol⁻¹

cBPs	E Coul/stdv
	nWCFP1	nWCFO	WCF
DG9-DG34	−9.1/2.6	−9.2/3.0	−9.3/2.5
DC10-DG33	−8.5/2.6	−8.2/2.7	−9.1/2.7
DG12-DC31	−8.5/2.6	−9.1/2.8	−8.9/2.7
DC13-DG30	−9.4/2.6	−9.4/3.0	−9/2.7
E Coul/stdv
UBP	nWCFP1	nWCFO	WCF
DS–DN	−9.0/2.1	−7.7/2.3	−1.7/0.9

---

UUD. Here DS and DN are positioned in the 3′ terminus of the DNA to study the dynamical properties of the UBP-incorporated DNA in solution employing AMBER force fields for 1 μs for three replicates. It has been observed that the UBP in a terminal position is more flexible in nature as observed from the RMSF values (Fig. S8, ESI†). This flexibility leads to form frequent mis-paired and frayed structures of the UBP (Fig. 7). RMSD values also suggest fluctuating characteristics of UUD structures. Notably, in this case nWCF structures are found to be predominant when the UBP is in the terminal position. Here also conformational change has been observed through both intra and inter-strand flipping. From the EDA analysis, it has been noticed that the Coulomb interaction between adjacent cBPs are smaller compared with the internal UBP systems, indicating flexibility of the UBP also perturb the neighbouring cBPs. It has been observed that DS and DN are stabilized through vdW interactions with associated E_vdW = −6.5 kcal mol⁻¹ indicating the interaction is significantly reduced than what is observed in nWCFP1 structures.

Fig. 7 Snapshots of different geometries of UUD form of DNA duplex during simulations through intra and inter-strand flipping.

Taken together, our results provide a picture of the effects of UBP-incorporated DNA in pure water. It should be noted that these effects may be different depending on several factors such as sequence context, salt concentration, and/or the orientation of glycosidic bonds. These effects are beyond the scope of the present work and will be investigated in the future.

Conclusions

We investigated structural aspects, conformational changes, and stability of DS–DN incorporated DG–DC base pair rich DNA duplex by considering both short and long forms of DNA duplex simulated with both polarizable AMOEBA and AMBER force fields. It was found from both force fields that unlike natural base pairs, UBP can persist as both WCF and nWCF conformers inside DNA duplex with flexible orientations and dynamical conformational change, which agrees with the previous simulated results^39,44 with AMBER and CHARMM forcefields. This could suggest that the flexibility of UBP inside DNA is not an artifact of the forcefields, but rather it is an intrinsic property of this DS–DN incorporated DNA. However, our simulation still may not explore the complete conformational surface due to insufficient sampling time. It is evident from our analysis that conformational orientation perturbs the stability of the neighboring cBP mostly for shorter fragments of DNA, which reflects on the populations of the conformers. In shorter DNA strands, nWCF conformers are predominant whereas equal distributions are noticed for long-DNA. Our simulated results also suggest the fluctuating nature of DS and DN in the terminal position.

Conflicts of interest

There are no conflicts to declare.

Acknowledgements

This work was supported by NIH R01GM108583 and Fulbright-Nehru postdoctoral Scholarship. Computing time from UNT CASCaM funded by NSF Grant No. CHE1531468 and OAC-2117247 is gratefully acknowledged. Tanay Debnath is a USIEF postdoctoral fellow (Fulbright-Nehru postdoctoral Scholarship). Computational time for this project was provided by the University of North Texas CASCaM CRUNTCh3 and 4 high-performance cluster and the University of Texas at Dallas, Cyberinfrastructure and Research Services, Ganymede and Titan HPC clusters and NCSA DELTA at University of Illinois Urbana-Champaign. The authors thank the reviewers for their constructive comments, which have helped improve our manuscript.

References

1. J. D. Watson and F. H. C. Crick, Nature, 1953, 171, 737–738.
2. T.-C. Tang, B. An, Y. Huang, S. Vasikaran, Y. Wang, X. Jiang, T. K. Lu and C. Zhong, Nat. Rev. Mater., 2021, 6, 332–350.
3. X. Yan, X. Liu, C. Zhao and G.-Q. Chen, Signal Transduction Targeted Ther., 2023, 8, 199.
4. A. Ambrogelly, S. Palioura and D. Söll, Nat. Chem. Biol., 2007, 3, 29–35.
5. E. C. Fischer, K. Hashimoto, Y. Zhang, A. W. Feldman, V. T. Dien, R. J. Karadeema, R. Adhikary, M. P. Ledbetter, R. Krishnamurthy and F. E. Romesberg, Nat. Chem. Biol., 2020, 16, 570–576.
6. M. Kimoto and I. Hirao, Chem. Soc. Rev., 2020, 49, 7602–7626.
7. J. Chelliserrykattil, H. Lu, A. H. F. Lee and E. T. Kool, ChemBioChem, 2008, 9, 2976–2980.
8. C. Kaul, M. Müller, M. Wagner, S. Schneider and T. Carell, Nat. Chem., 2011, 3, 794–800.
9. E. Meggers, P. L. Holland, W. B. Tolman, F. E. Romesberg and P. G. Schultz, J. Am. Chem. Soc., 2000, 122, 10714–10715.
10. N. Tarashima, Y. Komatsu, K. Furukawa and N. Minakawa, Chem. – Eur. J., 2015, 21, 10688–10695.
11. D. A. Malyshev, Y. J. Seo, P. Ordoukhanian and F. E. Romesberg, J. Am. Chem. Soc., 2009, 131, 14620–14621.
12. T. Lavergne, D. A. Malyshev and F. E. Romesberg, Chem. – Eur. J., 2012, 18, 1231–1239.
13. R. Yamashige, M. Kimoto, Y. Takezawa, A. Sato, T. Mitsui, S. Yokoyama and I. Hirao, Nucleic Acids Res., 2012, 40, 2793–2806.
14. Z. Yang, F. Chen, J. B. Alvarado and S. A. Benner, J. Am. Chem. Soc., 2011, 133, 15105–15112.
15. K. Hashimoto, E. C. Fischer and F. E. Romesberg, J. Am. Chem. Soc., 2021, 143, 8603–8607.
16. S. Hoshika, N. A. Leal, M.-J. Kim, M.-S. Kim, N. B. Karalkar, H.-J. Kim, A. M. Bates, N. E. Watkins, H. A. SantaLucia, A. J. Meyer, S. DasGupta, J. A. Piccirilli, A. D. Ellington, J. SantaLucia, M. M. Georgiadis and S. A. Benner, Science, 2019, 363, 884–887.
17. H. Wang, L. Wang, N. Ma, W. Zhu, B. Huo, A. Zhu and L. Li, ACS Synth. Biol., 2022, 11, 334–342.
18. T. P. Hettinger, Proc. Natl. Acad. Sci. U. S. A., 2017, 114, E6476–E6477.
19. N. R. Jena and P. Das, J. Biomol. Struct. Dyn., 2023, 41, 366–376.
20. I. Hirao, M. Kimoto and R. Yamashige, Acc. Chem. Res., 2012, 45, 2055–2065.
21. M. M. Georgiadis, I. Singh, W. F. Kellett, S. Hoshika, S. A. Benner and N. G. J. Richards, J. Am. Chem. Soc., 2015, 137, 6947–6955.
22. A. Marx and K. Betz, Chem. – Eur. J., 2020, 26, 3446–3463.
23. R. Dörrenhaus, P. K. Wagner and S. Kath-Schorr, Biol. Chem., 2023, 404, 883–896.
24. A. W. Feldman and F. E. Romesberg, Acc. Chem. Res., 2018, 51, 394–403.
25. D. A. Malyshev, K. Dhami, H. T. Quach, T. Lavergne, P. Ordoukhanian, A. Torkamani and F. E. Romesberg, Proc. Natl. Acad. Sci. U. S. A., 2012, 109, 12005–12010.
26. S. A. Mukba, P. K. Vlasov, P. M. Kolosov, E. Y. Shuvalova, T. V. Egorova and E. Z. Alkalaeva, Mol. Biol., 2020, 54, 475–484.
27. K. H. Lee, K. Hamashima, M. Kimoto and I. Hirao, Curr. Opin. Biotechnol, 2018, 51, 8–15.
28. A. W. Feldman, V. T. Dien, R. J. Karadeema, E. C. Fischer, Y. You, B. A. Anderson, R. Krishnamurthy, J. S. Chen, L. Li and F. E. Romesberg, J. Am. Chem. Soc., 2019, 141, 10644–10653.
29. K. Hamashima, M. Kimoto and I. Hirao, Curr. Opin. Chem. Biol., 2018, 46, 108–114.
30. K. Morihiro, Y. Moriyama, Y. Nemoto, H. Osumi and A. Okamoto, J. Am. Chem. Soc., 2021, 143, 14207–14217.
31. D. A. Malyshev and F. E. Romesberg, Angew. Chem., Int. Ed., 2015, 54, 11930–11944.
32. K. Betz, D. A. Malyshev, T. Lavergne, W. Welte, K. Diederichs, T. J. Dwyer, P. Ordoukhanian, F. E. Romesberg and A. Marx, Nat. Chem. Biol., 2012, 8, 612–614.
33. K. Duffy, S. Arangundy-Franklin and P. Holliger, BMC Biol., 2020, 18, 112.
34. V. T. Dien, M. Holcomb, A. W. Feldman, E. C. Fischer, T. J. Dwyer and F. E. Romesberg, J. Am. Chem. Soc., 2018, 140, 16115–16123.
35. D. A. Malyshev, K. Dhami, T. Lavergne, T. Chen, N. Dai, J. M. Foster, I. R. Corrêa and F. E. Romesberg, Nature, 2014, 509, 385–388.
36. N. G. J. Richards and M. M. Georgiadis, Acc. Chem. Res., 2017, 50, 1375–1382.
37. J. Oh, J. Shin, I. C. Unarta, W. Wang, A. W. Feldman, R. J. Karadeema, L. Xu, J. Xu, J. Chong, R. Krishnamurthy, X. Huang, F. E. Romesberg and D. Wang, Nat. Chem. Biol., 2021, 17, 906–914.
38. K. Betz, D. A. Malyshev, T. Lavergne, W. Welte, K. Diederichs, F. E. Romesberg and A. Marx, J. Am. Chem. Soc., 2013, 135, 18637–18643.
39. S. Jahiruddin and A. Datta, J. Phys. Chem. B, 2015, 119, 5839–5845.
40. L. Eberlein, F. R. Beierlein, N. J. R. van Eikema Hommes, A. Radadiya, J. Heil, S. A. Benner, T. Clark, S. M. Kast and N. G. J. Richards, J. Chem. Theory Comput., 2020, 16, 2766–2777.
41. Q. Wang, X.-Y. Xie, J. Han and G. Cui, J. Phys. Chem. B, 2017, 121, 10467–10478.
42. S. Jahiruddin, N. Mandal and A. Datta, ChemPhysChem, 2018, 19, 67–74.
43. I. Negi, P. Kathuria, P. Sharma and S. D. Wetmore, Phys. Chem. Chem. Phys., 2017, 19, 16365–16374.
44. R. Galindo-Murillo and J. Barroso-Flores, Phys. Chem. Chem. Phys., 2017, 19, 10571–10580.
45. R. Galindo-Murillo and J. Barroso-Flores, J. Biomol. Struct. Dyn., 2020, 38, 4098–4106.
46. M. J. Frisch, G. W. Trucks, H. B. Schlegel, G. E. Scuseria, M. A. Robb, J. R. Cheeseman, G. Scalmani, V. Barone, G. A. Petersson, H. Nakatsuji, X. Li, M. Caricato, A. V. Marenich, J. Bloino, B. G. Janesko, R. Gomperts, B. Mennucci, H. P. Hratchian, J. V. Ortiz, A. F. Izmaylov, J. L. Sonnenberg, D. Williams-Young, F. Ding, F. Lipparini, F. Egidi, J. Goings, B. Peng, A. Petrone, T. Henderson, D. Ranasinghe, V. G. Zakrzewski, J. Gao, N. Rega, G. Zheng, W. Liang, M. Hada, M. Ehara, K. Toyota, R. Fukuda, J. Hasegawa, M. Ishida, T. Nakajima, Y. Honda, O. Kitao, H. Nakai, T. Vreven, K. Throssell, J. A. Montgomery, Jr., J. E. Peralta, F. Ogliaro, M. J. Bearpark, J. J. Heyd, E. N. Brothers, K. N. Kudin, V. N. Staroverov, T. A. Keith, R. Kobayashi, J. Normand, K. Raghavachari, A. P. Rendell, J. C. Burant, S. S. Iyengar, J. Tomasi, M. Cossi, J. M. Millam, M. Klene, C. Adamo, R. Cammi, J. W. Ochterski, R. L. Martin, K. Morokuma, O. Farkas, J. B. Foresman and D. J. Fox, in Gaussian 16 Rev. A.03, Vol, Wallingford, CT, 2016.
47. J.-D. Chai and M. Head-Gordon, Phys. Chem. Chem. Phys., 2008, 10, 6615–6620.
48. R. Krishnan, J. S. Binkley, R. Seeger and J. A. Pople, J. Chem. Phys., 2008, 72, 650–654.
49. A. D. McLean and G. S. Chandler, J. Chem. Phys., 2008, 72, 5639–5648.
50. A. J. Misquitta, R. Podeszwa, B. Jeziorski and K. Szalewicz, J. Chem. Phys., 2005, 123(21), 214103.
51. R. M. Parrish, L. A. Burns, D. G. A. Smith, A. C. Simmonett, A. E. DePrince, III, E. G. Hohenstein, U. Bozkaya, A. Y. Sokolov, R. Di Remigio, R. M. Richard, J. F. Gonthier, A. M. James, H. R. McAlexander, A. Kumar, M. Saitow, X. Wang, B. P. Pritchard, P. Verma, H. F. Schaefer, III, K. Patkowski, R. A. King, E. F. Valeev, F. A. Evangelista, J. M. Turney, T. D. Crawford and C. D. Sherrill, J. Chem. Theory Comput., 2017, 13, 3185–3197.
52. J. Contreras-García, E. R. Johnson, S. Keinan, R. Chaudret, J.-P. Piquemal, D. N. Beratan and W. Yang, J. Chem. Theory Comput., 2011, 7, 625–632.
53. A. Grossfield, P. Ren and J. W. Ponder, J. Am. Chem. Soc., 2003, 125, 15671–15682.
54. P. Ren and J. W. Ponder, J. Phys. Chem. B, 2004, 108, 13427–13437.
55. P. Ren, C. Wu and J. W. Ponder, J. Chem. Theory Comput., 2011, 7, 3143–3161.
56. P. Ren and J. W. Ponder, J. Comput. Chem., 2002, 23, 1497–1506.
57. P. Ren and J. W. Ponder, J. Phys. Chem. B, 2003, 107, 5933–5947.
58. L. Lagardère, L.-H. Jolly, F. Lipparini, F. Aviat, B. Stamm, Z. F. Jing, M. Harger, H. Torabifard, G. A. Cisneros, M. J. Schnieders, N. Gresh, Y. Maday, P. Y. Ren, J. W. Ponder and J.-P. Piquemal, Chem. Sci., 2018, 9, 956–972.
59. A. J. Stone, J. Chem. Theory Comput., 2005, 1, 1128–1132.
60. L. Martínez, R. Andrade, E. G. Birgin and J. M. Martínez, J. Comput. Chem., 2009, 30, 2157–2164.
61. U. Essmann, L. Perera, M. L. Berkowitz, T. Darden, H. Lee and L. G. Pedersen, J. Chem. Phys., 1995, 103, 8577–8593.
62. R. Salomon-Ferrer, A. W. Götz, D. Poole, S. Le Grand and R. C. Walker, J. Chem. Theory Comput., 2013, 9, 3878–3888.
63. Tutorial-4: R.E.D. Server Development/PyRED (q4md-forcefieldtools.org).
64. C. Tian, K. Kasavajhala, K. A. A. Belfon, L. Raguette, H. Huang, A. N. Migues, J. Bickel, Y. Wang, J. Pincay, Q. Wu and C. Simmerling, J. Chem. Theory Comput., 2020, 16, 528–552.
65. M. Zgarbová, J. Šponer, M. Otyepka, T. E. Cheatham, III, R. Galindo-Murillo and P. Jurečka, J. Chem. Theory Comput., 2015, 11, 5723–5736.
66. P. Mark and L. Nilsson, J. Phys. Chem. A, 2001, 105, 9954–9960.
67. C. Schafmeister, W. Ross and V. Romanovski, LEAP, University of California, San Francisco, 1995.
68. D. A. Case, K. Belfon, I. Y. Ben-Shalom, S. R. Brozell, D. S. Cerutti, T. E. Cheatham, III, V. W. D. Cruzeiro, T. A. Darden, R. E. Duke, G. Giambasu, M. K. Gilson, H. Gohlke, A. W. Goetz,R. Harris, S. Izadi, S. A. Izmailov, K. Kasavajhala, A. Kovalenko, R. Krasny, T. Kurtzman, T. S. Lee, S. LeGrand, P. Li, C. Lin, J. Liu, T. Luchko, R. Luo, V. Man, K. M. Merz, Y. Miao, O. Mikhailovskii, G. Monard, H. Nguyen, A. Onufriev, F. Pan, S. Pantano, R. Qi, D. R. Roe, A. Roitberg, C. Sagui, S. Schott-Verdugo, J. Shen, C. L. Simmerling, N. R. Skrynnikov, J. Smith, J. Swails, R. C. Walker, J. Wang, L. Wilson, R. M. Wolf, X. Wu, Y. Xiong, Y. Xue, D. M. York and P. A. Kollman, AMBER 2020, University of California, San Francisco, 2020.
69. R. Zwanzig, J. Stat. Phys., 1973, 9, 215–220.
70. D. T. Gillespie, J. Chem. Phys., 2000, 113, 297–306.
71. P. J. Steinbach and B. R. Brooks, J. Comput. Chem., 1994, 15, 667–683.
72. D. R. Roe and T. E. Cheatham, III, J. Chem. Theory Comput., 2013, 9, 3084–3095.
73. S. v d Walt, S. C. Colbert and G. Varoquaux, Comput. Sci. Eng., 2011, 13, 22–30.
74. J. D. Hunter, Comput. Sci. Eng., 2007, 9, 90–95.
75. W. McKinney in Data Structures for Statistical Computing in Python, Proceedings of the 9th Python in Science Conference, Austin, TX, 2010, vol. 445, pp. 51–56.
76. E. Leddin, Cisneros Research Group and G. Andres Cisneros. ( 2021). Cisneros Research/AMBER-EDA: First Release (v0.1). Zenodo. https://zenodo.org/records/4469902.

---

Footnote

† Electronic supplementary information (ESI) available: Additional details of MD, EDA (PDF). Additional ESI for the initial coordinates and parameters for all of the studied systems (ESI-1.zip) (ZIP). See DOI: https://doi.org/10.1039/d3cp05572h

---