Toward a full characterization of nucleic acid components in aqueous solution: simulations of nucleosides

Comparison between molecular dynamics potentials for simulation of graphene-based nanomaterials for biomedical applications

OBJECTIVE

This article provides a substantial review of recent research and comparison on molecular dynamics potentials to determine which are most suitable for simulating the phenomena in graphene-based nanomaterials (GBNs).

SIGNIFICANCE

GBNs gain significant attention due to their remarkable properties and potential applications, notably in nanomedicine. However, the physical and chemical characteristics toward macromolecules that justify their nanomedical applications are not yet fully understood. The molecular interaction through molecular dynamic simulation offers the benefits for simulating inorganic molecules like GBNs, with necessary adjustments to account for physical and chemical interactions, or thermodynamic conditions.

METHOD

In this review, we explore various molecular dynamics potentials (force fields) used to simulate interactions and phenomena in graphene-based nanomaterials. Additionally, we offer a brief overview of the benefits and drawbacks of each force fields that available for analysis to assess which one is suitable to study the molecular interaction of graphene-based nanomaterials.

RESULT

We identify and compare various molecular dynamics potentials that available for analyzing GBNs, providing insights into their suitability for simulating specific phenomena in graphene-based nanomaterials. The specification of each force fields and its purpose can be used for further application of molecular dynamics simulation on GBNs.

CONCLUSION

GBNs hold significant promise for applications like nanomedicine, but their physical and chemical properties must be thoroughly studied for safe clinical use. Molecular dynamics simulations, using either reactive or non-reactive MD potentials depending on the expected chemical changes, are essential for accurately modeling these properties, requiring careful selection based on the specific application.

Revealing the excited-state dynamics of cytidine and the role of excited-state proton transfer process

The high photostability of DNAs and RNAs is inextricably related to the photochemical and photophysical properties of their building blocks, nucleobases and nucleosides, which can dissipate the absorbed UV light energy in a harmless manner. The deactivation mechanism of the nucleosides, especially the decay pathways of cytidine (Cyd), has been a matter of intense debate. In the current study, we employ high-level electronic structure calculations combined with excited state non-adiabatic dynamic simulations to provide a clear picture of the excited state deactivation of Cyd in both gas phase and aqueous solution. In both environments, a barrierless decay path driven by the ring-puckering motion and a relaxation channel with a small energy barrier driven by the elongation motion of CO bond are assigned to <200 fs and sub-picosecond decay time component, respectively. The presence of ribose group has a subtle effect on the dynamic behavior of Cyd in gas phase as the ribose-to-base hydrogen/proton transfer process is energetically inaccessible with a sizable energy barrier of about 1.4 eV. However, this energy barrier is significantly reduced in water, especially when an explicit water molecule is present. Therefore, we argue that the long-lived decay channel found in aqueous solution could be assigned to the Cyd-water intermolecular hydrogen/proton transfer process. The present study postulates a novel scenario toward deep understanding the intrinsic photostability of DNAs and RNAs and provides solid evidence to disclose the long history debate of cytidine excited-state decay mechanism, especially for the assignment of experimentally observed time components.

Data-informed reparameterization of modified RNA and the effect of explicit water models: application to pseudouridine and derivatives

Pseudouridine is one of the most abundant post-transcriptional modifications in RNA. We have previously shown that the FF99-derived parameters for pseudouridine and some of its naturally occurring derivatives in the AMBER distribution either alone or in combination with the revised γ torsion parameters (parmbsc0) failed to reproduce their conformational characteristics observed experimentally (Deb et al. in J Chem Inf Model 54:1129–1142, 2014; Deb et al. in J Comput Chem 37:1576–1588, 2016; Dutta et al. in J Chem Inf Model 60:4995–5002, 2020). However, the application of the recommended bsc0 correction did lead to an improvement in the description not only of the distribution in the γ torsional space but also of the sugar pucker distributions. In an earlier study, we examined the transferability of the revised glycosidic torsion parameters (χ_IDRP) for Ψ to its derivatives. We noticed that although these parameters in combination with the AMBER FF99-derived parameters and the revised γ torsional parameters resulted in conformational properties of these residues that were in better agreement with experimental observations, the sugar pucker distributions were still not reproduced accurately. Here we report a new set of partial atomic charges for pseudouridine, 1-methylpseudouridine, 3-methylpseudouridine and 2′-O-methylpseudouridine and a new set of glycosidic torsional parameters (χ_ND) based on chosen glycosidic torsional profiles that most closely corresponded to the NMR data for conformational propensities and studied their effect on the conformational distributions using REMD simulations at the individual nucleoside level. We have also studied the effect of the choice of water model on the conformational characteristics of these modified nucleosides. Our observations suggest that the current revised set of parameters and partial atomic charges describe the sugar pucker distributions for these residues more accurately and that the choice of a suitable water model is important for the accurate description of their conformational properties. We have further validated the revised sets of parameters by studying the effect of substitution of uridine with pseudouridine within single stranded RNA oligonucleotides on their conformational and hydration characteristics.

Quantum Vibrational Spectroscopy of Explicitly Solvated Thymidine in Semiclassical Approximation

In this paper, we demonstrate the possibility to perform spectroscopy simulations of solvated biological species taking into consideration quantum effects and explicit solvation. We achieve this goal by interfacing our recently developed divide-and-conquer approach for semiclassical initial value representation molecular dynamics with the polarizable AMOEBABIO18 force field. The method is applied to the study of solvation of the thymidine nucleoside in two different polar solvents, water and N,N-dimethylformamide. Such systems are made of up to 2476 atoms. Experimental evidence concerning the different behavior of thymidine in the two solvents is well reproduced by our study, even though quantitative estimates are hampered by the limited accuracy of the classical force field employed. Overall, this study shows that semiclassically approximate quantum dynamical studies of explicitly solvated biological systems are both computationally affordable and insightful.

