Molecular dynamics simulations and free energy calculations of base flipping in dsRNA

Spontaneous base flipping helps drive Nsp15's preferences in double stranded RNA substrates

Coronaviruses evade detection by the host immune system with the help of the endoribonuclease Nsp15, which regulates levels of viral double stranded RNA by cleaving 3' of uridine (U). While prior structural data shows that to cleave double stranded RNA, Nsp15's target U must be flipped out of the helix, it is not yet understood whether Nsp15 initiates flipping or captures spontaneously flipped bases. We address this gap by designing fluorinated double stranded RNA substrates that allow us to directly relate a U's sequence context to both its tendency to spontaneously flip and its susceptibility to cleavage by Nsp15. Through a combination of nuclease assays, 19F NMR spectroscopy, mass spectrometry, and single particle cryo-EM, we determine that Nsp15 acts most efficiently on unpaired Us, particularly those that are already flipped. Across sequence contexts, we find Nsp15's cleavage efficiency to be directly related to that U's tendency to spontaneously flip. Overall, our findings unify previous characterizations of Nsp15's cleavage preferences, and suggest that activity of Nsp15 during infection is partially driven by bulged or otherwise relatively accessible Us that appear at strategic positions in the viral RNA.

Inosine and its methyl derivatives: Occurrence, biogenesis, and function in RNA

Inosine is one of the most common post-transcriptional modifications. Since its discovery, it has been noted for its ability to contribute to non-Watson-Crick interactions within RNA. Rapidly accumulating evidence points to the widespread generation of inosine through hydrolytic deamination of adenosine to inosine by different classes of adenosine deaminases. Three naturally occurring methyl derivatives of inosine, i.e., 1-methylinosine, 2'-O-methylinosine and 1,2'-O-dimethylinosine are currently reported in RNA modification databases. These modifications are expected to lead to changes in the structure, folding, dynamics, stability and functions of RNA. The importance of the modifications is indicated by the strong conservation of the modifying enzymes across organisms. The structure, binding and catalytic mechanism of the adenosine deaminases have been well-studied, but the underlying mechanism of the catalytic reaction is not very clear yet. Here we extensively review the existing data on the occurrence, biogenesis and functions of inosine and its methyl derivatives in RNA. We also included the structural and thermodynamic aspects of these modifications in our review to provide a detailed and integrated discussion on the consequences of A-to-I editing in RNA and the contribution of different structural and thermodynamic studies in understanding its role in RNA. We also highlight the importance of further studies for a better understanding of the mechanisms of the different classes of deamination reactions. Further investigation of the structural and thermodynamic consequences and functions of these modifications in RNA should provide more useful information about their role in different diseases.

Role of conformational heterogeneity in ligand recognition by viral RNA molecules

Ribonucleic acid (RNA) molecules are known to undergo conformational changes in response to various environmental stimuli including temperature, pH, and ligands. In particular, viral RNA molecules are a key example of conformationally adapting molecules that have evolved to switch between many functional conformations. The transactivation response element (TAR) RNA from the type-1 human immunodeficiency virus (HIV-1) is a viral RNA molecule that is being increasingly explored as a potential therapeutic target due to its role in the viral replication process. In this work, we have studied the dynamics in TAR RNA in apo and liganded states by performing explicit-solvent molecular dynamics (MD) simulations initiated with 27 distinct structures. We determined that the TAR RNA structure is significantly stabilized on ligand binding with especially decreased fluctuations in its two helices. This rigidity is further coupled with the decreased flipping of bulge nucleotides, which were observed to flip more frequently in the absence of ligands. We found that initially-distinct structures of TAR RNA converged to similar conformations on removing ligands. We also report that conformational dynamics in unliganded TAR structures leads to the formation of binding pockets capable of accommodating ligands of various sizes.

Reaction Coordinate and Thermodynamics of Base Flipping in RNA

Base flipping is a key biophysical event involved in recognition of various ligands by ribonucleic acid (RNA) molecules. However, the mechanism of base flipping in RNA remains poorly understood, in part due to the lack of atomistic details on complex rearrangements in neighboring bases. In this work, we applied transition path sampling (TPS) methods to study base flipping in a double-stranded RNA (dsRNA) molecule that is known to interact with RNA-editing enzymes through this mechanism. We obtained an ensemble of 1000 transition trajectories to describe the base-flipping process. We used the likelihood maximization method to determine the refined reaction coordinate (RC) consisting of two collective variables (CVs), a distance and a dihedral angle between nucleotides that form stacking interactions with the flipping base. The free energy profile projected along the refined RC revealed three minima, two corresponding to the initial and final states and one for a metastable state. We suggest that the metastable state likely represents a wobbled conformation of nucleobases observed in NMR studies that is often characterized as the flipped state. The analyses of reactive trajectories further revealed that the base flipping is coupled to a global conformational change in a stem-loop of dsRNA.

