Toward true DNA base-stacking energies: MP2, CCSD(T), and complete basis set calculations

Computational approaches to predicting the impact of novel bases on RNA structure and stability

The use of computational modeling techniques to gain insight into nucleobase interactions has been a challenging endeavor to date. Accurate treatment requires the tackling of many challenges but also holds the promise of great rewards. The development of effective computational approaches to predict the binding affinities of nucleobases and analogues can, for example, streamline the process of developing novel nucleobase modifications, which should facilitate the development of new RNAi-based therapeutics. This brief review focuses on available computational approaches to predicting base pairing affinity in RNA-based contexts such as nucleobase-nucleobase interactions in duplexes and nucleobase-protein interactions. The challenges associated with such modeling along with potential future directions for the field are highlighted.

Higher order structural effects stabilizing the reverse Watson-Crick Guanine-Cytosine base pair in functional RNAs

The G:C reverse Watson-Crick (W:W trans) base pair, also known as Levitt base pair in the context of tRNAs, is a structurally and functionally important base pair that contributes to tertiary interactions joining distant domains in functional RNA molecules and also participates in metabolite binding in riboswitches. We previously indicated that the isolated G:C W:W trans base pair is a rather unstable geometry, and that dicationic metal binding to the Guanine base or posttranscriptional modification of the Guanine can increase its stability. Herein, we extend our survey and report on other H-bonding interactions that can increase the stability of this base pair. To this aim, we performed a bioinformatics search of the PDB to locate all the occurencies of G:C trans base pairs. Interestingly, 66% of the G:C trans base pairs in the PDB are engaged in additional H-bonding interactions with other bases, the RNA backbone or structured water molecules. High level quantum mechanical calculations on a data set of representative crystal structures were performed to shed light on the structural stability and energetics of the various crystallographic motifs. This analysis was extended to the binding of the preQ1 metabolite to a preQ1-II riboswitch.

The intrinsic conformational features of amino acids from a protein coil library and their applications in force field development

The local conformational (φ, ψ, χ) preferences of amino acid residues remain an active research area, which are important for the development of protein force fields. In this perspective article, we first summarize spectroscopic studies of alanine-based short peptides in aqueous solution. While most studies indicate a preference for the P(II) conformation in the unfolded state over α and β conformations, significant variations are also observed. A statistical analysis from various coil libraries of high-resolution protein structures is then summarized, which gives a more coherent view of the local conformational features. The φ, ψ, χ distributions of the 20 amino acids have been obtained from a protein coil library, considering both backbone and side-chain conformational preferences. The intrinsic side-chain χ(1) rotamer preference and χ(1)-dependent Ramachandran plot can be generally understood by combining the interaction of the side-chain Cγ/Oγ atom with two neighboring backbone peptide groups. Current all-atom force fields such as AMBER ff99sb-ILDN, ff03 and OPLS-AA/L do not reproduce these distributions well. A method has been developed by combining the φ, ψ plot of alanine with the influence of side-chain χ(1) rotamers to derive the local conformational features of various amino acids. It has been further applied to improve the OPLS-AA force field. The modified force field (OPLS-AA/C) reproduces experimental (3)J coupling constants for various short peptides quite well. It also better reproduces the temperature-dependence of the helix-coil transition for alanine-based peptides. The new force field can fold a series of peptides and proteins with various secondary structures to their experimental structures. MD simulations of several globular proteins using the improved force field give significantly less deviation (RMSD) to experimental structures. The results indicate that the local conformational features from coil libraries are valuable for the development of balanced protein force fields.

Orbital-based insights into parallel-displaced and twisted conformations in π-π interactions