The reorganization energy of compounds upon binding to proteins, from dynamic and solvated bound and unbound states

The intramolecular reorganization energy (ΔE_Reorg) of compounds upon binding to proteins is a component of the binding free energy, which has long received particular attention, for fundamental and practical reasons. Understanding ΔE_Reorg would benefit the science of molecular recognition and drug design. For instance, the tolerable strain energy of compounds upon binding has been elusive. Prior studies found some large ΔE_Reorg values (e.g. > 10 kcal/mol), received with skepticism since they imply excessive opposition to binding. Indeed, estimating ΔE_Reorg is technically difficult. Typically, ΔE_Reorg has been approached by taking two energy-minimized conformers representing the bound and unbound states, and subtracting their conformational energy. This is a drastic oversimplification, liable to conformational collapse of the unbound conformer. Instead, the present work applies extensive molecular dynamics (MD) and the modern OPLS3 force-field to simulate compounds bound and unbound states, in explicit solvent under physically relevant conditions. The thermalized unbound compounds populate multiple conformations, not reducible to one or a few energy-minimized conformers. The intramolecular energies in the bound and unbound states were averaged over pertinent conformational ensembles, and the reorganization enthalpy upon binding (ΔH_Reorg) deduced by subtraction. This was applied to 76 systems, including 43 approved drugs, carefully selected for i) the quality of the bioactive X-ray structures and ii) the diversity of the chemotypes, their properties and protein targets. It yielded comparatively low ΔH_Reorg values (median = 1.4 kcal/mol, mean = 3.0 kcal/mol). A new finding is the observation of negative ΔH_Reorg values. Indeed, reorganization energies do not have to oppose binding, e.g. when intramolecular interactions stabilize preferentially the bound state. Conversely, even with competing water molecules, intramolecular interactions can occur predominantly in the unbound compound, and be replaced by intermolecular counterparts upon protein binding. Such disruption of intramolecular interactions upon binding gives rise to occasional larger ΔH_Reorg values. Such counterintuitive larger ΔH_Reorg values may be rationalized as a redistribution of interactions upon binding, qualitatively compatible with binding.

Molecular Dynamics Simulation of the Conformational Preferences of Pseudouridine Derivatives: Improving the Distribution in the Glycosidic Torsion Space

There are only four derivatives of pseudouridine (Ψ) that are known to occur naturally in RNA as post-transcriptional modifications. We have studied the conformational consequences of pseudouridylation and further modifications using replica exchange molecular dynamics simulations at the nucleoside level, and the simulated conformational preferences were compared with the available experimental (NMR) data. We found that the existing AMBER FF99-derived parameters for these nucleosides did not reproduce the observed experimental features and while the recommended bsc0 correction could be combined with these parameters leading to an improvement in the description of sugar pucker distributions, the χ_OL3 correction could not be applied to these nucleosides as such because of base isomerization. On the other hand, the revised χ torsion parameters (χ_IDRP) for Ψ developed earlier by us (Deb, I., J. Comput. Chem., 2016, 37, 1576-1588) in combination with the AMBER provided parameters and the revised γ torsion parameters generated conformational distributions, which generally were in better agreement with the experimental data. A significant shift of the distribution of base orientation toward the syn conformation was observed with our revised parameter sets compared to the large excess of anti conformation predicted by the FF99 parameters. Overall, our observations indicated that our revised set of parameters (χ_IDRP) for Ψ were also able to generate conformational distributions for all of the derivatives of Ψ in better agreement with the experimental data.

Energy windows for computed compound conformers: covering artefacts or truly large reorganization energies?

The generation of 3D conformers of small molecules underpins most computational drug discovery. Thus, the conformer quality is critical and depends on their energetics. A key parameter is the empirical conformational energy window (ΔEw), since only conformers within ΔEw are retained. However, ΔEw values in use appear unrealistically large. We analyze the factors pertaining to the conformer energetics and ΔEw. We argue that more attention must be focused on the problem of collapsed low-energy conformers. That is due to artificial intramolecular stabilization and occurs even with continuum solvation. Consequently, the conformational energy of extended bioactive structures is artefactually increased, which inflates ΔEw. Thus, this Perspective highlights the issues arising from low-energy conformers and suggests improvements via empirical or physics-based strategies.

Structural insight into antisense gapmer-RNA oligomer duplexes through molecular dynamics simulations

There is an extensive research carrying out on antisense technology and the molecules entering into clinical trials are increasing rapidly. Phosphorothioate (PS) is a chemical modification in which nonbridged oxygen is replaced with a sulfur, consequently providing resistance against nuclease activity. The 2'-4' conformationally restricted nucleoside has the structural features of both 2'-O-methoxy ethyl RNA (MOE), which shows good toxicity profile, and locked nucleic acid (LNA), which shows good binding affinity towards the target RNA. These modifications have been studied and suggested that they can be a potential therapeutic agents in antisense therapy. Mipomersen (ISIS 301012), which contains the novel nucleoside modification has been used to target to apolipoprotein (Apo B), which reduces LDL cholesterol by 6-41%. In this study, classical molecular dynamics (MD) simulations were performed on six different antisense gapmer/target-RNA oligomer duplexes (LNA-PS-LNA/RNA, RcMOE-PS-RcMOE/RNA, ScMOE-PS-ScMOE/RNA, MOE-PS-MOE/RNA, PS-DNA/RNA and DNA/RNA) to investigate the structural dynamics, stability and solvation properties. The LNA, MOE nucleotides present in respective duplexes are showing the structure of A-form and the PS-DNA nucleotides resemble the structure of B-form helix with respect to some of the helical parameters. Free energy calculations suggest that the oligomer, which contains LNA binds to the RNA strongly than other modifications as shown in experimental results. The MOE modified nucleotide, which although had a lower binding affinity but higher solvent accessible surface area (SASA) compared to the other modifications, may be influencing the toxicity and hence may be used it in Mipomersen, the second antisense molecule which is approved by FDA. Communicated by Ramaswamy H. Sarma.