Free Energy Landscape and Conformational Kinetics of Hoogsteen Base Pairing in DNA vs. RNA

Genetic information is encoded in the DNA double helix, which, in its physiological milieu, is characterized by the iconical Watson-Crick nucleo-base pairing. Recent NMR relaxation experiments revealed the transient presence of an alternative, Hoogsteen (HG) base pairing pattern in naked DNA duplexes, and estimated its relative stability and lifetime. In contrast with DNA, such structures were not observed in RNA duplexes. Understanding HG base pairing is important because the underlying "breathing" motion between the two conformations can significantly modulate protein binding. However, a detailed mechanistic insight into the transition pathways and kinetics is still missing. We performed enhanced sampling simulation (with combined metadynamics and adaptive force-bias method) and Markov state modeling to obtain accurate free energy, kinetics, and the intermediates in the transition pathway between Watson-Crick and HG base pairs for both naked B-DNA and A-RNA duplexes. The Markov state model constructed from our unbiased MD simulation data revealed previously unknown complex extrahelical intermediates in the seemingly simple process of base flipping in B-DNA. Extending our calculation to A-RNA, for which HG base pairing is not observed experimentally, resulted in relatively unstable, single-hydrogen-bonded, distorted Hoogsteen-like bases. Unlike B-DNA, the transition pathway primarily involved base paired and intrahelical intermediates with transition timescales much longer than that of B-DNA. The seemingly obvious flip-over reaction coordinate (i.e., the glycosidic torsion angle) is unable to resolve the intermediates. Instead, a multidimensional picture involving backbone dihedral angles and distance between hydrogen bond donor and acceptor atoms is required to gain insight into the molecular mechanism.

Exploring the binding mode and thermodynamics of inverse agonists against estrogen-related receptor alpha

Since estrogen-related receptor alpha (ERRα), one of three estrogen-related receptors, displays constitutively active transcriptional activities and important implications in both physiological and pathological processes of breast cancers, ERRα was recently recognized as a new target to fight breast cancers, and regulating the activity of ERRα with inverse agonists has thus become a promising new therapeutic strategy. A few inverse agonists cyclohexylmethyl-(1-p-tolyl-1H-indol-3-ylmethyl)-amine (compound 1), thiadiazoacrylamide (XCT790), and 1-(2,5-diethoxy-benzyl)-3-phenyl-area analogues (compounds 2 and 3) were reported to be capable of targeting ERRα. However, the detailed mechanism by which the inverse agonists deactivate ERRα remains unclear, especially in the aspects of quantitative binding and hot spot residues. Therefore, to gain insights into the interaction modes between inverse agonists and ERRα ligand binding domain, all-atom molecular dynamics (MD) simulations were firstly carried out for the complexes of inverse agonists and ERRα. The binding free energies were then calculated with MM-PBSA method to quantitatively discuss the binding of the inverse agonists with ERRα. The binding affinities were finally decomposed to per-residue contributions to identify the hot spot residues as well as assess their role in the binding mechanism. MD simulations show that the inverse agonists stretch downwards into the ERRα ligand binding pocket (LBP) formed by H3 and H11 helices, and upon the binding H12 adopts a well-defined position in the coactivator groove, where PGC-1α binds to ERRα. Binding energy analysis indicates that compound 3 and XCT790 bind more tightly to ERRα than compounds 1 and 2, and the energy difference mainly results from the contribution of van der Waals interaction. Both binding mode analysis and affinity decomposition per-residue indicate that compound 1, XCT790, and compound 3 have similar binding spectra to ERRα, primarily interacting with the residues of H3, H5, H6/H7 loop, and H11 helix, while compound 2 lacks a significant interaction with the H5 region. The hot spot residues significantly binding to the three inverse agonists in common include Leu324, Phe328, Phe382, Leu398, Phe495, and Leu500. It is essential for an effective inverse agonist to strongly bind with the aromatic ring cluster consisting of Phe328(H3), Phe495(H11), and Phe382(H5/H6 loop) as well as Leu500.

All-atom MD simulations were for the first time carried out for the complexes of inverse agonists and ERRα, and their binding free energies were also calculated with MM-PBSA to quantitatively discuss the binding of the inverse agonists with ERRα.

Differentiation and classification of RNA motifs using small molecule-based pattern recognition

Understanding how to design small molecules that target coding and non-coding RNA has the potential to exponentially increase the number of therapeutically-relevant druggable targets, which are currently mostly proteins. However, there is limited information on the principles at the basis of RNA recognition. In this chapter, we describe a pattern-based technique that can be used for the simultaneous elucidation of RNA motifs and small molecule features for RNA selective recognition, termed Pattern Recognition of RNA by Small Molecules (PRRSM). We provide protocols for the computational design and synthetic preparation of an RNA training set as well as how to perform the assay in plate reader format. Furthermore, we provide details on how to perform and interpret the statistical analysis and indicate possible future extensions of the technique. By combining insights into characteristics of the small molecules and of the RNA that leads to differentiation, PRRSM promises to accelerate the elucidation of the determinants at the basis of RNA recognition.

