Thermodynamic Insights of Base Flipping in TNA Duplex: Force Fields, Salt Concentrations, and Free-Energy Simulation Methods
Zhaoxi Sun, John Z. H. Zhang
Abstract
Open AccessCCS ChemistryRESEARCH ARTICLE1 Feb 2021Thermodynamic Insights of Base Flipping in TNA Duplex: Force Fields, Salt Concentrations, and Free-Energy Simulation Methods Zhaoxi Sun and John Z.H. Zhang Zhaoxi Sun *Corresponding author: E-mail Address: [email protected] E-mail Address: [email protected] State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062 and John Z.H. Zhang *Corresponding author: E-mail Address: [email protected] E-mail Address: [email protected] State Key Laboratory of Precision Spectroscopy, School of Chemistry and Molecular Engineering, East China Normal University, Shanghai 200062 NYU-ECNU Center for Computational Chemistry at NYU Shanghai, Shanghai 200062 Department of Chemistry, New York University, NY, NY 10003 https://doi.org/10.31635/ccschem.020.202000202 SectionsSupplemental MaterialAboutAbstractPDF ToolsAdd to favoritesDownload CitationsTrack Citations ShareFacebookTwitterLinked InEmail Threofuranosyl nucleic acid (TNA) is an analogue of DNA with a shift in the internucleotide linkages from the wild-type 5'-to-3' direction to 3'-to-2.' This alteration leads to higher chemical stability, less reactive groups, and lower conformational flexibility. Experimental observations indicate that these characteristic changes are attributable to a minimal perturbation of the interaction network, but the thermodynamic stability of the duplex remains unaltered in the TNA mutation. We applied the equilibrium and nonequilibrium free-energy simulations employing three popular assisted model building with energy refinement (AMBER) force fields for nucleotides to investigate this mutation-dependent behavior in the base flipping from T (DNA) residue to the T-to-TFT mutation (TNA) computationally. The force fields were performed similarly, as described in the base-paired state. However, after exploring the high-energy regions with free-energy simulations, we observed that these three force fields behaved differently. Previous reports conclude that the net-neutral and excess-salt simulations provided similar results. Nonetheless, our free-energy simulation indicated that the presence of excess salt affected the thermodynamic stability. The free-energy barrier along the base-flipping pathway was generally elevated upon the addition of excess salts, but the relative height of the free-energy barriers in DNA and TNA duplexes did not change significantly. This phenomenon emphasizes the importance of adding sufficient salts in the simulation scheme to reproduce the experimental condition. Download figure Download PowerPoint Introduction The capacity of the Watson–Crick (WC) base pair is observed in various alternative sugars (e.g., hexopyranoses, pentopyranoses, and tetrofuranoses).1–8 One of the simplest analogues of DNA is threofuranosyl nucleic acid (TNA), in which the internucleotide linkages are shifted from the wild-type 5'-to-3' direction to 3'-to-2' one.9,10 The generated TNA has been shown to have the following benefits: (1) It is able to form stable duplexes with DNA, RNA, and TNA,4,7,11–16 which enables it to transfer information to DNA or RNA. (2) The enzymatic and nonenzymatic polymerization of the TNA nucleotide is possible due to this cross-pairing with other polynucleotides.17–23 (3) TNA is often assumed to be a good predecessor of RNA,24,25 but the chemical stability of TNA is higher, and the number of reactive groups is smaller, which leads to fewer side reactions, and thus, more faithful copying. (4) Regioselectivity is no longer a problem due to the presence of only two hydroxyl groups in each sugar moiety of the TNA polymer.24 (6) A direct outcome of the change of the nucleotide linkage in TNA is the alternation of the number of covalent bonds between consecutive phosphates, that is, from 6 in DNA to 5 in TNA, which reduces the conformational flexibility, possibly by making TNA more suitable for information storage than DNA. Favorably, the thermodynamic stability of the duplex and the interaction network are often unaltered upon the TNA mutation.2 For instance, the stacking interactions are virtually unchanged in TNA-modified Dickerson–Drew dodecamer (DDD). The O4' atoms in the tetrose sugars also share remarkably similar interaction network with those of the deoxyribose O4' atoms in the DDD duplex.26 The unique functional role of nucleotides makes them one of the key targets in understanding biological processes. The genetic code deposited in the molecule hides inside their duplex and triplex structures, and thus, inaccessible to other biomolecules. When the duplex is activated due to the change in the surrounding conditions or the intrinsic fluctuation of the duplex and the base flips outward to be exposed to solvent, enabling the bases to interact with enzymes triggered biological processes such as modifications of the nucleotide sequence.27,28 Computer simulation is a powerful tool to study the dynamics of DNA-related molecules and atomic details.29–36 Technically, molecular dynamics (MD) simulations usually are unable to model base flipping due to a vast difference between the timescale of the base-flipping event and the time step for integrating the equations of motion. Specifically, the free-energy penalty of flipping a base is about 10 kcal/mol.37–42 As a result, the timescale of base flipping is about millisecond (ms).43 By contrast, to ensure numerical stability of integrating the equations of motion, the time step in MD simulations could not exceed several femtoseconds (fs). Such a huge gap between the timescale accessible in MD simulations and that of base flipping makes it hard to acquire converged statistics in computer simulations. Thus, to obtain converged and statistically meaningful data, equilibrium-enhanced sampling techniques such as umbrella sampling44–48 and replica-exchange methods49–55 are often employed. For instance, an umbrella sampling is applied to study the protonation-dependent behavior of base flipping (P-DBBF) in RNA systems, the P-DBBF was found to couple with a syn-to-anti transformation of a guanine (G) group in the G-adenine (A) mismatch.39 Although the nonequilibrium technique of steered MD (SMD)56–63 is less frequently employed to investigate the thermodynamics in base flipping, they do have significant potentials and showed similar performance in the construction of the free-energy landscape of base flipping, compared with equilibrium-enhanced sampling techniques.38,64 Successful applications have been reported in recent years that include, the variation of base-flipping free-energy landscapes upon the sulfur (S) substitution at G groups in G-cytosine (C) base pairs and G-thymine (T) mismatches.64 However, the existing publications often focus on the mutation-dependent behaviors of the nucleotide systems such as synthetic nucleobase pairs and naturally occurring mutations.39,64–73 The behavior of the duplex could also be influenced by the alteration in the nucleotide sugar moiety, which is less frequently studied. Therefore, we selected TNA as a polynucleotide duplex model to investigate the effect of the alteration in the sugar moiety on the thermodynamics of the duplex. The thermodynamic profiles along the base-flipping pathway in wild-type DNA and the mutant TNA duplexes were constructed from equilibrium and nonequilibrium free-energy simulations. Comparisons of the performances of the two different free-energy simulation methods are provided. All-atom descriptions of the DNA systems have been more popular in the computational community in recent years due to the increase in computation power. The mostly used derivatives of assisted model building with energy refinement (AMBER) force fields are bsc174 and OL15,75–77 both of which are based on the earlier AMBER force field named parm9978 with its modification bsc0 counterpart.79 The dihedral parameter is the most challenging part in the development of nucleic acid force fields. Recent established force fields focus mostly on the refinement of this dihedral term or angle torsion potentials. The physical meaning of the torsion potentials is not well-defined; however, pieces of evidence suggest that it is closely related to the contributions from the electronic structures used as the final step in tuning the force field. Recent benchmark studies on force fields showed that these two last-generation force fields performed similarly in describing the structural ensemble of canonical structures of DNA in long brute-force simulations, and the conformational ensemble described by the last-generation AMBER force fields was found to be comparable with nuclear magnetic resonance (NMR)-derived ones.40 However, in our recent work based on flipping in AT tracts, we observed that the free-energy barriers along the base-flipping pathway were significantly different under different force fields.38 Additionally, the base flipping of a mutated DNA duplex also exhibited this behavior,64 but both last-generation AMBER force fields were able to provide accurate quantitative results for the substitution-induced variation of the free-energy barrier along the base-flipping pathway. Accordingly, in our current work, we performed simulations with all of the three AMBER force fields, including the oldest bsc0 modification and the newer OL15 and bsc1 modifications to assess their thermodynamic abilities in nucleic acid systems (i.e., DNA and TNA duplexes). We envisaged that the ion concentration in the solvated system might also influence the thermodynamic and kinetic behaviors of the different systems regarding nucleic acid simulations. Since nucleic acid systems are polyanions with negative charges at the outer phosphate groups, neutralization often occurs by the addition of Na+ cations. Thus, by adding excess salt to achieve physiological ion concentrations of 0.1 M or 0.15 M, we obtained a scheme for excess-salt simulations in which the structural dynamics of the nucleic acid systems could not be altered significantly at these excess-salt concentrations.80–83 Meanwhile, in the absence of excess salts, these observations were noted from unbiased simulations, where the phase space regions explored were close to the canonical structure. Moreover, in our enhanced sampling simulations, some high-free-energy regions were visited, which might make the conclusions inapplicable. Therefore, we aimed at performing simulations of the reaction systems under the two indicated salt concentrations to investigate the excess-salt effect. Experimental Methods Preparation of the reaction system The original 2'-deoxy-T (T) in DNA could be mutated to (l)-alpha-threofuranosyl-T (TFT) in TNA, as previously described.26 The structures of the DDD and its T-to-TFT mutant were obtained from the crystal structure data derived from the pdbid 1N1O.26 The crystal structure and illustrations of the flipping processes are depicted in Figure 1a, with the middle compound structures showing the difference between the T residue in DNA and the TFT residue in the TNA revealed by the mutation of the canonical 5'-3' DNA backbone is mutated to the 3'-2' TNA backbone. The seventh base residue is at the center of the duplex, and thus, its flipping is likely influenced at a minimal level by the terminal fraying effect, which led us to study the flipping of this base residue. We used the high-quality atomic charge scheme, AM1-BCC,84,85 a reliable amber force field program, AMBER14SB,86 and the general AMBER force field (GAFF)87 in comparative MD simulations study of DNA and TNA to parameterize the mutated T residue, which is named as TFT in our current work. We studied the flipping of the seventh residue of the canonical T mutation to TFT by investigating the variations in thermodynamics of the DNA<->TNA mutation. Three AMBER force fields include the OL15 combination of modifications of the AMBER force field,75–77 another last-generation AMBER force field of bsc1,74 and the oldest bsc0 force field.79 As the two last-generation AMBER force fields were built based on the oldest bsc0, from a thermodynamic perspective, our present study represents an assessment of the improvement made by the modifications incorporated into the last-generation force fields. The whole system was solvated with the MD simulations water molecule model, TIP3P.88,89 The simulation box (i.e. the solvated DNA/TNA) was replicated in whole space by periodic boundary conditions. Then we added nonpolarizable spherical counter ions of Na+ for neutralization, as well as the addition of excess salts (Na–Cl ion pairs) and performed simulations to study the variations in the thermodynamic profiles. We noted that a maximum concentration of excess-salt ions of 90 mM was consistent with the conditions of the experimental study.26 Details regarding the DNA<->TNA MD system preparation are included in the Supporting Information. Figure 1 | (a) An illustration of the flipping of the 7th residue of T/TFT in the middle of the TNA/DNA duplex and the comparison between the T residue in DNA and the TFT residue in TNA. (b) The definition of the collective variable (CV) used to describe the flipping event. The center of masses (COMs) of the four groups (heavy atoms only) in circle, that is, the flipping base T7/TFT7, the sugar moiety of the flipping base, the sugar moiety of T8 (the 3' side of the flipping base), and the bases of the base pair T8-A17, were used to define the dihedral. Download figure Download PowerPoint Free-energy simulation We used the equilibrium-enhanced sampling method of umbrella sampling and our nonequilibrium stratification method to bias the sampling.38 The reaction coordinate or the collective variable (CV) had a slow degree of freedom, important for describing the process of interest. In this case, the flipping dihedral, defined by four centers of masses (COMs), was used as the CV. Notably, the pseudo-dihedral was defined by the flipping base T7/TFT7, the sugar moiety of the flipping base, the sugar moiety of T8 (the 3' side of the flipping base), and the bases of the base pair T8-A17. Only heavy atoms were included in this sampling method. The definition is illustrated in Figure 1b. Similar CVs have been applied widely to several DNA and RNA systems.37,39,64 Thus, the data were reweighted with a variational free-energy profile (vFEP)90 to recover the expectations of observables in the original, unbiased ensemble. The convergence check of the umbrella sampling simulations was performed with the block averaging method, a widely used method to check the convergence of the simulation. In this method, the whole production run was divided into several time blocks to monitor the convergence behavior of the free-energy landscape. When similar or identical PMF estimates were obtained from the different time blocks, a convergence was reached. Hence, the nonequilibrated time blocks could be identified and omitted, and the well-equilibrated time blocks were used later for comparisons. The nonequilibrium stratification method38,64 was applied as the time-dependent biasing potential to drive the system from one conformational state to another. The whole pulling process, ranging from 0° to 360 °, was divided into a series of smaller segments. Then bidirectional pulling was initiated from these equilibrated structures to accumulate the microscopic nonequilibrium works. We employed the Crooks' Equation (CE)91 or Bennett Acceptance Ratio (BAR)92 to reweight the statistics. Further, we tested a new convergence criterion for bidirectional reweighting proposed originally in our previous work. In all of our current simulations, the SHAKE93 algorithm, used for integrating the equations of motion degrees of freedom of coordinates, was applied to examine bond-length constraints for bonds involving hydrogen atoms in all molecules to minimize the fluctuation of chemical bonds.94 Langevin dynamics95 with the collision frequency of 4 ps−1 were employed for temperature regulation at 300 K. Isotropic position scaling along with the Berendsen barostat was employed to regulate the pressure. A time step of 1 fs was used to integrate the equations of motion. A cutoff point of 10 Å for nonbonded interactions in the real space was applied, and the long-range electrostatic forces were determined with the particle mesh Ewald (PME) method.96 The MD simulations were performed with the AMBER97 16 suite, and all other analyses were performed with homemade codes. The full computational details are provided in the Supporting Information. Results and Discussion Equilibrium simulations The first thing to check in free-energy simulations is the convergence, thus, we used block averaging to check it. In our previous base-flipping simulations, 1 ns, 2 ns, or 4 ns time blocks were used to estimate the free-energy profiles, which were estimated to be close to the minimum sampling time for each umbrella window attained in the converged sampling results.38,64 To make the convergence as reliable as possible, in the current case, we used 4 ns time blocks, with the results shown in Supporting Information Figures S1 (OL15), S2 (bsc1), and S3 (bsc0). Note that in our previous study, the statistical error obtained from the bootstrap analysis was much smaller than the systematic error.38 Accordingly, the convergence was hindered mainly by bias elimination rather than variance minimization. Thus, in the current study, we did not apply the statistical error in the free-energy profile. We observed that in the first several time blocks, the free-energy profiles obtained were lower due to the gradual equilibration of the system in each umbrella window. However, in the last several time blocks, the free-energy profiles were higher, and fluctuations were observed. Based on these outcomes, we use the overall results obtained from these last time blocks to calculate the potential of mean forces (PMFs). The summary of the time blocks used is provided in Supporting Information Table S2. As different systems require different lengths of equilibration, the total simulation times for different systems described with varying fields of force under different salt concentrations differed, and the resulting statistics are summarized in Supporting Information Table S1. An empirical inference made from our statistical analyses is that in base-flipping simulations, to obtain converged estimates of the free-energy profile from umbrella sampling simulations, 2 µs simulation time was often required for the full performance of the whole system, and simpler systems with the existence of smaller fluctuations, the simulation time could be even shorter (e.g., 1 µs).38,64 Thus, any result obtained from sub-µs simulations should be treated with care. Nonequilibrium simulation techniques The nonequilibrium stratification method, coupled with the statistically optimal bidirectional estimator, is very robust in constructing the thermodynamic profiles in base flipping.38,64 As a result, we checked its applicability in the current base-flipping cases. The OL15 force field was used as the illustrative discussion since the behaviors of the other force fields were similar. The first critical parameter to check is the autocorrelation of the observed mechanical bias, viz, the autocorrelation of the CV, or the flipping dihedral. We acquired statistical inefficiency in the equilibrium ensemble calculated during the initial configurational sampling procedure, as shown in Supporting Information Figure S4a. This statistical inefficiency had a lower bound of 2 ps, derived from the sampling interval of 2 ps. We observed minor fluctuations in the value of the statistical inefficiency, although the typical value obtained was 2 ps. Therefore, we utilized this value in our later calculations of the sampling time of the nonequilibrium free-energy simulation. Previously, we observed the correlation between the convergence of the free-energy simulation and the relative size of the dimensionless standard deviation (SD), as well as the overlap scalar. We found that if the SD is smaller than the overlap scalar, convergence is reached, and the convergence determined by this criterion agreed with the time-invariant behavior of the free-energy profile and the monotonically decreasing behavior of the SD profile during Therefore, the free-energy we had a at the size required for convergence, determined by the Supporting Information Figure the comparison between the dimensionless SD profile and the profile. We that a SD profile with the with such as or the criterion to be and the convergence was we that the convergence required We to make a more reliable of the size required for convergence by the on the free-energy profile in Supporting Information Figure In each new were added to the data and that the simulation the PMF should unchanged or at a minimal level with the addition of new We found that in the first the free-energy barrier was and some fluctuations were observed. the the PMF did not change significantly with and thus, the convergence was with a total of The time of the SD profile also on the convergence behavior of the the simulation the in all configurational should monotonically with reported phenomenon is that the variance or SD is often more than the energy Therefore, we that the time of the SD profile could provide information that could not be obtained from other statistics (e.g., As shown in the time SD profile in Supporting Information Figure it was that the first several fluctuations in the SD and it was the that an important point of all to monotonically with sampling to were reached. Therefore, the time of the SD profile also revealed that were required to the convergence in nonequilibrium Therefore, all of the three statistics indicated that should be the minimum size for convergence to which this value was used in the later for the nonequilibrium The results obtained under the bsc0 force field are summarized in Supporting Information Figures and We noted that the convergence behavior of the simulation under the bsc0 force field was a than that under However, to of the performance of the nonequilibrium method, we estimated the minimum size as Equilibrium nonequilibrium simulations the convergence behavior of the equilibrium and nonequilibrium free-energy simulations, we to check the of the free-energy estimates obtained from these two In Figure 2 and Supporting Information Figure we direct of the free-energy profiles constructed from these two free-energy methods for the wild-type DNA duplex, and the TNA mutant described the OL15 and bsc0 force fields. We an between the free-energy profiles of the two which indicated the and applicability of the nonequilibrium method in constructing the thermodynamic information in the base-flipping process of DNA-related Figure 2 | the OL15 force the between the free-energy profiles of base flipping in (a) the wild-type DNA duplex and (b) the TNA mutant constructed from nonequilibrium stratification and equilibrium umbrella sampling with the variational free-energy profile Download figure Download PowerPoint The however, provided no information about the of the method. Therefore, we calculated the simulation time required for convergence in each method and summarized the statistics in Supporting Information Table We found that the equilibrium umbrella sampling 4 ns umbrella window and a total of to convergence, which an overall simulation time of As for the nonequilibrium method, as indicated we that (i.e., 10 were for convergence to at a required stratification to a total simulation time of indicated that the nonequilibrium method could be much than the equilibrium method in constructing the thermodynamic profiles of base in our previous work, we observed that the nonequilibrium method could be than the equilibrium However, these observations did not indicate that the nonequilibrium method should be as were some in the of the For instance, the umbrella window the stratification was which was the potential of the Further, in the nonequilibrium simulations, the equilibration time in each umbrella window that in each configurational state was included in the likely to also influence the of each method. Thus, a should be that the nonequilibrium techniques could be used to obtain the thermodynamic profiles with similar and as the equilibrium method. variation of the thermodynamic stability of the base pair The difference between DNA and TNA duplexes in the sugar The internucleotide linkages in the wild-type DNA system are to in TNA, is a shift in the to a The number of covalent bonds the consecutive is, thus, from 6 in DNA to 5 in Experimental observations have indicated that DNA and TNA duplexes are of similar thermodynamic Thus, we we could reproduce this behavior in computational In Figure we present between the free-energy profiles along the base-flipping pathway in DNA and TNA described with the three AMBER force fields indicated earlier in the net-neutral and excess-salt simulations. the free-energy profiles shown in Figure we that the of the in both DNA and TNA in different force fields under different salt concentrations were the which indicated that various force fields described the canonical structures in a similar Accordingly, in the brute-force the dynamics in the base-paired state were similar in different force fields. It is that the system only the brute-force method, the phase space explored was close to the canonical structure. Although different force fields have good performance in with long unbiased they might have different behaviors in regions with higher which could be of in the simulation of biological processes. Therefore, we on the height of the free-energy barrier and the of the free-energy profile. Figure | between the free-energy profiles of the wild-type DNA and the TNA mutant described the force field under different ion concentrations (i.e., net-neutral and excess-salt (a)