Detailed Analysis of 17β-Estradiol-Aptamer Interactions: A Molecular Dynamics Simulation Study

Micro-pollutants such as 17β-Estradiol (E2) have been detected in different water resources and their negative effects on the environment and organisms have been observed. Aptamers are established as a possible detection tool, but the underlying ligand binding is largely unexplored. In this study, a previously described 35-mer E2-specific aptamer was used to analyse the binding characteristics between E2 and the aptamer with a MD simulation in an aqueous medium. Because there is no 3D structure information available for this aptamer, it was modeled using coarse-grained modeling method. The E2 ligand was positioned inside a potential binding area of the predicted aptamer structure, the complex was used for an 25 ns MD simulation, and the interactions were examined for each time step. We identified E2-specific bases within the interior loop of the aptamer and also demonstrated the influence of frequently underestimated water-mediated hydrogen bonds. The study contributes to the understanding of the behavior of ligands binding with aptamer structure in an aqueous solution. The developed workflow allows generating and examining further appealing ligand-aptamer complexes.

Modeling p K Shift in DNA Triplexes Containing Locked Nucleic Acids

The protonation states for nucleic acid bases are difficult to assess experimentally. In the context of DNA triplex, the protonation state of cytidine in the third strand is particularly important, because it needs to be protonated in order to form Hoogsteen hydrogen bonds. A sugar modification, locked nucleic acid (LNA), is widely used in triplex forming oligonucleotides to target sites in the human genome. In this study, the parameters for LNA are developed in line with the CHARMM nucleic acid force field and validated toward the available structural experimental data. In conjunction, two computational methods were used to calculate the protonation state of the third strand cytidine in various DNA triplex environments: λ-dynamics and multiple pH regime. Both approaches predict p K of this cytidine shifted above physiological pH when cytidine is in the third strand in a triplex environment. Both methods show an upshift due to cytidine methylation, and a small downshift when the sugar configuration is locked. The predicted p K values for cytidine in DNA triplex environment can inform the design of better-binding oligonucleotides.

A combined variable temperature 600 MHz NMR/MD study of the calcium release agent cyclic adenosine diphosphate ribose (cADPR): Structure, conformational analysis, and thermodynamics of the conformational equilibria

A combined variable temperature 600 MHz NMR/molecular dynamics study of the Ca²⁺-release agent cyclic adenosine 5'-diphosphate ribose (cADPR) was conducted. In addition to elucidating the major and minor orientations of the conformationally flexible furanose rings, γ- (C4'-C5'), and β- (C5'-O5') bonds, the thermodynamics (ΔH^o, ΔS^o) associated with each of these conformational equilibria were determined. Both furanose rings were biased towards a south conformation (64-74%) and both β-bonds heavily favored trans conformations. The R-ring γ-bond was found to exist almost exclusively as the γ⁺ conformer, whereas the A-ring γ-bond was a mixture of the γ⁺ and γ^t conformers, with the trans conformer being slightly favored. Enthalpic factors accounted for most of the observed conformational preferences, although the R-ring furanose exists as its major conformation based solely on entropic factors. There was excellent agreement between the NMR and MD results, particularly with regard to the conformer identities, but the MD showed a bias towards γ⁺ conformers. The MD results showed that both N-glycosidic χ-bonds are exclusively syn. Collectively the data allowed for the construction of a model for cADPR in which many of the conformationally flexible units in fact effectively adopt single orientations and where most of the conformational diversity resides in its A-ring furanose and γ-bond.

The free energy of locking a ring: Changing a deoxyribonucleoside to a locked nucleic acid

Locked nucleic acid (LNA), a modified nucleoside which contains a bridging group across the ribose ring, improves the stability of DNA/RNA duplexes significantly, and therefore is of interest in biotechnology and gene therapy applications. In this study, we investigate the free energy change between LNA and DNA nucleosides. The transformation requires the breaking of the bridging group across the ribose ring, a problematic transformation in free energy calculations. To address this, we have developed a 3-step (easy to implement) and a 1-step protocol (more efficient, but more complicated to setup), for single and dual topologies in classical molecular dynamics simulations, using the Bennett Acceptance Ratio method to calculate the free energy. We validate the approach on the solvation free energy difference for the nucleosides thymidine, cytosine, and 5-methyl-cytosine. © 2017 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

Rapid communication capturing the destabilizing effect of dihydrouridine through molecular simulations

The structural effects of the commonly occurring modified nucleoside dihydrouridine (D) observed experimentally in model oligonucleotides include a strong destabilization of the C3'-endo sugar conformation of D, the disruption of stacking interactions of neighboring residues with D and a possible destabilization of the C3'-endo sugar pucker of the 5'-neighboring nucleoside. Our simulations with a combination of a set of parameters for modified RNA residues with the recently developed AMBER FF99χ force field having reoptimized glycosidic torsion angle parameters for standard nucleosides was found to reproduce the destabilizing effect of dihydrouridine better than with the AMBER FF99 force field for nucleic acids for which the parameters for the modified residues were originally developed.

Towards understanding the unbound state of drug compounds: Implications for the intramolecular reorganization energy upon binding