BAR-based optimum adaptive steered MD for configurational sampling

The equilibrium and nonequilibrium adaptive alchemical free energy simulation methods optimum Bennett's acceptance ratio and optimum crooks' equation (OCE), based on the statistically optimal bidirectional reweighting estimator named Bennett's Acceptance Ratio or Crooks' equation, perform initial sampling in the staging alchemical transformation and then determine the importance rank of different states via the time-derivative of the variance. The method is proven to give speedups compared with the equal time rule. In the current work, we extend the time derivative of variance guided adaptive sampling method to the configurational space, falling in the term of steered MD (SMD). The SMD approach biasing physically meaningful collective variable (CV) such as one dihedral or one distance to pulling the system from one conformational state to another. By minimizing the variance of the free energy differences along the pathway in an optimized way, a new type of adaptive SMD (ASMD) is introduced. As exhibits in the alchemical case, this adaptive sampling method outperforms the traditional equal-time SMD in nonequilibrium stratification. Also, the method gives much more efficient calculation of potential of mean force than the selection criterion-based ASMD scheme, which is proven to be more efficient than traditional SMD. The OCE workflow is periodicity-of-CV dependent while ASMD is not. The performance is demonstrated in a dihedral flipping case and two distance pulling cases, accounting for periodic and nonperiodic CVs, respectively. © 2019 Wiley Periodicals, Inc.

Improvement of RNA Simulations with Torsional Revisions of the AMBER Force Field

Our current knowledge on the unique roles of RNA in cells makes it vital to investigate the properties of RNA systems using computational methods because of the potential pharmaceutical applications. With the continuous advancement of computer technology, it is now possible to study RNA folding. Molecular mechanics calculations are useful in discovering the structural and thermodynamic properties of RNA systems. Yet, the predictions depend on the quality of the RNA force field, which is a set of parameters describing the potential energy of the system. Torsional parameters are one of the terms in a force field that can be revised using physics-based approaches. This chapter focuses on improvements provided by revisions of torsional parameters of the AMBER (Assisted Model Building with Energy Refinement) RNA force field. The theory behind torsional revisions and re-parameterization of several RNA torsions is briefly described. Applications of the revised torsional parameters to study RNA nucleosides, single-stranded RNA tetramers, and RNA repeat expansions are described in detail. It is concluded that RNA force fields require constant revisions and should be benchmarked against diverse RNA systems such as single strands and internal loops in order to test their qualities.

RNA Structural Dynamics As Captured by Molecular Simulations: A Comprehensive Overview

With both catalytic and genetic functions, ribonucleic acid (RNA) is perhaps the most pluripotent chemical species in molecular biology, and its functions are intimately linked to its structure and dynamics. Computer simulations, and in particular atomistic molecular dynamics (MD), allow structural dynamics of biomolecular systems to be investigated with unprecedented temporal and spatial resolution. We here provide a comprehensive overview of the fast-developing field of MD simulations of RNA molecules. We begin with an in-depth, evaluatory coverage of the most fundamental methodological challenges that set the basis for the future development of the field, in particular, the current developments and inherent physical limitations of the atomistic force fields and the recent advances in a broad spectrum of enhanced sampling methods. We also survey the closely related field of coarse-grained modeling of RNA systems. After dealing with the methodological aspects, we provide an exhaustive overview of the available RNA simulation literature, ranging from studies of the smallest RNA oligonucleotides to investigations of the entire ribosome. Our review encompasses tetranucleotides, tetraloops, a number of small RNA motifs, A-helix RNA, kissing-loop complexes, the TAR RNA element, the decoding center and other important regions of the ribosome, as well as assorted others systems. Extended sections are devoted to RNA-ion interactions, ribozymes, riboswitches, and protein/RNA complexes. Our overview is written for as broad of an audience as possible, aiming to provide a much-needed interdisciplinary bridge between computation and experiment, together with a perspective on the future of the field.

A Phenotypic Screen for Functional Mutants of Human Adenosine Deaminase Acting on RNA 1

Adenosine deaminases acting on RNA (ADARs) are RNA-editing enzymes responsible for the conversion of adenosine to inosine at specific locations in cellular RNAs. ADAR1 and ADAR2 are two members of the family that have been shown to be catalytically active. Earlier, we reported a phenotypic screen for the study of human ADAR2 using budding yeast S. cerevisiae as the host system. While this screen has been successfully applied to the study of ADAR2, it failed with ADAR1. Here, we report a new reporter that uses a novel editing substrate and is suitable for the study of ADAR1. We screened plasmid libraries with randomized codons for two important residues in human ADAR1 (G1007 and E1008). The screening results combined with in vitro deamination assays led to the identification of mutants that are more active than the wild type protein. Furthermore, a screen of the ADAR1 E1008X library with a reporter construct bearing an A•G mismatch at the editing site suggests one role for the residue at position 1008 is to sense the identity of the base pairing partner for the editing site adenosine. This work has provided a starting point for future in vitro evolution studies of ADAR1 and led to new insight into ADAR's editing site selectivity.