Dispersion and electrostatics are known to stabilize π-π interactions, but the preference for parallel-displaced (PD) and/or twisted (TW) over sandwiched (S) conformations is not well understood. Orbital interactions are generally believed to play little to no role in π-stacking. However, orbital analysis of the dimers of benzene, pyridine, cytosine and several polyaromatic hydrocarbons demonstrates that PD and/or TW structures convert one or more π-type dimer MOs with out-of-phase or antibonding inter-ring character at the S stack to in-phase or bonding in the PD/TW stack. This change in dimer MO character can be described in terms of a qualitative stack bond order (SBO) defined as the difference between the number of occupied in-phase/bonding and out-of-phase/antibonding inter-ring π-type MOs. The concept of an SBO is introduced here in analogy to the bond order in molecular orbital theory. Thus, whereas the SBO of the S structure is zero, parallel displacement or twisting the stack results in a non-zero SBO and overall bonding character. The shift in bonding/antibonding character found at optimal PD/TW structures maximizes the inter-ring density, as measured by intermolecular Wiberg bond indices (WBIs). Values of WBIs calculated as a function of the parallel-displacement are found to correlate with the dispersion and other contributions to the π-π interaction energy determined by the highly accurate density-fitting DFT symmetry adapted perturbation theory (DF-DFT-SAPT) method. These DF-DFT-SAPT calculations also suggest that the dispersion and other contributions are maximized at the PD conformation rather than the S when conducted on a potential energy curve where the inter-ring distance is optimized at fixed slip distances. From these results of this study, we conclude that descriptions of the qualitative manner in which orbitals interact within π-stacking interactions can supplement high-level calculations of the interaction energy and provide an intuitive tool for applications to crystal design, molecular recognition and other fields where non-covalent interactions are important.

A dynamic structural model of expanded RNA CAG repeats: a refined X-ray structure and computational investigations using molecular dynamics and umbrella sampling simulations