There has been an explosion of structural information for pharmaceutical compounds bound to biological targets, but the conformations and dynamics of compounds free in solution are poorly characterized, if at all. Yet, knowledge of the unbound state is essential to understand the fundamentals of molecular recognition, including the much debated conformational intramolecular reorganization energy of a compound upon binding (ΔEReorg). Also, dependable observation of the unbound compounds is important for ligand-based drug discovery, e.g. with pharmacophore modelling. Here, these questions are addressed with long (⩾0.5μs) state-of-the-art molecular dynamics (MD) simulations of 26 compounds (including 7 approved drugs) unbound in explicit solvent. These compounds were selected to be chemically diverse, with a range of flexibility, and good quality bioactive X-ray structures. The MD-simulated free compounds are compared to their bioactive structure and conformers generated with ad hoc sampling in vacuo or with implicit generalized Born (GB) aqueous solvation models. The GB conformational models clearly depart from those obtained in explicit solvent, and suffer from conformational collapse almost as severe as in vacuo. Thus, the global energy minima in vacuo or with GB are not suitable representations of the unbound state, which can instead be extensively sampled by MD simulations. Many, but not all, MD-simulated compounds displayed some structural similarity to their bioactive structure, supporting the notion of conformational pre-organization for binding. The ligand-protein complexes were also simulated in explicit solvent, to estimate ΔEReorg as an enthalpic difference ΔHReorg between the intramolecular energies in the bound and unbound states. This fresh approach yielded ΔHReorg values⩽6kcal/mol for 18 out of 26 compounds. For three particularly polar compounds 15⩽ΔHReorg⩽20kcal/mol, supporting the notion that ΔHReorg can be substantial. Those large ΔHReorg values correspond to a redistribution of electrostatic interactions upon binding. Overall, the study illustrates how MD simulations offer a promising avenue to characterize the unbound state of medicinal compounds.

Structural Flexibility and Conformation Features of Cyclic Dinucleotides in Aqueous Solutions

Cyclic dinucleotides are able to trigger the innate immune system by activating STING. It was found that the binding affinity of asymmetric 2'3'-cGAMP to symmetric dimer of STING is 3 orders of magnitude higher than that of the symmetric 3'3'-cyclic dinucleotides. Such a phenomenon has not been understood yet. Here we show that the subtle changes in phosphodiester linkage of CDNs lead to their distinct structural properties which correspond to the varied binding affinities. 2'-5' and/or 3'-5' linked CDNs adopt specific while different types of ribose puckers and backbone conformations. That ribose conformations and base types have different propensities for anti or syn glycosidic conformations further affects the overall flexibility of CDNs. The counterbalance between backbone ring tension and electrostatic repulsion, both affected by the ring size, also contributes to the different flexibility of CDNs. Our calculations reveal that the free energy cost for 2'3'-cGAMP to adopt the STING-bound structure is smaller than that for 3'3'-cGAMP and cyclic-di-GMP. These findings may serve as a reference for design of CDN-analogues as vaccine adjuvants. Moreover, the cyclization pattern of CDNs closely related to their physiological roles suggests the importance of understanding structural properties in the study of protein-ligand interactions.

Additive CHARMM force field for naturally occurring modified ribonucleotides

More than 100 naturally occurring modified nucleotides have been found in RNA molecules, in particular in tRNAs. We have determined molecular mechanics force field parameters compatible with the CHARMM36 all‐atom additive force field for all these modifications using the CHARMM force field parametrization strategy. Emphasis was placed on fine tuning of the partial atomic charges and torsion angle parameters. Quantum mechanics calculations on model compounds provided the initial set of target data, and extensive molecular dynamics simulations of nucleotides and oligonucleotides in aqueous solutions were used for further refinement against experimental data. The presented parameters will allow for computational studies of a wide range of RNAs containing modified nucleotides, including the ribosome and transfer RNAs. © 2016 The Authors. Journal of Computational Chemistry Published by Wiley Periodicals, Inc.

A Two-Metal-Ion-Mediated Conformational Switching Pathway for HDV Ribozyme Activation

RNA enzymes serve as a potentially powerful platform from which to design catalysts and engineer new biotechnology. A fundamental understanding of these systems provides insight to guide design. The hepatitis delta virus ribozyme (HDVr) is a small, self-cleaving RNA motif widely distributed in nature, that has served as a paradigm for understanding basic principles of RNA catalysis. Nevertheless, questions remain regarding the precise roles of divalent metal ions and key nucleotides in catalysis. In an effort to establish a reaction mechanism model consistent with available experimental data, we utilize molecular dynamics simulations to explore different conformations and metal ion binding modes along the HDVr reaction path. Building upon recent crystallographic data, our results provide a dynamic model of the HDVr reaction mechanism involving a conformational switch between multiple non-canonical G25:U20 base pair conformations in the active site. These local nucleobase dynamics play an important role in catalysis by modulating the metal binding environments of two Mg²⁺ ions that support catalysis at different steps of the reaction pathway. The first ion plays a structural role by inducing a base pair flip necessary to obtain the catalytic fold in which C75 moves towards to the scissile phosphate in the active site. Ejection of this ion then permits a second ion to bind elsewhere in the active site and facilitate nucleophile activation. The simulations collectively describe a mechanistic scenario that is consistent with currently available experimental data from crystallography, phosphorothioate substitutions, and chemical probing studies. Avenues for further experimental verification are suggested.

Excited-state deactivation in 8-oxo-deoxyguanosine: comparison between anionic and neutral forms

Ab initio explorations of excited-state potential-energy surfaces show that a radiationless deactivation mechanism via intramolecular excited-state proton transfer is available in neutral 8-oxo-deoxyguanosine, whereas it is not available in the anionic form.

Excited-state hydrogen atom abstraction initiates the photochemistry of β-2'-deoxycytidine

Understanding the effects of ultraviolet radiation on nucleotides in solution is an important step towards a comprehensive description of the photochemistry of nucleic acids and their constituents. Apart from having implications for mutagenesis and DNA photoprotection mechanisms, the photochemistry of cytidines is a central element in UV-assisted syntheses of pyrimidine nucleotides under prebiotically plausible conditions. In this contribution, we present UV-irradiation experiments of β-2'-deoxycytidine in aqueous solution involving H-D exchange followed by NMR spectroscopic analysis of the photoproducts. We further elucidate the outcome of these experiments by means of high-level quantum chemical calculations. In particular, we show that prolonged UV-irradiation of cytidine may lead to H-C1' hydrogen atom abstraction by the carbonyl oxygen atom of cytosine. This process may enable photoanomerisation and nucleobase loss, two previously unexplained photoreactions observed in pyrimidine nucleotides.