New Insights into the Biological Role of Mammalian ADARs; the RNA Editing Proteins

The ADAR proteins deaminate adenosine to inosine in double-stranded RNA which is one of the most abundant modifications present in mammalian RNA. Inosine can have a profound effect on the RNAs that are edited, not only changing the base-pairing properties, but can also result in recoding, as inosine behaves as if it were guanosine. In mammals there are three ADAR proteins and two ADAR-related proteins (ADAD) expressed. All have a very similar modular structure; however, both their expression and biological function differ significantly. Only two of the ADAR proteins have enzymatic activity. However, both ADAR and ADAD proteins possess the ability to bind double-strand RNA. Mutations in ADARs have been associated with many diseases ranging from cancer, innate immunity to neurological disorders. Here, we will discuss in detail the domain structure of mammalian ADARs, the effects of RNA editing, and the role of ADARs in human diseases.

Triple helical DNA in a duplex context and base pair opening

It is fundamental to explore in atomic detail the behavior of DNA triple helices as a means to understand the role they might play in vivo and to better engineer their use in genetic technologies, such as antigene therapy. To this aim we have performed atomistic simulations of a purine-rich antiparallel triple helix stretch of 10 base triplets flanked by canonical Watson–Crick double helices. At the same time we have explored the thermodynamic behavior of a flipping Watson–Crick base pair in the context of the triple and double helix. The third strand can be accommodated in a B-like duplex conformation. Upon binding, the double helix changes shape, and becomes more rigid. The triple-helical region increases its major groove width mainly by oversliding in the negative direction. The resulting conformations are somewhere between the A and B conformations with base pairs remaining almost perpendicular to the helical axis. The neighboring duplex regions maintain a B DNA conformation. Base pair opening in the duplex regions is more probable than in the triplex and binding of the Hoogsteen strand does not influence base pair breathing in the neighboring duplex region.

Cooperative gene regulation by microRNA pairs and their identification using a computational workflow

MicroRNAs (miRNAs) are an integral part of gene regulation at the post-transcriptional level. Recently, it has been shown that pairs of miRNAs can repress the translation of a target mRNA in a cooperative manner, which leads to an enhanced effectiveness and specificity in target repression. However, it remains unclear which miRNA pairs can synergize and which genes are target of cooperative miRNA regulation. In this paper, we present a computational workflow for the prediction and analysis of cooperating miRNAs and their mutual target genes, which we refer to as RNA triplexes. The workflow integrates methods of miRNA target prediction; triplex structure analysis; molecular dynamics simulations and mathematical modeling for a reliable prediction of functional RNA triplexes and target repression efficiency. In a case study we analyzed the human genome and identified several thousand targets of cooperative gene regulation. Our results suggest that miRNA cooperativity is a frequent mechanism for an enhanced target repression by pairs of miRNAs facilitating distinctive and fine-tuned target gene expression patterns. Human RNA triplexes predicted and characterized in this study are organized in a web resource at www.sbi.uni-rostock.de/triplexrna/.

Hidden Conformation Events in DNA Base Extrusions: A Generalized Ensemble Path Optimization and Equilibrium Simulation Study

DNA base extrusion is a crucial component of many biomolecular processes. Elucidating how bases are selectively extruded from the interiors of double-strand DNAs is pivotal to accurately understanding and efficiently sampling this general type of conformational transitions. In this work, the on-the-path random walk (OTPRW) method, which is the first generalized ensemble sampling scheme designed for finite-temperature-string path optimizations, was improved and applied to obtain the minimum free energy path (MFEP) and the free energy profile of a classical B-DNA major-groove base extrusion pathway. Along the MFEP, an intermediate state and the corresponding transition state were located and characterized. The MFEP result suggests that a base-plane-elongation event rather than the commonly focused base-flipping event is dominant in the transition state formation portion of the pathway; and the energetic penalty at the transition state is mainly introduced by the stretching of the Watson-Crick base pair. Moreover to facilitate the essential base-plane-elongation dynamics, the surrounding environment of the flipped base needs to be intimately involved. Further taking the advantage of the extended-dynamics nature of the OTPRW Hamiltonian, an equilibrium generalized ensemble simulation was performed along the optimized path; and based on the collected samples, several base-flipping (opening) angle collective variables were evaluated. In consistence with the MFEP result, the collective variable analysis result reveals that none of these commonly employed flipping (opening) angles alone can adequately represent the base extrusion pathway, especially in the pre-transition-state portion. As further revealed by the collective variable analysis, the base-pairing partner of the extrusion target undergoes a series of in-plane rotations to facilitate the base-plane-elongation dynamics. A base-plane rotation angle is identified to be a possible reaction coordinate to represent these in-plane rotations. Notably, these in-plane rotation motions may play a pivotal role in determining the base extrusion selectivity.