One class of functionally important RNA is repeating transcripts that cause disease through various mechanisms. For example, expanded CAG repeats can cause Huntington's and other disease through translation of toxic proteins. Herein, a crystal structure of r[5'UUGGGC(CAG)3GUCC]2, a model of CAG expanded transcripts, refined to 1.65 Å resolution is disclosed that shows both anti-anti and syn-anti orientations for 1 × 1 nucleotide AA internal loops. Molecular dynamics (MD) simulations using AMBER force field in explicit solvent were run for over 500 ns on the model systems r(5'GCGCAGCGC)2 (MS1) and r(5'CCGCAGCGG)2 (MS2). In these MD simulations, both anti-anti and syn-anti AA base pairs appear to be stable. While anti-anti AA base pairs were dynamic and sampled multiple anti-anti conformations, no syn-anti ↔ anti-anti transformations were observed. Umbrella sampling simulations were run on MS2, and a 2D free energy surface was created to extract transformation pathways. In addition, an explicit solvent MD simulation over 800 ns was run on r[5'GGGC(CAG)3GUCC]2, which closely represents the refined crystal structure. One of the terminal AA base pairs (syn-anti conformation), transformed to anti-anti conformation. The pathway followed in this transformation was the one predicted by umbrella sampling simulations. Further analysis showed a binding pocket near AA base pairs in syn-anti conformations. Computational results combined with the refined crystal structure show that global minimum conformation of 1 × 1 nucleotide AA internal loops in r(CAG) repeats is anti-anti but can adopt syn-anti depending on the environment. These results are important to understand RNA dynamic-function relationships and to develop small molecules that target RNA dynamic ensembles.

Impact of geometry optimization on base-base stacking interaction energies in the canonical A- and B-forms of DNA

Base stacking is known to make an important contribution to the stability of DNA and RNA, and accordingly, significant efforts are ongoing to calculate stacking energies using ab initio quantum mechanical methods. To date, impressive improvements have been made in the model chemistries used to perform stacking energy calculations, including extensions that include robust treatments of electron correlation with extended basis sets, as required to treat interactions where dispersion makes a significant contribution. However, those efforts typically use rigid monomer geometries when calculating the interaction energies. To overcome this, in the present work, we describe a novel internal coordinate definition that allows the relative, intermolecular orientation of stacked base monomers to be constrained during geometry optimizations while allowing full optimization of the intramolecular degrees of freedom. Use of the novel reference frame to calculate the impact of full geometry optimization versus constraining the bases to be planar on base monomer stacking energies, combined with density-fitted, spin-component scaling MP2 treatment of electron correlation, shows that full optimization makes the average stacking energy more favorable by -3.4 and -1.5 kcal/mol for the canonical A and B conformations of the 16 5' to 3' base stacked monomers. Thus, treatment of geometry optimization impacts the stacking energies to an extent similar to or greater than the impact of current state of the art increases in the rigor of the model chemistry itself used to treat base stacking. Results also indicate that stacking favors the B-form of DNA, though the average difference versus the A-form decreases from -2.6 to -0.6 kcal/mol when the intramolecular geometry is allowed to fully relax. However, stacking involving cytosine is shown to favor the A-form of DNA, with that contribution generally larger in the fully optimized bases. The present results show the importance of allowing geometry optimization, as well as properly treating the appropriate model chemistry, in studies of nucleic acid base stacking.

Infrared and electron spin resonance spectral studies of some copper purine and pyrimidine complexes

Copper guanine and barbital complexes were prepared and characterized by elemental analyses and spectral measurements. The data typified the formation of stoichiometries 1:1 (M:L) with possible Cu-Cu interaction "association". The complexes are with different geometries: square planar, square pyramidal and tetrahedral. The mode of bonding was identified by IR spectra. EPR spectra of the powdered complexes were recorded at X band at the room temperature. Different ESR parameters were calculated and discussed: g(//), g(⊥), A(//), [g], G, F, K, α(2). Molecular modeling techniques and quantum chemical methods have been performed for copper complexes to correlate the chemical structures of the complexes with their physical molecular properties. Bond lengths, bond orders, bond angles, dihedral angles, close contact, dipole moment (μ), sum of the total negative charge (STNC), electronegativity (χ), chemical potential (Pi), global hardness (η), softness (σ), the highest occupied molecular orbital energy (E(HOMO)), the lowest unoccupied molecular orbital energy (E(LUMO)) and the energy gap (ΔE) were calculated using PM3 semi-empirical and Molecular Mechanics (MM+) methods. The study displays a good correlation between the theoretical and experimental data which confirms the reliability of the quantum chemical methods.

Quantum-mechanical analysis of the energetic contributions to π stacking in nucleic acids versus rise, twist, and slide

Symmetry-adapted perturbation theory (SAPT) is applied to pairs of hydrogen-bonded nucleobases to obtain the energetic components of base stacking (electrostatic, exchange-repulsion, induction/polarization, and London dispersion interactions) and how they vary as a function of the helical parameters Rise, Twist, and Slide. Computed average values of Rise and Twist agree well with experimental data for B-form DNA from the Nucleic Acids Database, even though the model computations omitted the backbone atoms (suggesting that the backbone in B-form DNA is compatible with having the bases adopt their ideal stacking geometries). London dispersion forces are the most important attractive component in base stacking, followed by electrostatic interactions. At values of Rise typical of those in DNA (3.36 Å), the electrostatic contribution is nearly always attractive, providing further evidence for the importance of charge-penetration effects in π-π interactions (a term neglected in classical force fields). Comparison of the computed stacking energies with those from model complexes made of the "parent" nucleobases purine and 2-pyrimidone indicates that chemical substituents in DNA and RNA account for 20-40% of the base-stacking energy. A lack of correspondence between the SAPT results and experiment for Slide in RNA base-pair steps suggests that the backbone plays a larger role in determining stacking geometries in RNA than in B-form DNA. In comparisons of base-pair steps with thymine versus uracil, the thymine methyl group tends to enhance the strength of the stacking interaction through a combination of dispersion and electrosatic interactions.

Guanine base stacking in G-quadruplex nucleic acids

G-quadruplexes constitute a class of nucleic acid structures defined by stacked guanine tetrads (or G-tetrads) with guanine bases from neighboring tetrads stacking with one another within the G-tetrad core. Individual G-quadruplexes can also stack with one another at their G-tetrad interface leading to higher-order structures as observed in telomeric repeat-containing DNA and RNA. In this study, we investigate how guanine base stacking influences the stability of G-quadruplexes and their stacked higher-order structures. A structural survey of the Protein Data Bank is conducted to characterize experimentally observed guanine base stacking geometries within the core of G-quadruplexes and at the interface between stacked G-quadruplex structures. We couple this survey with a systematic computational examination of stacked G-tetrad energy landscapes using quantum mechanical computations. Energy calculations of stacked G-tetrads reveal large energy differences of up to 12 kcal/mol between experimentally observed geometries at the interface of stacked G-quadruplexes. Energy landscapes are also computed using an AMBER molecular mechanics description of stacking energy and are shown to agree quite well with quantum mechanical calculated landscapes. Molecular dynamics simulations provide a structural explanation for the experimentally observed preference of parallel G-quadruplexes to stack in a 5′–5′ manner based on different accessible tetrad stacking modes at the stacking interfaces of 5′–5′ and 3′–3′ stacked G-quadruplexes.

A reinvestigation of the dimer of para-benzoquinone and pyrimidine with MP2, CCSD(T), and DFT using functionals including those designed to describe dispersion

We reevaluate the interaction of pyridine and p-benzoquinone using functionals designed to treat dispersion. We compare the relative energies of four different structures: stacked, T-shaped (identified for the first time), and two planar H-bonded geometries using these functionals (B97-D, ωB97x-D, M05, M05-2X, M06, M06L, and M06-2X), other functionals (PBE1PBE, B3LYP, X3LYP), MP2, and CCSD(T) using basis sets as large as cc-pVTZ. The functionals designed to treat dispersion behave erratically as the predictions of the most stable structure vary considerably. MP2 predicts the experimentally observed structure (H-bonded) to be the least stable, while single-point CCSD(T) at the MP2 optimized geometry correctly predicts the observed structure to be the most stable. We have confirmed the assignment of the experimental structure using new calculations of the vibrational frequency shifts previously used to identify the structure. The MP2/cc-pVTZ vibrational calculations are in excellent agreement with the observations. All methods used to calculate the energies provide vibrational shifts that agree with the observed structure even though most do not predict this structure to be most stable. The implications for evaluating possible π-stacking in biologically important systems are discussed.

Computational Evaluation of Intermolecular Interactions of a Universal Base 3-Nitropyrrole in Stacked Dimers and DNA Duplexes

The stacking interactions between a universal base of 3-nitropyrrole (3NP) and four canonical nucleobases were studied by means of ab initio molecular orbital calculations. The stabilities of the complexes are comparable to those of the stacked dimers of canonical bases reported previously. The detailed analysis of the interaction energies revealed the importance of the dipole-dipole interaction included in the Hartree-Fock terms to determine the geometry dependence of the stacking energies. It was also clarified that the dispersion energies included in the electron-correlation terms were essential to obtain adequate stabilities. The contribution of the nitro group was evaluated by the comparative studies of pyrrole and 3NP. The increased molecular dipole moment and surface are expected to account for the enhancement of the stability of the stacked dimers containing 3NP. The force field parameters required for calculation of the molecular mechanics of 3NP were obtained for 3NP on the basis of these molecular orbital calculations. The energy-minimized structures obtained by the molecular mechanics calculations of 3NP accorded with those obtained by the molecular orbital calculations described above. A DNA duplex structure containing 3NP-A, 3NP-T, or 3NP-C was calculated by use of these force field parameters. In the case of 3NP-A, the computationally calculated structure was in good agreement with that previously determined by use of (1)H-NMR except for the orientation of the nitro group.

Basis set dependence of higher-order correlation effects in π-type interactions

The basis set dependence of higher-order correlation effects on π-type interaction energies was examined by scanning the potential energy surfaces of five dimer systems. The dimers of acetylene (H-C≡C-H), diacetylene (H-C≡C-C≡C-H), cyanogen (N≡C-C≡N), diphosphorous (P≡P), and 1,4-diphosphabutadiyne (P≡C-C≡P) were studied in three different configurations: cross, parallel-displaced, and t-shaped. More than 800 potential energy curves (PECs) were generated by computing the interaction energies for all 15 dimer configurations over a range of intermolecular distances with the MP2, coupled-cluster single double (CCSD), and coupled-cluster single double triple (CCSD(T)) methods in conjunction with 21 basis sets ranging from a small 6-31G*(0.25) split-valence basis set to a large aug-cc-pVQZ correlation consistent basis set. Standard extrapolation techniques were also used to construct MP2, CCSD, and CCSD(T) complete basis set (CBS) limit PECs as well as CBS limit higher-order correlation corrections based on the differences between CCSD(T) and MP2 interaction energies, denoted δ(MP2)(CCSD(T)), and the corresponding differences between CCSD(T) and CCSD interactions energies, denoted δ(CCSD)(CCSD(T)). Double-ζ basis sets struggled to reproduce the former but provided quite reasonable descriptions of the latter as long as diffuse functions were included. The aug-cc-pVDZ basis deviated from the δ(CCSD)(CCSD(T)) CBS limit by only 0.06 kcal mol(-1) on average and never by more than 0.24 kcal mol(-1), whereas the corresponding deviations were approximately twice that for the δ(MP2)(CCSD(T)) term. While triple-ζ basis sets typically improved results, only aug-cc-pVTZ provided appreciable improvement over utilizing the aug-cc-pVDZ basis set to compute δ(CCSD)(CCSD(T)). Counterpoise (CP) corrections were also applied to all double- and triple-ζ basis sets, but they rarely yielded a better description of these higher-order correlation effects. CP corrections only consistently improved results when the aug-cc-pVDZ basis set was used to compute δ(MP2)(CCSD(T)), yielding mean and maximum absolute deviations from the CBS values of 0.10 and 0.39 kcal mol(-1), respectively, for all five dimer systems.

Magnitude and nature of carbohydrate-aromatic interactions in fucose-phenol and fucose-indole complexes: CCSD(T) level interaction energy calculations

The CH/π contact structures of the fucose-phenol and fucose-indole complexes and the stabilization energies by formation of the complexes (E(form)) were studied by ab initio molecular orbital calculations. The three types of interactions (CH/π and OH/π interactions and OH/O hydrogen bonds) were compared and evaluated in a single molecular system and at the same level of theory. The E(form) calculated for the most stable CH/π contact structure of the fucose-phenol complex at the CCSD(T) level (-4.9 kcal/mol) is close to that for the most stable CH/π contact structure of the fucose-benzene complex (-4.5 kcal/mol). On the other hand the most stable CH/π contact structure of the fucose-indole complex has substantially larger E(form) (-6.5 kcal/mol). The dispersion interaction is the major source of the attraction in the CH/π contact structures of the fucose-phenol and fucose-indole complexes as in the case of the fucose-benzene complex. The electrostatic interactions in the CH/π contact structures are small (less than 1.5 kcal/mol). The nature of the interactions between the nonpolar surface of the carbohydrate and aromatic rings is completely different from that of the conventional hydrogen bonds where the electrostatic interaction is the major source of the attraction. The distributed multipole analysis and DFT-SATP analysis show that the dispersion interactions in the CH/π contact structure of fucose-indole complex are substantially larger than those in the CH/π contact structures of fucose-benzene and fucose-phenol complexes. The large dispersion interactions are responsible for the large E(form) for the fucose-indole complex.

Benchmarking AMBER force fields for RNA: comparisons to NMR spectra for single-stranded r(GACC) are improved by revised χ torsions

Accurately modeling unpaired regions of RNA is important for predicting structure, dynamics, and thermodynamics of folded RNA. Comparisons between NMR data and molecular dynamics simulations provide a test of force fields used for modeling. Here, NMR spectroscopy, including NOESY, (1)H-(31)P HETCOR, DQF-COSY, and TOCSY, was used to determine conformational preferences for single-stranded GACC RNA. The spectra are consistent with a conformational ensemble containing major and minor A-form-like structures. In a series of 50 ns molecular dynamics (MD) simulations with the AMBER99 force field in explicit solvent, initial A-form-like structures rapidly evolve to disordered conformations. A set of 50 ns simulations with revised χ torsions (AMBER99χ force field) gives two primary conformations, consistent with the NMR spectra. A single 1.9 μs MD simulation with the AMBER99χ force field showed that the major and minor conformations are retained for almost 68% of the time in the first 700 ns, with multiple transformations from A-form to non-A-form conformations. For the rest of the simulation, random-coil structures and a stable non-A-form conformation inconsistent with NMR spectra were seen. Evidently, the AMBER99χ force field improves structural predictions for single-stranded GACC RNA compared to the AMBER99 force field, but further force field improvements are needed.