Interstrand cross-link and bioconjugate formation in RNA from a modified nucleotide

RNA oligonucleotides containing a phenyl selenide derivative of 5-methyluridine were chemically synthesized by solid-phase synthesis. The phenyl selenide is rapidly converted to an electrophilic, allylic phenyl seleneate under mild oxidative conditions. The phenyl seleneate yields interstrand cross-links when part of a duplex and is useful for synthesizing oligonucleotide conjugates. Formation of the latter is illustrated by reaction of an oligonucleotide containing the phenyl selenide with amino acids in the presence of mild oxidant. The products formed are analogous to those observed in tRNA that are believed to be formed posttranslationally via a biosynthetic intermediate that is chemically homologous to the phenyl seleneate.

Balancing the interactions of ions, water, and DNA in the Drude polarizable force field

Recently we presented a first-generation all-atom Drude polarizable force field for DNA based on the classical Drude oscillator model, focusing on optimization of key dihedral angles followed by extensive validation of the force field parameters. Presently, we describe the procedure for balancing the electrostatic interactions between ions, water, and DNA as required for development of the Drude force field for DNA. The proper balance of these interactions is shown to impact DNA stability and subtler conformational properties, including the conformational equilibrium between the BI and BII states, and the A and B forms of DNA. The parametrization efforts were simultaneously guided by gas-phase quantum mechanics (QM) data on small model compounds and condensed-phase experimental data on the hydration and osmotic properties of biologically relevant ions and their solutions, as well as theoretical predictions for ionic distribution around DNA oligomer. In addition, fine-tuning of the internal base parameters was performed to obtain the final DNA model. Notably, the Drude model is shown to more accurately reproduce counterion condensation theory predictions of DNA charge neutralization by the condensed ions as compared to the CHARMM36 additive DNA force field, indicating an improved physical description of the forces dictating the ionic solvation of DNA due to the explicit treatment of electronic polarizability. In combination with the polarizable DNA force field, the availability of Drude polarizable parameters for proteins, lipids, and carbohydrates will allow for simulation studies of heterogeneous biological systems.

Conformational preferences of modified uridines: comparison of AMBER derived force fields

The widespread occurrence of modified residues in RNA sequences necessitates development of accurate parameters for these modifications for reliable modeling of RNA structure and dynamics. A comprehensive set of parameters for the 107 naturally occurring RNA modifications was proposed by Aduri et al. (J. Chem. Theory Comput. 2007, 3, 1464-1475) for the AMBER FF99 force field. In this work, we tested these parameters on a set of modified uridine residues, namely, dihydrouridine, 2-thiouridine, 4-thiouridine, pseudouridine, and uridine-5-oxyacetic acid, by performing molecular dynamics and replica exchange molecular dynamics simulations of these nucleosides. Although our simulations using the FF99 force field did not, in general, reproduce the experimentally observed conformational characteristics well, combination of the parameter set with recent revisions of the FF99 force field for RNA showed noticeable improvement for some of the nucleosides.

The free energy landscape of pseudorotation in 3'-5' and 2'-5' linked nucleic acids

The five-membered furanose ring is a central component of the chemical structure of biological nucleic acids. The conformations of the furanose ring can be analytically described using the concept of pseudorotation, and for RNA and DNA they are dominated by the C2'-endo and C3'-endo conformers. While the free energy difference between these two conformers can be inferred from NMR measurements, a free energy landscape of the complete pseudorotation cycle of nucleic acids in solution has remained elusive. Here, we describe a new free energy calculation method for molecular dynamics (MD) simulations using the two pseudorotation parameters directly as the collective variables. To validate our approach, we calculated the free energy surface of ribose pseudorotation in guanosine and 2'-deoxyguanosine. The calculated free energy landscape reveals not only the relative stability of the different pseudorotation conformers, but also the main transition path between the stable conformations. Applying this method to a standard A-form RNA duplex uncovered the expected minimum at the C3'-endo state. However, at a 2'-5' linkage, the minimum shifts to the C2'-endo conformation. The free energy of the C3'-endo conformation is 3 kcal/mol higher due to a weaker hydrogen bond and a reduced base stacking interaction. Unrestrained MD simulations suggest that the conversion from C3'-endo to C2'-endo and vice versa is on the nanosecond and microsecond time scale, respectively. These calculations suggest that 2'-5' linkages may enable folded RNAs to sample a wider spectrum of their pseudorotation conformations.

Simulating DNA by molecular dynamics: aims, methods, and validation

The structure and dynamics of the B-DNA double helix involves subtle sequence-dependent effects which are decisive for its function, but difficult to characterize. These structural and dynamic effects can be addressed by simulations of DNA sequences in explicit solvent. Here, we present and discuss the state-of-art of B-DNA molecular dynamics simulations with the major force fields in use today. We explain why a critical analysis of the MD trajectories is required to assess their reliability, and estimate the value and limitations of these models. Overall, simulations of DNA bear great promise towards deciphering the structural and physical subtleties of this biopolymer, where much remains to be understood.

The DNA and RNA sugar-phosphate backbone emerges as the key player. An overview of quantum-chemical, structural biology and simulation studies

Knowledge of geometrical and physico-chemical properties of the sugar-phosphate backbone substantially contributes to the comprehension of the structural dynamics, function and evolution of nucleic acids. We provide a side by side overview of structural biology/bioinformatics, quantum chemical and molecular mechanical/simulation studies of the nucleic acids backbone. We highlight main features, advantages and limitations of these techniques, with a special emphasis given to their synergy. The present status of the research is then illustrated by selected examples which include classification of DNA and RNA backbone families, benchmark structure-energy quantum chemical calculations, parameterization of the dihedral space of simulation force fields, incorporation of arsenate into DNA, sugar-phosphate backbone self-cleavage in small RNA enzymes, and intricate geometries of the backbone in recurrent RNA building blocks. Although not apparent from the current literature showing limited overlaps between the QM, simulation and bioinformatics studies of the nucleic acids backbone, there in fact should be a major cooperative interaction between these three approaches in studies of the sugar-phosphate backbone.