Partial base flipping is sufficient for strand slippage near DNA duplex termini

Strand slippage is a structural mechanism by which insertion-deletion (indel) mutations are introduced during replication by polymerases. Three-dimensional atomic-resolution structural pathways are still not known for the decades-old template slippage description. The dynamic nature of the process and the higher energy intermediates involved increase the difficulty of studying these processes experimentally. In the present study, restrained and unrestrained molecular dynamics simulations, carried out using multiple nucleic acid force fields, are used to demonstrate that partial base-flipping can be sufficient for strand slippage at DNA duplex termini. Such strand slippage can occur in either strand, i.e. near either the 3' or the 5' terminus of a DNA strand, which suggests that similar structural flipping mechanisms can cause both primer and template slippage. In the repetitive mutation hot-spot sequence studied, non-canonical base-pairing with exposed DNA groove atoms of a neighboring G:C base-pair stabilizes a partially flipped state of the cytosine. For its base-pair partner guanine, a similar partially flipped metastable intermediate was not detected, and the propensity for sustained slippage was also found to be lower. This illustrates that a relatively small metastable DNA structural distortion in polymerase active sites could allow single base insertion or deletion mutations to occur, and stringent DNA groove molecular recognition may be required to maintain intrinsic DNA polymerase fidelity. The implications of a close relationship between base-pair dissociation, base unstacking, and strand slippage are discussed in the context of sequence dependence of indel mutations.

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.

Sequence-dependent structural dynamics of primate adenosine-to-inosine editing substrates

Humans have the highest level of adenosine-to-inosine (A-to-I) editing amongst primates, yet the reasons for this difference remain unclear. Sequence analysis of the Alu Sg elements (A-to-I RNA substrates) corresponding to the Nup50 gene in human, chimp, and rhesus reveals subtle sequence variations surrounding the edit sites. We have developed three constructs that represent human (HuAp5), chimp (ChAp5), and rhesus (RhAp5) Nup50 Alu Sg A-to-I editing substrates. Here, 2-aminopurine (2-Ap) was substituted for edited adenosine (A5) so as to monitor the fluorescence intensity with respect to temperature. UV and steady-state fluorescence (SSF) T(M) plots indicate that local and global unfolding are coincident, with the human construct displaying a T(M) of approximately 70°C, compared to 60°C for chimp and 54°C for rhesus. However, time-resolved fluorescence (TRF) resolves three different fluorescence lifetimes that we assign to folded, intermediate(s), and unfolded states. The TRF data fit well to a two-intermediate model, whereby both intermediates (M, J) are in equilibrium with each other, and the folded/unfolded states. Our model suggests that, at 37°C, human state J and the folded state will be the most heavily populated in comparison to the other primate constructs. In order for adenosine deaminase acting on RNA (ADAR) to efficiently dock, a stable duplex must be present that corresponds to the human construct, globally. Next, the enzyme must "flip out" the base of interest to facilitate the A-to-I conversion; a nucleotide in an intermediate-like position would enhance this conformational change. Our experiments demonstrate that subtle variations in RNA sequence might contribute to the high A-to-I editing levels found in humans.

Mechanistic insights into editing-site specificity of ADARs

Adenosine deaminases that act on RNA (ADARs) deaminate adenosines in dsRNA to produce inosines. ADARs are essential in mammals and are particularly important in the nervous system. Altered levels of adenosine-to-inosine (A-to-I) editing are observed in several diseases. The extent to which an adenosine is edited depends on sequence context. Human ADAR2 (hADAR2) has 5' and 3' neighbor preferences, but which amino acids mediate these preferences, and by what mechanism, is unknown. We performed a screen in yeast to identify mutations in the hADAR2 catalytic domain that allow editing of an adenosine within a disfavored triplet. Binding affinity, catalytic rate, base flipping, and preferences were monitored to understand the effects of the mutations on ADAR reactivity. Our data provide information on the amino acids that affect preferences and point to a conserved loop as being of key importance. Unexpectedly, our data suggest that hADAR2's preferences derive from differential base flipping rather than from direct recognition of neighboring bases. Our studies set the stage for understanding the basis of altered editing levels in disease and for developing therapeutic reagents.

Recent Developments and Applications of the CHARMM force fields