Optimization of the CHARMM additive force field for DNA: Improved treatment of the BI/BII conformational equilibrium

The B-form of DNA can populate two different backbone conformations: BI and BII, defined by the difference between the torsion angles ε and ζ (BI = ε-ζ < 0 and BII = ε-ζ > 0). BI is the most populated state, but the population of the BII state, which is sequence dependent, is significant and accumulating evidence shows that BII affects the overall structure of DNA, and thus influences protein-DNA recognition. This work presents a reparametrization of the CHARMM27 additive nucleic acid force field to increase the sampling of the BII form in MD simulations of DNA. In addition, minor modifications of sugar puckering were introduced to facilitate sampling of the A form of DNA under the appropriate environmental conditions. Parameter optimization was guided by quantum mechanical data on model compounds, followed by calculations on several DNA duplexes in the condensed phase. The selected optimized parameters were then validated against a number of DNA duplexes, with the most extensive tests performed on the EcoRI dodecamer, including comparative calculations using the Amber Parm99bsc0 force field. The new CHARMM model better reproduces experimentally observed sampling of the BII conformation, including sampling as a function of sequence. In addition, the model reproduces the A form of the 1ZF1 duplex in 75 % ethanol, and yields a stable Z-DNA conformation of duplex (GTACGTAC) in its crystal environment. The resulting model, in combination with a recent reoptimization of the CHARMM27 force field for RNA, will be referred to as CHARMM36.

What is the conformation of physiologically-active dinucleoside polyphosphates in solution? Conformational analysis of free dinucleoside polyphosphates by NMR and molecular dynamics simulations

Dinucleoside polyphosphates, or dinucleotides (Np(n)N'; N, N' = A, U, G, C; n = 2-7), are naturally occurring ubiquitous physiologically active compounds. Despite the interest in dinucleotides, and the relevance of their conformation to their biological function, the conformation of dinucleotides has been insufficiently studied. Therefore, here we performed conformational analysis of a series of Np(n)N' Na(+) salts (N = A, G, U, C; N' = A, G, U, C; n = 2-5) by various NMR techniques. All studied dinucleotides, except for Up(4/5)U, formed intramolecular base stacking interactions in aqueous solutions as indicated by NMR. The conformation around the glycosidic angle in Np(n)N's was found to be anti/high anti and the preferred conformation around the C4'-C5', C5'-O5' bonds was found to be gauche-gauche (gg). The ribose moiety in Np(n)N's showed a small preference for the S conformation, but when attached to cytosine the ribose ring preferred the N conformation. However, no predominant conformation was observed for the ribose moiety in any of the dinucleotides. Molecular dynamics simulations of Ap(2)A and Ap(4)A Na(+) salts supported the experimental results. In addition, three modes of base-stacking were found for Ap(2/4)A: α-α, β-β and α-β, which exist in equilibrium, while none is dominant. We conclude that natural, free Np(n)N's (n = 2-5) at physiological pH exist mostly in a folded (stacked), rather than extended conformation, in several interconverting stacking modes. Intramolecular base stacking of Np(n)N's does not alter the conformation of each of the nucleotide moieties, which remains the same as that of the mononucleotides in solution.

Intrinsic flexibility of B-DNA: the experimental TRX scale

B-DNA flexibility, crucial for DNA–protein recognition, is sequence dependent. Free DNA in solution would in principle be the best reference state to uncover the relation between base sequences and their intrinsic flexibility; however, this has long been hampered by a lack of suitable experimental data. We investigated this relationship by compiling and analyzing a large dataset of NMR ³¹P chemical shifts in solution. These measurements reflect the BI ↔ BII equilibrium in DNA, intimately correlated to helicoidal descriptors of the curvature, winding and groove dimensions. Comparing the ten complementary DNA dinucleotide steps indicates that some steps are much more flexible than others. This malleability is primarily controlled at the dinucleotide level, modulated by the tetranucleotide environment. Our analyses provide an experimental scale called TRX that quantifies the intrinsic flexibility of the ten dinucleotide steps in terms of Twist, Roll, and X-disp (base pair displacement). Applying the TRX scale to DNA sequences optimized for nucleosome formation reveals a 10 base-pair periodic alternation of stiff and flexible regions. Thus, DNA flexibility captured by the TRX scale is relevant to nucleosome formation, suggesting that this scale may be of general interest to better understand protein-DNA recognition.

CHARMM: the biomolecular simulation program

CHARMM (Chemistry at HARvard Molecular Mechanics) is a highly versatile and widely used molecular simulation program. It has been developed over the last three decades with a primary focus on molecules of biological interest, including proteins, peptides, lipids, nucleic acids, carbohydrates, and small molecule ligands, as they occur in solution, crystals, and membrane environments. For the study of such systems, the program provides a large suite of computational tools that include numerous conformational and path sampling methods, free energy estimators, molecular minimization, dynamics, and analysis techniques, and model-building capabilities. The CHARMM program is applicable to problems involving a much broader class of many-particle systems. Calculations with CHARMM can be performed using a number of different energy functions and models, from mixed quantum mechanical-molecular mechanical force fields, to all-atom classical potential energy functions with explicit solvent and various boundary conditions, to implicit solvent and membrane models. The program has been ported to numerous platforms in both serial and parallel architectures. This article provides an overview of the program as it exists today with an emphasis on developments since the publication of the original CHARMM article in 1983.

Conformational sampling of druglike molecules with MOE and catalyst: implications for pharmacophore modeling and virtual screening

Computational conformational sampling is integral to small molecule pharmaceutical research, for detailed conformational analysis and high-throughput 3D library enumeration. These two regimes were tested in details for the general-purpose modeling program MOE, using its three conformational sampling methods, i.e. systematic search, stochastic search, and Conformation Import. The tests include i) identification of the global energy minimum, ii) reproduction of the bioactive conformation, iii) measures of conformational coverage with 3D descriptors, and iv) compute times. The bioactive conformers are from a new set of 256 diverse, druglike, protein-bound ligands compiled and analyzed with particular care. The MOE results are compared to those obtained from the established program Catalyst. Key parameters controlling the conformational coverage were varied systematically. Coverage and diversity of the conformational space were characterized with unique pharmacophore triplets or quadruplets. Overall, the protocols in both MOE and Catalyst performed well for their intended tasks. MOE performed at least as well as Catalyst for high-throughput library generation and detailed conformational modeling. This work provides a guide and specific recommendations regarding the usage of conformational sampling tools in MOE.

Cytosine ribose flexibility in DNA: a combined NMR 13C spin relaxation and molecular dynamics simulation study

Using (13)C spin relaxation NMR in combination with molecular dynamic (MD) simulations, we characterized internal motions within double-stranded DNA on the pico- to nano-second time scale. We found that the C-H vectors in all cytosine ribose moieties within the Dickerson-Drew dodecamer (5'-CGCGAATTCGCG-3') are subject to high amplitude motions, while the other nucleotides are essentially rigid. MD simulations showed that repuckering is a likely motional model for the cytosine ribose moiety. Repuckering occurs with a time constant of around 100 ps. Knowledge of DNA dynamics will contribute to our understanding of the recognition specificity of DNA-binding proteins such as cytosine methyltransferase.

Collective headgroup conformational transition in twisted micellar superstructures

Predictions on amphiphilic self-assemblies traditionally rely on considerations on molecular shape and charge of the surfactant. In the case of functional surfactants a more sophisticated toolbox becomes necessary to design amphiphiles encoding chemical functionalities that provide additional responsive properties to their self-assemblies. Here we report on a comprehensive and combined structural-spectroscopic characterization of 1,2-dilauroyl-phosphatidyl-adenosine (DLPA) micelles in phosphate buffer. The temperature dependence, more precisely the thermal history of the sample, is explicitly taken into account. The experimental data, supplemented with MD simulations, indicate the presence of two possible states at room temperature, characterized by distinctly different structural properties that depend on the thermal history of the sample. The twisted superstructures, produced by aging DLPA micelles through intermicellar assembly of locally cylindrical aggregates at room temperature, collapse upon warming at 35 °C, yielding aligned filaments and/or wormlike structures. The initial superstructures cannot be recovered by thermal inversion. The reason for this behaviour is that the thermal activation causes a redistribution of syn-anti conformations of adenosine headgroups, as indicated by spectroscopic results (NMR, CD, FTIR), which is then collectively frozen thanks to molecular constraints present in the aggregate.

Insights into structure, dynamics and hydration of locked nucleic acid (LNA) strand-based duplexes from molecular dynamics simulations

Locked nucleic acid (LNA) is a chemically modified nucleic acid with its sugar ring locked in an RNA-like (C3′-endo) conformation. LNAs show extraordinary thermal stabilities when hybridized with DNA, RNA or LNA itself. We performed molecular dynamics simulations on five isosequential duplexes (LNA–DNA, LNA–LNA, LNA–RNA, RNA–DNA and RNA–RNA) in order to characterize their structure, dynamics and hydration. Structurally, the LNA–DNA and LNA–RNA duplexes are found to be similar to regular RNA–DNA and RNA–RNA duplexes, whereas the LNA–LNA duplex is found to have its helix partly unwound and does not resemble RNA–RNA duplex in a number of properties. Duplexes with an LNA strand have on average longer interstrand phosphate distances compared to RNA–DNA and RNA–RNA duplexes. Furthermore, intrastrand phosphate distances in LNA strands are found to be shorter than in DNA and slightly shorter than in RNA. In case of induced sugar puckering, LNA is found to tune the sugar puckers in partner DNA strand toward C3′-endo conformations more efficiently than RNA. The LNA–LNA duplex has lesser backbone flexibility compared to the RNA–RNA duplex. Finally, LNA is less hydrated compared to DNA or RNA but is found to have a well-organized water structure.

Coupling of the guanosine glycosidic bond conformation and the ribonucleotide cleavage reaction: implications for barnase catalysis

To examine the possible relationship of guanine-dependent GpA conformations with ribonucleotide cleavage, two potential of mean force (PMF) calculations were performed in aqueous solution. In the first calculation, the guanosine glycosidic (Gchi) angle was used as the reaction coordinate, and computations were performed on two GpA ionic species: protonated (neutral) or deprotonated (negatively charged) guanosine ribose O2 '. Similar energetic profiles featuring two minima corresponding to the anti and syn Gchi regions were obtained for both ionic forms. For both forms the anti conformation was more stable than the syn, and barriers of approximately 4 kcal/mol were obtained for the anti --> syn transition. Structural analysis showed a remarkable sensitivity of the phosphate moiety to the conformation of the Gchi angle, suggesting a possible connection between this conformation and the mechanism of ribonucleotide cleavage. This hypothesis was confirmed by the second PMF calculations, for which the O2 '--P distance for the deprotonated GpA was used as reaction coordinate. The computations were performed from two selected starting points: the anti and syn minima determined in the first PMF study of the deprotonated guanosine ribose O2'. The simulations revealed that the O2 ' attack along the syn Gchi was more favorable than that along the anti Gchi: energetically, significantly lower barriers were obtained in the syn than in the anti conformation for the O--P bond formation; structurally, a lesser O2 '--P initial distance, and a better suited orientation for an in-line attack was observed in the syn relative to the anti conformation. These results are consistent with the catalytically competent conformation of barnase-ribonucleotide complex, which requires a guanine syn conformation of the substrate to enable abstraction of the ribose H2 ' proton by the general base Glu73, thereby suggesting a coupling between the reactive substrate conformation and enzyme structure and mechanism.