Empirical force fields commonly used to describe the condensed phase properties of complex systems such as biological macromolecules are continuously being updated. Improvements in quantum mechanical (QM) methods used to generate target data, availability of new experimental target data, incorporation of new classes of compounds and new theoretical developments (eg. polarizable methods) make force-field development a dynamic domain of research. Accordingly, a number of improvements and extensions of the CHARMM force fields have occurred over the years. The objective of the present review is to provide an up-to-date overview of the CHARMM force fields. A limited presentation on the historical aspects of force fields will be given, including underlying methodologies and principles, along with a brief description of the strategies used for parameter development. This is followed by information on the CHARMM additive and polarizable force fields, including examples of recent applications of those force fields.

Molecular Mechanics Investigation of an Adenine-Adenine Non-Canonical Pair Conformational Change

Conformational changes are important in RNA for binding and catalysis and understanding these changes is important for understanding how RNA functions. Computational techniques using all-atom molecular models can be used to characterize conformational changes in RNA. These techniques are applied to an RNA conformational change involving a single base pair within a nine base pair RNA duplex. The Adenine-Adenine (AA) non-canonical pair in the sequence 5'GGUGAAGGCU3' paired with 3'PCCGAAGCCG5', where P is Purine, undergoes conformational exchange between two conformations on the timescale of tens of microseconds, as demonstrated in a previous NMR solution structure [Chen, G., et al., Biochemistry, 2006. 45: 6889-903]. The more populated, major, conformation was estimated to be 0.5 to 1.3 kcal/mol more stable at 30 °C than the less populated, minor, conformation. Both conformations are trans-Hoogsteen/sugar edge pairs, where the interacting edges on the adenines change with the conformational change. Targeted Molecular Dynamics (TMD) and Nudged Elastic Band (NEB) were used to model the pathway between the major and minor conformations using the AMBER software package. The adenines were predicted to change conformation via intermediates in which they are stacked as opposed to hydrogen-bonded. The predicted pathways can be described by an improper dihedral angle reaction coordinate. Umbrella sampling along the reaction coordinate was performed to model the free energy profile for the conformational change using a total of 1800 ns of sampling. Although the barrier height between the major and minor conformations was reasonable, the free energy difference between the major and minor conformations was the opposite of that expected based on the NMR experiments. Variations in the force field applied did not improve the misrepresentation of the free energies of the major and minor conformations. As an alternative, the Molecular Mechanics Poisson-Boltzmann Surface Area (MMPBSA) approximation was applied to predict free energy differences between the two conformations using a total of 800 ns of sampling. MM-PBSA also incorrectly predicted the major conformation to be higher in free energy than the minor conformation.

Impact of 2'-hydroxyl sampling on the conformational properties of RNA: update of the CHARMM all-atom additive force field for RNA

Here, we present an update of the CHARMM27 all-atom additive force field for nucleic acids that improves the treatment of RNA molecules. The original CHARMM27 force field parameters exhibit enhanced Watson-Crick base pair opening which is not consistent with experiment, whereas analysis of molecular dynamics (MD) simulations show the 2'-hydroxyl moiety to almost exclusively sample the O3' orientation. Quantum mechanical (QM) studies of RNA related model compounds indicate the energy minimum associated with the O3' orientation to be too favorable, consistent with the MD results. Optimization of the dihedral parameters dictating the energy of the 2'-hydroxyl proton targeting the QM data yielded several parameter sets, which sample both the base and O3' orientations of the 2'-hydroxyl to varying degrees. Selection of the final dihedral parameters was based on reproduction of hydration behavior as related to a survey of crystallographic data and better agreement with experimental NMR J-coupling values. Application of the model, designated CHARMM36, to a collection of canonical and noncanonical RNA molecules reveals overall improved agreement with a range of experimental observables as compared to CHARMM27. The results also indicate the sensitivity of the conformational heterogeneity of RNA to the orientation of the 2'-hydroxyl moiety and support a model whereby the 2'-hydroxyl can enhance the probability of conformational transitions in RNA.

ADAR proteins: structure and catalytic mechanism

Since the discovery of the adenosine deaminase (ADA) acting on RNA (ADAR) family of proteins in 1988 (Bass and Weintraub, Cell 55:1089-1098, 1988) (Wagner et al. Proc Natl Acad Sci U S A 86:2647-2651, 1989), we have learned much about their structure and catalytic mechanism. However, much about these enzymes is still unknown, particularly regarding the selective recognition and processing of specific adenosines within substrate RNAs. While a crystal structure of the catalytic domain of human ADAR2 has been solved, we still lack structural data for an ADAR catalytic domain bound to RNA, and we lack any structural data for other ADARs. However, by analyzing the structural data that is available along with similarities to other deaminases, mutagenesis and other biochemical experiments, we have been able to advance the understanding of how these fascinating enzymes function.

On the molecular basis of uracil recognition in DNA: comparative study of T-A versus U-A structure, dynamics and open base pair kinetics

Uracil (U) can be found in DNA as a mismatch paired either to adenine (A) or to guanine (G). Removal of U from DNA is performed by a class of enzymes known as uracil–DNA–glycosylases (UDG). Recent studies suggest that recognition of U–A and U–G mismatches by UDG takes place via an extra-helical mechanism. In this work, we use molecular dynamics simulations to analyze the structure, dynamics and open base pair kinetics of U–A base pairs relative to their natural T–A counterpart in 12 dodecamers. Our results show that the presence of U does not alter the local conformation of B-DNA. Breathing dynamics and base pair closing kinetics are only weakly dependent on the presence of U versus T, with open T–A and U–A pairs lifetimes in the nanosecond timescale. Additionally, we observed spontaneous base flipping in U–A pairs. We analyze the structure and dynamics for this event and compare the results to available crystallographic data of open base pair conformations. Our results are in agreement with both structural and kinetic data derived from NMR imino proton exchange measurements, providing the first detailed description at the molecular level of elusive events such as spontaneous base pair opening and flipping in mismatched U–A sequences in DNA. Based on these results, we propose that base pair flipping can occur spontaneously at room temperature via a 3-step mechanism with an open base pair intermediate. Implications for the molecular basis of U recognition by UDG are discussed.

Expanded RNA-binding activities of mammalian Argonaute 2

Mammalian Argonaute 2 (Ago2) protein associates with microRNAs (miRNAs) or small interfering RNAs (siRNAs) forming RNA-induced silencing complexes (RISCs/miRNPs). In the present work, we characterize the RNA-binding and nucleolytic activity of recombinant mouse Ago2. Our studies show that recombinant mouse Ago2 binds efficiently to miRNAs forming active RISC. Surprisingly, we find that recombinant mouse Ago2 forms active RISC using pre-miRNAs or long unstructured single stranded RNAs as guides. Furthermore, we demonstrate that, in vivo, endogenous human Ago2 binds directly to pre-miRNAs independently of Dicer, and that Ago2:pre-miRNA complexes are found both in the cytoplasm and in the nucleus of human cells.

Multiple determinants within iron-responsive elements dictate iron regulatory protein binding and regulatory hierarchy

Iron regulatory proteins (IRPs) are iron-regulated RNA binding proteins that, along with iron-responsive elements (IREs), control the translation of a diverse set of mRNA with 5' IRE. Dysregulation of IRP action causes disease with etiology that may reflect differential control of IRE-containing mRNA. IREs are defined by a conserved stem-loop structure including a midstem bulge at C8 and a terminal CAGUGH sequence that forms an AGU pseudo-triloop and N19 bulge. C8 and the pseudo-triloop nucleotides make the majority of the 22 identified bonds with IRP1. We show that IRP1 binds 5' IREs in a hierarchy extending over a ninefold range of affinities that encompasses changes in IRE binding affinity observed with human L-ferritin IRE mutants. The limits of this IRE binding hierarchy are predicted to arise due to small differences in binding energy (e.g., equivalent to one H-bond). We demonstrate that multiple regions of the IRE stem not predicted to contact IRP1 help establish the binding hierarchy with the sequence and structure of the C8 region displaying a major role. In contrast, base-pairing and stacking in the upper stem region proximal to the terminal loop had a minor role. Unexpectedly, an N20 bulge compensated for the lack of an N19 bulge, suggesting the existence of novel IREs. Taken together, we suggest that a regulatory binding hierarchy is established through the impact of the IRE stem on the strength, not the number, of bonds between C8 or pseudo-triloop nucleotides and IRP1 or through their impact on an induced fit mechanism of binding.

An experimentally guided umbrella sampling protocol for biomolecules

We present a simple method for utilizing experimental data to improve the efficiency of numerical calculations of free energy profiles from molecular dynamics simulations. The method involves umbrella sampling simulations with restraining potentials based on a known approximate estimate of the free energy profile derived solely from experimental data. The use of the experimental data results in optimal restraining potentials, guides the simulation along relevant pathways, and decreases overall computational time. In demonstration of the method, two systems are showcased. First, guided, unguided (regular) umbrella sampling simulations and exhaustive sampling simulations are compared to each other in the calculation of the free energy profile for the distance between the ends of a pentapeptide. The guided simulation use restraints based on a simulated "experimental" potential of mean force of the end-to-end distance that would be measured by fluorescence resonance energy transfer (obtained from exhaustive sampling). Statistical analysis shows a dramatic improvement in efficiency for a 5 window guided umbrella sampling over 5 and 17 window unguided umbrella sampling simulations. Moreover, the form of the potential of mean force for the guided simulations evolves, as one approaches convergence, along the same milestones as the extensive simulations, but exponentially faster. Second, the method is further validated by replicating the forced unfolding pathway of the titin I27 domain using guiding umbrella sampling potentials determined from actual single molecule pulling data. Comparison with unguided umbrella sampling reveals that the use of guided sampling encourages unfolding simulations to converge faster to a forced unfolding pathway that agrees with previous results and produces a more accurate potential of mean force.