The structure of LNA:DNA hybrids from molecular dynamics simulations: the effect of locked nucleotides

Locked nucleic acids (LNAs) exhibit a modified sugar fragment that is restrained to the C3'-endo conformation. LNA-containing duplexes are rather stable and have a more rigid structure than DNA duplexes, with a propensity for A-conformation of the double helix. To gain detailed insight into the local structure of LNA-modified DNA oligomers (as a foundation for subsequent exploration of the electron-transfer capabilities of such modified duplexes), we carried out molecular dynamics simulations on a set of LNA:DNA 9-mer duplexes and analyzed the resulting structures in terms of base step parameters and the conformations of the sugar residues. The perturbation introduced by a single locked nucleotide was found to be fairly localized, extending mostly to the first neighboring base pairs; such duplexes featured a B-type helix. With increasing degree of LNA modification the structure gradually changed; the duplex with one complete LNA strand assumed a typical A-DNA structure. The relative populations of the sugar conformations agreed very well with NMR data, lending credibility to the validity of the computational protocol.

How many conformers determine the thymidine low-temperature matrix infrared spectrum? DFT and MP2 quantum chemical study

A comprehensive conformational analysis of isolated 2'-beta-deoxy-thymidine (T), canonical DNA nucleoside, has been performed by means of ab initio calculations at the MP2/6-311++G(d,p)//DFT B3LYP/6-31G(d,p) level of theory. At 298.15 K, all 92 conformers of isolated dT are within a 7.49 kcal/mol Gibbs energy range. Syn orientation for the base and South (S) conformers for the sugar dominate at this temperature: syn/anti = 61.6%:38.4% and S/N = 74.5%:25.5%. However, at 420 K, the majority of conformers contain anti base and the population of North (N) sugars increases: syn/anti = 38.0%:62.0% and S/N = 59.5%:40.5%. The whole conformational parameters (P, chi, gamma, delta, beta, epsilon, nu max) were analyzed as well as the energies of the OH...O intramolecular H-bonds on the basis of nu(OH) stretching vibrations. Convolution of calculated IR spectra of all of the T conformers appears consistent with its low-temperature matrix spectrum (Ivanov et al. Low Temp. Phys. 2003, 29, 809). The maximal discrepancy in frequencies between calculated and experimental spectra is less than 1%. A conclusion was made that for reliable reconstruction of the isolated nucleoside IR spectrum the quasi whole set of conformers should be taken into consideration. In essence, this result opens up a possibility to reconstruct IR spectra of isolated nucleosides at physiological temperatures with rather satisfactory probability.

A posteriori design of crystal contacts to improve the X-ray diffraction properties of a small RNA enzyme

The hairpin ribozyme is a small catalytic RNA comprising two helix-loop-helix domains linked by a four-way helical junction (4WJ). In its most basic form, each domain can be formed independently and reconstituted without a 4WJ to yield an active enzyme. The production of such minimal junctionless hairpin ribozymes is achievable by chemical synthesis, which has allowed structures to be determined for numerous nucleotide variants. However, abasic and other destabilizing core modifications hinder crystallization. This investigation describes the use of a dangling 5'-U to form an intermolecular U.U mismatch, as well as the use of synthetic linkers to tether the loop A and B domains, including (i) a three-carbon propyl linker (C3L) and (ii) a nine-atom triethylene glycol linker (S9L). Both linker constructs demonstrated similar enzymatic activity, but S9L constructs yielded crystals that diffracted to 2.65 A resolution or better. In contrast, C3L variants diffracted to 3.35 A and exhibited a 15 A expansion of the c axis. Crystal packing of the C3L construct showed a paucity of 6(1) contacts, which comprise numerous backbone to 2'-OH hydrogen bonds in junctionless and S9L complexes. Significantly, the crystal packing in minimal structures mimics stabilizing features observed in the 4WJ hairpin ribozyme structure. The results demonstrate how knowledge-based design can be used to improve diffraction and overcome otherwise destabilizing defects.

Transition state analysis of acid-catalyzed dAMP hydrolysis

Multiple kinetic isotope effects (KIEs) on deoxyadenosine monophosphate (dAMP) hydrolysis in 0.1 M HCl were used to determine the transition state (TS) structure and probe its intrinsic reactivity. The experimental KIEs revealed a stepwise (SN1) mechanism, with a discrete oxacarbenium ion intermediate. This is the first direct evidence for the deoxyribosyl oxacarbenium ion in solution. In 50% methanol/0.1 M HCl the products were deoxyribose 5-phosphate (dRMP) and alpha- and beta-methyl dRMP. The alpha-Me-dRMP/beta-Me-dRMP ratio was 8.5:1. Assuming that a free oxacarbenium ion is equally susceptible to nucleophilic attack on either face, this indicated that approximately 20% proceeded through a solvent-separated ion pair complex, or free oxacarbenium ion, a DN+AN mechanism, while approximately 80% of the reaction proceeded through a contact ion pair complex. The oxacarbenium ion lifetime was estimated at 10(-11)-10(-10) s. Computational transition states were found for ANDN, DN*AN, DN*AN, and DN+AN mechanisms using hybrid density functional theory calculations. After taking into account 20% of DN+AN, there was an excellent match of calculated to experimental KIEs for 80% of the reaction having a DN*AN mechanism. That is, C-N bond cleavage is reversible, with dAMP and the {oxacarbenium ion*adenine} complex in equilibrium. The first irreversible step is water attack on the oxacarbenium ion. The calculated 1'-14C KIE for a stepwise mechanism with irreversible C-N bond cleavage (DN*AN) was 1.052, in the range previously associated only with ANDN transition states, and close to the calculated ANDN value, 1.059. The 1'-14C KIE was strongly dependent on the adenine protonation state.