Molecular dynamics simulations of Factor Xa: insight into conformational transition of its binding subsites

Protein flexibility and conformational diversity is well known to be a key characteristic of the function of many proteins. Human blood coagulation proteins have multiple substrates, and various protein-protein interactions are required for the smooth functioning of the coagulation cascade to maintain blood hemostasis. To address how a protein may cope with multiple interactions with its structurally diverse substrates and the accompanied structural changes that may drive these changes, we studied human Factor X. We employed 20 ns of molecular dynamics (MD) and steered molecular dynamics (SMD) simulations on two different conformational forms of Factor X, open and closed, and observed an interchangeable conformational transition from one to another. This work also demonstrates the roles of various aromatic residues involved in aromatic-aromatic interactions, which make this dynamic transition possible.

Recent advances in the study of nucleic acid flexibility by molecular dynamics

The recent use of molecular dynamics (MD) simulations to study flexibility of nucleic acids has been reviewed from an analysis of the publications appearing in the past two years (from 2005 till date). Despite the existence of some unsolved problems in the methodologies, these years have been witness to major advances in the field. Based on a critical review of the most recent contributions, excitement exists on the expected evolution of the field in the next years.

Coupling of fast and slow modes in the reaction pathway of the minimal hammerhead ribozyme cleavage

By employing classical molecular dynamics, correlation analysis of coupling between slow and fast dynamical modes, and free energy (umbrella) sampling using classical as well as mixed quantum mechanics molecular mechanics force fields, we uncover a possible pathway for phosphoryl transfer in the self-cleaving reaction of the minimal hammerhead ribozyme. The significance of this pathway is that it initiates from the minimal hammerhead crystal structure and describes the reaction landscape as a conformational rearrangement followed by a covalent transformation. The delineated mechanism is catalyzed by two metal (Mg(2+)) ions, proceeds via an in-line-attack by CYT 17 O2' on the scissile phosphorous (ADE 1.1 P), and is therefore consistent with the experimentally observed inversion configuration. According to the delineated mechanism, the coupling between slow modes involving the hammerhead backbone with fast modes in the cleavage site appears to be crucial for setting up the in-line nucleophilic attack.

Molecular dynamics simulations of RNA: an in silico single molecule approach

RNA molecules are now known to be involved in the processing of genetic information at all levels, taking on a wide variety of central roles in the cell. Understanding how RNA molecules carry out their biological functions will require an understanding of structure and dynamics at the atomistic level, which can be significantly improved by combining computational simulation with experiment. This review provides a critical survey of the state of molecular dynamics (MD) simulations of RNA, including a discussion of important current limitations of the technique and examples of its successful application. Several types of simulations are discussed in detail, including those of structured RNA molecules and their interactions with the surrounding solvent and ions, catalytic RNAs, and RNA-small molecule and RNA-protein complexes. Increased cooperation between theorists and experimentalists will allow expanded judicious use of MD simulations to complement conceptually related single molecule experiments. Such cooperation will open the door to a fundamental understanding of the structure-function relationships in diverse and complex RNA molecules. .

Conformational transitions in RNA single uridine and adenosine bulge structures: a molecular dynamics free energy simulation study

Extra unmatched nucleotides (single base bulges) are common structural motifs in folded RNA molecules and can participate in RNA-ligand binding and RNA tertiary structure formation. Often these processes are associated with conformational transitions in the bulge region such as flipping out of the bulge base from an intrahelical stacked toward a looped out state. Knowledge of the flexibility of bulge structures and energetics of conformational transitions is an important prerequisite to better understand the function of this RNA motif. Molecular dynamics simulations were performed on single uridine and adenosine bulge nucleotides at the center of eight basepair RNA molecules and indicated larger flexibility of the bulge bases compared to basepaired regions. The umbrella sampling method was applied to study the bulge base looping out process and accompanying conformational and free energy changes. Looping out toward the major groove resulted in partial disruption of adjacent basepairs and was found to be less favorable compared to looping out toward the minor groove. For both uridine and adenosine bulges, a positive free energy change for full looping out was obtained which was approximately 1.5 kcal mol-1 higher in the case of the adenosine compared to the uridine bulge system. The simulations also indicated stable partially looped out states with the bulge bases located in the RNA minor groove and forming base triples with 5'-neighboring basepairs. In the case of the uridine bulge this state was more stable than the intrahelical stacked bulge structure. Induced looping out toward the minor groove involved crossing of an energy barrier of approximately 3.5 kcal mol-1 before reaching the base triple state. A continuum solvent analysis of intermediate bulge states indicated that electrostatic interactions stabilize looped out and base triple states, whereas van der Waals interactions and nonpolar contributions favor the stacked bulge conformation.