Free-energy landscape of RNA hairpins constructed via dihedral angle principal component analysis

Energy landscapes of homopolymeric RNAs revealed by deep unsupervised learning

Conformational dynamics of RNA plays important roles in a variety of cellular functions such as transcriptional regulation, catalysis, scaffolding, and sensing. Recently, RNAs with low-complexity sequences have been shown to phase separate and form condensate phases similar to lowcomplexity protein domains. The affinity for phase separation and the material characteristics of RNA condensates are strongly dependent on sequence composition and patterning. We hypothesize that differences in the affinities for RNA phase separation can be uncovered by studying sequence-dependent conformational dynamics of single RNA chains. To this end, we have employed atomistic simulations and deep dimensionality reduction techniques to map temperature-dependent conformational free energy landscapes for 20 base-long homopolymeric RNA sequences: poly(U), poly(G), poly(C), and poly(A). The energy landscapes of homopolymeric RNAs reveal a plethora of metastable states with qualitatively different populations stemming from differences in base chemistry. Through detailed analysis of base, phosphate, and sugar interactions, we show that experimentally observed temperature-driven shifts in metastable state populations align with experiments on RNA phase transitions. Specifically, we find that the thermodynamics of unfolding of homopolymeric RNA follows the poly(G) > poly(A) > poly(C) > poly(U) order of stability, mirroring the propensity of RNA to form condensates. To conclude, this work shows that at least for homopolymeric RNA sequences the single-chain conformational dynamics contains sufficient information for predicting and quantifying condensate forming affinities of RNAs. Thus, we anticipate that atomically detailed studies of temeprature -dependent energy landscapes of RNAs will be a useful guide for understanding the propensity of various RNA molecules to form condensates.

Log-periodic oscillations as real-time signatures of hierarchical dynamics in proteins

The time-dependent relaxation of a dynamical system may exhibit a power-law behavior that is superimposed by log-periodic oscillations. D. Sornette [Phys. Rep. 297, 239 (1998)] showed that this behavior can be explained by a discrete scale invariance of the system, which is associated with discrete and equidistant timescales on a logarithmic scale. Examples include such diverse fields as financial crashes, random diffusion, and quantum topological materials. Recent time-resolved experiments and molecular dynamics simulations suggest that discrete scale invariance may also apply to hierarchical dynamics in proteins, where several fast local conformational changes are a prerequisite for a slow global transition to occur. Employing entropy-based timescale analysis and Markov state modeling to a simple one-dimensional hierarchical model and biomolecular simulation data, it is found that hierarchical systems quite generally give rise to logarithmically spaced discrete timescales. By introducing a one-dimensional reaction coordinate that collectively accounts for the hierarchically coupled degrees of freedom, the free energy landscape exhibits a characteristic staircase shape with two metastable end states, which causes the log-periodic time evolution of the system. The period of the log-oscillations reflects the effective roughness of the energy landscape and can, in simple cases, be interpreted in terms of the barriers of the staircase landscape.

Classification of MLH1 Missense VUS Using Protein Structure-Based Deep Learning-Ramachandran Plot-Molecular Dynamics Simulations Method

Pathogenic variation in DNA mismatch repair (MMR) gene MLH1 is associated with Lynch syndrome (LS), an autosomal dominant hereditary cancer. Of the 3798 MLH1 germline variants collected in the ClinVar database, 38.7% (1469) were missense variants, of which 81.6% (1199) were classified as Variants of Uncertain Significance (VUS) due to the lack of functional evidence. Further determination of the impact of VUS on MLH1 function is important for the VUS carriers to take preventive action. We recently developed a protein structure-based method named "Deep Learning-Ramachandran Plot-Molecular Dynamics Simulation (DL-RP-MDS)" to evaluate the deleteriousness of MLH1 missense VUS. The method extracts protein structural information by using the Ramachandran plot-molecular dynamics simulation (RP-MDS) method, then combines the variation data with an unsupervised learning model composed of auto-encoder and neural network classifier to identify the variants causing significant change in protein structure. In this report, we applied the method to classify 447 MLH1 missense VUS. We predicted 126/447 (28.2%) MLH1 missense VUS were deleterious. Our study demonstrates that DL-RP-MDS is able to classify the missense VUS based solely on their impact on protein structure.

Exploring the viral protease inhibitor space driven by consensus scoring-based virtual screening

Dengue fever presents a major health concern, and the lack of an effective vaccine or definite therapeutic regimen has led the research community to identify safe-by-design potential targets for drug discovery. Since the association of the NS2B co-factor with the protease domain of NS3 is imperative for the catalytic activity of the enzyme complex, inhibitors blocking their interaction could provide an alternative strategy to combat the dengue virus. In this context, the present study is aimed at exploring computer-assisted modeling of significant physicochemical features required for the inhibition of the dengue virus protease complex. First of all, alanine scanning was utilized to map hot spot residues critical for the association of the two subunits, NS2B and NS3pro, by studying their energy profiles. Then, consensus score-based virtual screening was performed to search through the commercially available chemical datasets. After screening, 1,575 small molecules were moved forward into docking studies to investigate their interactions with crucial interfacial residues (i.e., Tyr23, Lys26, Phe46, and Leu58), with only 233 molecules passing that stage. The top 30 molecules were selected based on a detailed profile of intermolecular interactions. After that, the top five molecules were selected for detailed mechanistic studies via molecular dynamics simulations followed by subsequent binding free energy calculations, principal component analysis in conjunction with free energy landscape. To the best of our knowledge, this is the first systematic and comprehensive investigation to identify protein-protein interaction blockers against the target protein at such a large scale, using integrated computational tools. Our results highlight the enhanced stability and good binding affinities towards the target protein of these compounds, which might act as new scaffolds for NS2B-NS3 protease inhibition. Future studies will be directed to explore the detailed atomistic-based structural and energetic framework of the mutation-induced affinity change between the protease domain of the DENV-2 NS3 protein and its cofactor NS2B. The detailed insight in turn might suggest precise and focused targeted points for the structure-based drug design but the computational cost may be a challenge.

Graphical abstract

Supplementary Information

The online version contains supplementary material available at 10.1007/s40203-023-00174-0.

Selecting Features for Markov Modeling: A Case Study on HP35

Markov state models represent a popular means to interpret molecular dynamics trajectories in terms of memoryless transitions between metastable conformational states. To provide a mechanistic understanding of the considered biomolecular process, these states should reflect structurally distinct conformations and ensure a time scale separation between fast intrastate and slow interstate dynamics. Adopting the folding of villin headpiece (HP35) as a well-established model problem, here we discuss the selection of suitable input coordinates or "features", such as backbone dihedral angles and interresidue distances. We show that dihedral angles account accurately for the structure of the native energy basin of HP35, while the unfolded region of the free energy landscape and the folding process are best described by tertiary contacts of the protein. To construct a contact-based model, we consider various ways to define and select contact distances and introduce a low-pass filtering of the feature trajectory as well as a correlation-based characterization of states. Relying on input data that faithfully account for the mechanistic origin of the studied process, the states of the resulting Markov model are clearly discriminated by the features, describe consistently the hierarchical structure of the free energy landscape, and─as a consequence─correctly reproduce the slow time scales of the process.

Integration of deep learning with Ramachandran plot molecular dynamics simulation for genetic variant classification

Functional classification of genetic variants is a key for their clinical applications in patient care. However, abundant variant data generated by the next-generation DNA sequencing technologies limit the use of experimental methods for their classification. Here, we developed a protein structure and deep learning (DL)-based system for genetic variant classification, DL-RP-MDS, which comprises two principles: 1) Extracting protein structural and thermodynamics information using the Ramachandran plot-molecular dynamics simulation (RP-MDS) method, 2) combining those data with an unsupervised learning model of auto-encoder and a neural network classifier to identify the statistical significance patterns of the structural changes. We observed that DL-RP-MDS provided higher specificity than over 20 widely used in silico methods in classifying the variants of three DNA damage repair genes: TP53, MLH1, and MSH2. DL-RP-MDS offers a powerful platform for high-throughput genetic variant classification. The software and online application are available at https://genemutation.fhs.um.edu.mo/DL-RP-MDS/.

On the Uses of PCA to Characterise Molecular Dynamics Simulations of Biological Macromolecules: Basics and Tips for an Effective Use

Principal Component Analysis (PCA) is a procedure widely used to examine data collected from molecular dynamics simulations of biological macromolecules. It allows for greatly reducing the dimensionality of their configurational space, facilitating further qualitative and quantitative analysis. Its simplicity and relatively low computational cost explain its extended use. However, a judicious implementation of PCA requires the knowledge of its theoretical grounds as well as its weaknesses and capabilities. In this article, we review these issues and discuss several strategies developed over the last years to mitigate the main PCA flaws and enhance the reproducibility of its results.

Unsupervised Learning Methods for Molecular Simulation Data

Unsupervised learning is becoming an essential tool to analyze the increasingly large amounts of data produced by atomistic and molecular simulations, in material science, solid state physics, biophysics, and biochemistry. In this Review, we provide a comprehensive overview of the methods of unsupervised learning that have been most commonly used to investigate simulation data and indicate likely directions for further developments in the field. In particular, we discuss feature representation of molecular systems and present state-of-the-art algorithms of dimensionality reduction, density estimation, and clustering, and kinetic models. We divide our discussion into self-contained sections, each discussing a specific method. In each section, we briefly touch upon the mathematical and algorithmic foundations of the method, highlight its strengths and limitations, and describe the specific ways in which it has been used-or can be used-to analyze molecular simulation data.

Exact Topology of the Dynamic Probability Surface of an Activated Process by Persistent Homology

To gain insight into the reaction mechanism of activated processes, we introduce an exact approach for quantifying the topology of high-dimensional probability surfaces of the underlying dynamic processes. Instead of Morse indexes, we study the homology groups of a sequence of superlevel sets of the probability surface over high-dimensional configuration spaces using persistent homology. For alanine-dipeptide isomerization, a prototype of activated processes, we identify locations of probability peaks and connecting ridges, along with measures of their global prominence. Instead of a saddle point, the transition state ensemble (TSE) of conformations is at the most prominent probability peak after reactants/products, when proper reaction coordinates are included. Intuition-based models, even those exhibiting a double-well, fail to capture the dynamics of the activated process. Peak occurrence, prominence, and locations can be distorted upon subspace projection. While principal component analysis accounts for conformational variance, it inflates the complexity of the surface topology and destroys the dynamic properties of the topological features. In contrast, TSE emerges naturally as the most prominent peak beyond the reactant/product basins, when projected to a subspace of minimum dimension containing the reaction coordinates. Our approach is general and can be applied to investigate the topology of high-dimensional probability surfaces of other activated processes.

Molecular Mechanism of Small-Molecule Inhibitors in Blocking the PD-1/PD-L1 Pathway through PD-L1 Dimerization

Programmed cell death-1 (PD-1), which is a molecule involved in the inhibitory signal in the immune system and is important due to blocking of the interactions between PD-1 and programmed cell death ligand-1 (PD-L1), has emerged as a promising immunotherapy for treating cancer. In this work, molecular dynamics simulations were performed on complex systems consisting of the PD-L1 dimer with (S)-BMS-200, (R)-BMS-200 and (MOD)-BMS-200 (i.e., S, R and MOD systems) to systematically evaluate the inhibitory mechanism of BMS-200-related small-molecule inhibitors in detail. Among them, (MOD)-BMS-200 was modified from the original (S)-BMS-200 by replacing the hydroxyl group with a carbonyl to remove its chirality. Binding free energy analysis indicates that BMS-200-related inhibitors can promote the dimerization of PD-L1. Meanwhile, no significant differences were observed between the S and MOD systems, though the R system exhibited a slightly higher energy. Residue energy decomposition, nonbonded interaction, and contact number analyses show that the inhibitors mainly bind with the C, F and G regions of the PD-L1 dimer, while nonpolar interactions of key residues Ile54, Tyr56, Met115, Ala121 and Tyr123 on both PD-L1 monomers are the dominant binding-related stability factors. Furthermore, compared with (S)-BMS-200, (R)-BMS-200 is more likely to form hydrogen bonds with charged residues. Finally, free energy landscape and protein-protein interaction analyses show that the key residues of the PD-L1 dimer undergo remarkable conformational changes induced by (S)-BMS-200, which boosts its intimate interactions. This systematic investigation provides a comprehensive molecular insight into the ligand recognition process, which will benefit the design of new small-molecule inhibitors targeting PD-L1 for use in anticancer therapy.

Molecular Dynamics Simulations Unveil the Basis of the Sequential Binding of RsmE to the Noncoding RNA RsmZ

CsrA/RsmE are dimeric proteins that bind to targeted mRNAs repressing translation. This mechanism modulates several metabolic pathways and allows bacteria to efficiently adjust their responses to environmental changes. In turn, small RNAs (sRNA) such as CsrB or RsmZ, restore translation by sequestering CsrA/RsmE dimers. Thus, these molecules act in tandem as a gene-expression regulatory system. Recently, a combined NMR-EPR approach solved the structure of part of RsmZ of Pseudomonas fluorescens, attached to three RsmE dimers. The study demonstrated that RsmE assembles onto RsmZ following a specific sequential order. The reasons underlying this peculiar behavior are still unclear. Here, we present a molecular dynamics analysis that explores the conformational diversity of RsmZ and RsmZ-RsmE complexes. The results reveal a clear pattern regarding the exposure of the alternative GGA binding motifs of RsmZ. This pattern is tuned by the attachment of RsmE dimers. Altogether, the observations provide a simple and convincing explanation for the order observed in the sequestration of RsmE dimers. Typical structures for RsmZ and RsmZ-RsmE complexes have been identified. Their characteristics concerning the exposure of the GGA sequences are presented and their most significant interactions are described.

Assessing Force Field Potential Energy Function Accuracy via a Multipolar Description of Atomic Electrostatic Interactions in RNA

Computational investigations of RNA properties often rely on a molecular mechanical approach to define molecular potential energy. Force fields for RNA typically employ a point charge model of electrostatics, which does not provide a realistic quantum-mechanical picture. In reality, electron distributions around nuclei are not spherically symmetric and are geometry dependent. A multipole expansion method which allows for incorporation of polarizability and anisotropy in a force field is described, and its applicability to modeling the behavior of RNA molecules is investigated. Transferability of the model, critical for force field development, is also investigated.

Perspective: Identification of collective variables and metastable states of protein dynamics

The statistical analysis of molecular dynamics simulations requires dimensionality reduction techniques, which yield a low-dimensional set of collective variables (CVs) {x _i } = x that in some sense describe the essential dynamics of the system. Considering the distribution P( x ) of the CVs, the primal goal of a statistical analysis is to detect the characteristic features of P( x ), in particular, its maxima and their connection paths. This is because these features characterize the low-energy regions and the energy barriers of the corresponding free energy landscape ΔG( x ) = -k _B T ln P( x ), and therefore amount to the metastable states and transition regions of the system. In this perspective, we outline a systematic strategy to identify CVs and metastable states, which subsequently can be employed to construct a Langevin or a Markov state model of the dynamics. In particular, we account for the still limited sampling typically achieved by molecular dynamics simulations, which in practice seriously limits the applicability of theories (e.g., assuming ergodicity) and black-box software tools (e.g., using redundant input coordinates). We show that it is essential to use internal (rather than Cartesian) input coordinates, employ dimensionality reduction methods that avoid rescaling errors (such as principal component analysis), and perform density based (rather than k-means-type) clustering. Finally, we briefly discuss a machine learning approach to dimensionality reduction, which highlights the essential internal coordinates of a system and may reveal hidden reaction mechanisms.

Principal component analysis on a torus: Theory and application to protein dynamics

A dimensionality reduction method for high-dimensional circular data is developed, which is based on a principal component analysis (PCA) of data points on a torus. Adopting a geometrical view of PCA, various distance measures on a torus are introduced and the associated problem of projecting data onto the principal subspaces is discussed. The main idea is that the (periodicity-induced) projection error can be minimized by transforming the data such that the maximal gap of the sampling is shifted to the periodic boundary. In a second step, the covariance matrix and its eigendecomposition can be computed in a standard manner. Adopting molecular dynamics simulations of two well-established biomolecular systems (Aib₉ and villin headpiece), the potential of the method to analyze the dynamics of backbone dihedral angles is demonstrated. The new approach allows for a robust and well-defined construction of metastable states and provides low-dimensional reaction coordinates that accurately describe the free energy landscape. Moreover, it offers a direct interpretation of covariances and principal components in terms of the angular variables. Apart from its application to PCA, the method of maximal gap shifting is general and can be applied to any other dimensionality reduction method for circular data.

Effect of Nucleotide State on the Protofilament Conformation of Tubulin Octamers

At the molecular level, the dynamic instability (random growth and shrinkage) of the microtubule (MT) is driven by the nucleotide state (GTP vs GDP) in the β subunit of the tubulin dimers at the MT cap. Here, we use large-scale molecular dynamics (MD) simulations and normal-mode analysis (NMA) to characterize the effect of a single GTP cap layer on tubulin octamers composed of two neighboring protofilaments (PFs). We utilize recently reported high-resolution structures of dynamic MTs to simulate a GDP octamer both with and without a single GTP cap layer. We perform multiple replicas of long-time atomistic MD simulations (3 replicas, 0.3 μs for each replica, 0.9 μs for each octamer system, and 1.8 μs total) of both octamers. We observe that a single GTP cap layer induces structural differences in neighboring PFs, finding that one PF possesses a gradual curvature, compared to the second PF which possesses a kinked conformation. This results in either curling or splaying between these PFs. We suggest that this is due to asymmetric strengths of longitudinal contacts between the two PFs. Furthermore, using NMA, we calculate mechanical properties of these octamer systems and find that octamer system with a single GTP cap layer possesses a lower flexural rigidity.

Systematic Dissociation Pathway Searches Guided by Principal Component Modes

We introduce a novel method, Pathway Search guided by Internal Motions (PSIM), that efficiently finds molecular dissociation pathways of a ligand-receptor system with guidance from principal component (PC) modes obtained from molecular dynamics (MD) simulations. Modeling ligand-receptor dissociation pathways can provide insights into molecular recognition and has practical applications, including understanding kinetic mechanisms and barriers to binding/unbinding as well as design of drugs with desired kinetic properties. PSIM uses PC modes in multilayer internal coordinates to identify natural molecular motions that guide the search for conformational switches and unbinding pathways. The new multilayer internal coordinates overcome problems with Cartesian and classical internal coordinates that fail to smoothly present dihedral rotation or generate nonphysical distortions. We used HIV-1 protease, which has large-scale flap motions, as an example protein to demonstrate use of the multilayer internal coordinates. We provide examples of algorithms and implementation of PSIM with alanine dipeptide and chemical host-guest systems, 2-naphthyl ethanol-β-cyclodextrin and tetramethylammonium-cryptophane complexes. Tetramethylammonium-cryptophane has slow binding/unbinding kinetics. Its residence time, the length to dissociate tetramethylammonium from the host, is ∼14 s from experiments, and PSIM revealed 4 dissociation pathways in approximately 150 CPU h. We also searched the releasing pathways for the product glyceraldehyde-3-phosphate from tryptophan synthase, and one complete dissociation pathway was constructed after running multiple search iterations in approximately 300 CPU h. With guidance by internal PC modes from MD simulations, the PSIM method has advantages over simulation-based methods to search for dissociation pathways of molecular systems with slow noncovalent kinetic behavior.

Equilibrium Denaturation and Preferential Interactions of an RNA Tetraloop with Urea

Urea is an important organic cosolute with implications in maintaining osmotic stress in cells and differentially stabilizing ensembles of folded biomolecules. We report an equilibrium study of urea-induced denaturation of a hyperstable RNA tetraloop through unbiased replica exchange molecular dynamics. We find that, in addition to destabilizing the folded state, urea smooths the RNA free energy landscape by destabilizing specific configurations, and forming favorable interactions with RNA nucleobases. A linear concentration-dependence of the free energy (m-value) is observed, in agreement with the results of other RNA hairpins and proteins. Additionally, analysis of the hydrogen-bonding and stacking interactions within RNA primarily show temperature-dependence, while interactions between RNA and urea primarily show concentration-dependence. Our findings provide valuable insight into the effects of urea on RNA folding and describe the thermodynamics of a basic RNA hairpin as a function of solution chemistry.

Protein Structure Classification and Loop Modeling Using Multiple Ramachandran Distributions

Recently, the study of protein structures using angular representations has attracted much attention among structural biologists. The main challenge is how to efficiently model the continuous conformational space of the protein structures based on the differences and similarities between different Ramachandran plots. Despite the presence of statistical methods for modeling angular data of proteins, there is still a substantial need for more sophisticated and faster statistical tools to model the large-scale circular datasets. To address this need, we have developed a nonparametric method for collective estimation of multiple bivariate density functions for a collection of populations of protein backbone angles. The proposed method takes into account the circular nature of the angular data using trigonometric spline which is more efficient compared to existing methods. This collective density estimation approach is widely applicable when there is a need to estimate multiple density functions from different populations with common features. Moreover, the coefficients of adaptive basis expansion for the fitted densities provide a low-dimensional representation that is useful for visualization, clustering, and classification of the densities. The proposed method provides a novel and unique perspective to two important and challenging problems in protein structure research: structure-based protein classification and angular-sampling-based protein loop structure prediction.

Lid domain plasticity and lipid flexibility modulate enzyme specificity in human monoacylglycerol lipase

Human monoacylglycerol lipase (MAGL) is a membrane-interacting enzyme that generates pro-inflammatory signaling molecules. For this reason, MAGL inhibition is a promising strategy to treat pain, cancer, and neuroinflammatory diseases. MAGL can hydrolyze monoacylglycerols bearing an acyl chain of different lengths and degrees of unsaturation, cleaving primarily the endocannabinoid 2-arachidonoylglycerol. Importantly, the enzymatic binding site of MAGL is confined by a 75-amino-acid-long, flexible cap domain, named 'lid domain', which is structurally similar to that found in several other lipases. However, it is unclear how lid domain plasticity affects catalysis in MAGL. By integrating extensive molecular dynamics simulations and free-energy calculations with mutagenesis and kinetic experiments, we here define a lid-domain-mediated mechanism for substrate selection and binding in MAGL catalysis. In particular, we clarify the key role of Phe159 and Ile179, two conserved residues within the lid domain, in regulating substrate specificity in MAGL. We conclude by proposing that other structurally related lipases may share this lid-domain-mediated mechanism for substrate specificity.

How to understand atomistic molecular dynamics simulations of RNA and protein-RNA complexes?

We provide a critical assessment of explicit-solvent atomistic molecular dynamics (MD) simulations of RNA and protein/RNA complexes, written primarily for non-specialists with an emphasis to explain the limitations of MD. MD simulations can be likened to hypothetical single-molecule experiments starting from single atomistic conformations and investigating genuine thermal sampling of the biomolecules. The main advantage of MD is the unlimited temporal and spatial resolution of positions of all atoms in the simulated systems. Fundamental limitations are the short physical time-scale of simulations, which can be partially alleviated by enhanced-sampling techniques, and the highly approximate atomistic force fields describing the simulated molecules. The applicability and present limitations of MD are demonstrated on studies of tetranucleotides, tetraloops, ribozymes, riboswitches and protein/RNA complexes. Wisely applied simulations respecting the approximations of the model can successfully complement structural and biochemical experiments. WIREs RNA 2017, 8:e1405. doi: 10.1002/wrna.1405 For further resources related to this article, please visit the WIREs website.

Distributions of experimental protein structures on coarse-grained free energy landscapes

Predicting conformational changes of proteins is needed in order to fully comprehend functional mechanisms. With the large number of available structures in sets of related proteins, it is now possible to directly visualize the clusters of conformations and their conformational transitions through the use of principal component analysis. The most striking observation about the distributions of the structures along the principal components is their highly non-uniform distributions. In this work, we use principal component analysis of experimental structures of 50 diverse proteins to extract the most important directions of their motions, sample structures along these directions, and estimate their free energy landscapes by combining knowledge-based potentials and entropy computed from elastic network models. When these resulting motions are visualized upon their coarse-grained free energy landscapes, the basis for conformational pathways becomes readily apparent. Using three well-studied proteins, T4 lysozyme, serum albumin, and sarco-endoplasmic reticular Ca(2+) adenosine triphosphatase (SERCA), as examples, we show that such free energy landscapes of conformational changes provide meaningful insights into the functional dynamics and suggest transition pathways between different conformational states. As a further example, we also show that Monte Carlo simulations on the coarse-grained landscape of HIV-1 protease can directly yield pathways for force-driven conformational changes.

Free-energy landscape of a hyperstable RNA tetraloop

We report the characterization of the energy landscape and the folding/unfolding thermodynamics of a hyperstable RNA tetraloop obtained through high-performance molecular dynamics simulations at microsecond timescales. Sampling of the configurational landscape is conducted using temperature replica exchange molecular dynamics over three isochores at high, ambient, and negative pressures to determine the thermodynamic stability and the free-energy landscape of the tetraloop. The simulations reveal reversible folding/unfolding transitions of the tetraloop into the canonical A-RNA conformation and the presence of two alternative configurations, including a left-handed Z-RNA conformation and a compact purine Triplet. Increasing hydrostatic pressure shows a stabilizing effect on the A-RNA conformation and a destabilization of the left-handed Z-RNA. Our results provide a comprehensive description of the folded free-energy landscape of a hyperstable RNA tetraloop and highlight the significant advances of all-atom molecular dynamics in describing the unbiased folding of a simple RNA secondary structure motif.

Factor Analysis of Conformations and NMR Signals of Rotaxanes: AIMD and Polarizable MD Simulations

The interlocked ⟨rod | ring⟩ structures of pseudorotaxanes and [2]rotaxanes are usually maintained by the complex hydrogen-bonding (H-bonding) network between the rod and ring. Ab initio molecular dynamics (AIMD) using generalized energy-based fragmentation approach and polarizable force field (polar FF)-based molecular dynamics (MD) simulations were performed to investigate the conformational changes of mechanically interlocked systems and to obtain the ensemble-averaged NMR chemical shifts. Factor analysis (FA) demonstrates that the ring H-donor (2,6 pyridinedicarboxamide group) plays an important role in the ring-rod recognition. In comparison to the conventional fixed-charge force field, the polarization effect is crucial to account for the H-bonding interactions in supramolecular systems. In the hybrid scheme, the polar FF-based MD simulations are used to generate different initial states for the AIMD simulations, which are able to give better prediction of ensemble-averaged NMR signals for chemically equivalent amide protons. The magnitude of the deshielding shift of NMR signal is correlated with the length of hydrogen bond. The polar FF model with variable charges shows that the dipole-dipole interactions between the flexible diethylene glycol chain of ring and polar solvents induce the upfield shifts of NMR signals of rod H-donors and the directional distribution of the neighboring CH3CN solvents.

Clustangles: An Open Library for Clustering Angular Data

Dihedral angles are good descriptors of the numerous conformations visited by large, flexible systems, but their analysis requires directional statistics. A single package including the various multivariate statistical methods for angular data that accounts for the distinct topology of such data does not exist. Here, we present a lightweight standalone, operating-system independent package called Clustangles to fill this gap. Clustangles will be useful in analyzing the ever-increasing number of structures in the Protein Data Bank and clustering the copious conformations from increasingly long molecular dynamics simulations.

Principal component analysis of molecular dynamics: on the use of Cartesian vs. internal coordinates

Principal component analysis of molecular dynamics simulations is a popular method to account for the essential dynamics of the system on a low-dimensional free energy landscape. Using Cartesian coordinates, first the translation and overall rotation need to be removed from the trajectory. Since the rotation depends via the moment of inertia on the molecule's structure, this separation is only straightforward for relatively rigid systems. Adopting millisecond molecular dynamics simulations of the folding of villin headpiece and the functional dynamics of BPTI provided by D. E. Shaw Research, it is demonstrated via a comparison of local and global rotational fitting that the structural dynamics of flexible molecules necessarily results in a mixing of overall and internal motion. Even for the small-amplitude functional motion of BPTI, the conformational distribution obtained from a Cartesian principal component analysis therefore reflects to some extend the dominant overall motion rather than the much smaller internal motion of the protein. Internal coordinates such as backbone dihedral angles, on the other hand, are found to yield correct and well-resolved energy landscapes for both examples. The virtues and shortcomings of the choice of various fitting schemes and coordinate sets as well as the generality of these results are discussed in some detail.

Effect of loop composition on the stability and folding kinetics of RNA hairpins with large loops

RNA hairpins are ubiquitous structural elements in biological RNAs, where they have the potential to regulate RNA folding and interactions with other molecules. There are established methods for predicting the thermodynamic stability of an RNA hairpin, but there are still relatively few detailed examinations of the kinetics of folding. Nonetheless, several recent studies indicate that hairpin folding does not proceed via a simple two-state model. Here, we monitor fluorescence from hairpins constructed as molecular beacons in ensemble, fluorescence correlation spectroscopy, and stopped-flow experiments to describe the folding of RNA hairpins with long (15 nucleotide) loops. Our results show that folding of these hairpins occurs through more than two states and that the mechanism of folding includes a fast intermediate phase observed on the tens of microseconds time scale and a slow phase, attributed to formation of the native folded hairpin loop and stem, observed on the milliseconds time scale. The composition of the RNA loop determines the time scale of intermediate and native folded states. Hairpins with a polyuracil loop sequence exhibit slower relaxation of the intermediate state and faster relaxation of the native folded state when compared to that of hairpins with cytosine or adenine in the loop. We hypothesize this composition dependence could be attributed to nucleobase stacking in cytosine and adenine containing regions of the loop, which would be absent in hairpins containing polyuracil loops. Such base stacking could destabilize the intermediate folds, thereby speeding the relaxation of the intermediate relative to similar sized hairpins with no base stacking in the loop. Likewise, the lower intermediate stability could prolong the relaxation of the native folded state.

Computational insights into the inhibition and destabilization of morin on the oligomer of full-length human islet amyloid polypeptide

In this work, we simulated the full-length human islet amyloid peptide (hIAPP) pentamer with and without morins to investigate the mechanism of inhibition and destabilization of this inhibitor on hIAPP oligomer, and identify its possible binding sites on hIAPP.

Molecular mechanism of the inhibition and remodeling of human islet amyloid polypeptide (hIAPP(1-37)) oligomer by resveratrol from molecular dynamics simulation

Natural polyphenols are one of the most actively investigated categories of amyloid inhibitors, and resveratrol has recently been reported to inhibit and remodel the human islet amyloid polypeptide (hIAPP) oligomers and fibrils. However, the exact mechanism of its action is still unknown, especially for the full-length hIAPP1-37. To this end, we performed all-atom molecular dynamics simulations for hIAPP1-37 pentamer with and without resveratrol. The obtained results show that the binding of resveratrol is able to cause remarkable conformational changes of hIAPP1-37 pentamer, in terms of secondary structures, order degree, and morphology. By clustering analysis, two possible binding sites of resveratrol on the hIAPP1-37 pentamer were found, located at the grooves of the top and bottom surfaces of β-sheet layer, respectively. After the binding free energy calculation and residue energy decomposition, it can be concluded that the bottom site is the more possible one, and that the nonpolar interactions act as the driving force for the binding of hIAPP1-37 to resveratrol. In addition, Arg11 is the most important residue for the binding of resveratrol. The full understanding of inhibitory mechanism of resveratrol on the hIAPP1-37 oligomer, and the identification of its binding sites on this protein are helpful for the future design and discovery of new amyloid inhibitors.

Elucidating molecular motion through structural and dynamic filters of energy-minimized conformer ensembles

Complex RNA structures are constructed from helical segments connected by flexible loops that move spontaneously and in response to binding of small molecule ligands and proteins. Understanding the conformational variability of RNA requires the characterization of the coupled time evolution of interconnected flexible domains. To elucidate the collective molecular motions and explore the conformational landscape of the HIV-1 TAR RNA, we describe a new methodology that utilizes energy-minimized structures generated by the program “Fragment Assembly of RNA with Full-Atom Refinement (FARFAR)”. We apply structural filters in the form of experimental residual dipolar couplings (RDCs) to select a subset of discrete energy-minimized conformers and carry out principal component analyses (PCA) to corroborate the choice of the filtered subset. We use this subset of structures to calculate solution T₁ and T_1ρ relaxation times for ¹³C spins in multiple residues in different domains of the molecule using two simulation protocols that we previously published. We match the experimental T₁ times to within 2% and the T_1ρ times to within less than 10% for helical residues. These results introduce a protocol to construct viable dynamic trajectories for RNA molecules that accord well with experimental NMR data and support the notion that the motions of the helical portions of this small RNA can be described by a relatively small number of discrete conformations exchanging over time scales longer than 1 μs.

Improvements in Markov State Model Construction Reveal Many Non-Native Interactions in the Folding of NTL9

Markov State Models (MSMs) provide an automated framework to investigate the dynamical properties of high-dimensional molecular simulations. These models can provide a human-comprehensible picture of the underlying process, and have been successfully used to study protein folding, protein aggregation, protein ligand binding, and other biophysical systems. The MSM requires the construction of a discrete state-space such that two points are in the same state if they can interconvert rapidly. In the following, we suggest an improved method, which utilizes second order Independent Components Analysis (also known as time-structure based Independent Components Analysis, or tICA), to construct the state-space. We apply this method to simulations of NTL9 (provided by Lindorff-Larsen et al. Science2011), and show that the MSM is an improvement over previously built models using conventional distance metrics. Additionally, the resulting model provides insight into the role of non-native contacts by revealing many slow timescales associated with compact, non-native states.

Space warping order parameters and symmetry: application to multiscale simulation of macromolecular assemblies

Coarse-grained features of macromolecular assemblies are understood via a set of order parameters (OPs) constructed in terms of their all-atom configuration. OPs are shown to be slowly changing in time and capture the large-scale spatial features of macromolecular assemblies. The relationship of these variables to the classic notion of OPs based on symmetry breaking phase transitions is discussed. OPs based on space warping transformations are analyzed in detail as they naturally provide a connection between overall structure of an assembly and all-atom configuration. These OPs serve as the basis of a multiscale analysis that yields Langevin equations for OP dynamics. In this context, the characteristics of OPs and PCA modes are compared. The OPs enable efficient all-atom multiscale simulations of the dynamics of macromolecular assemblies in response to changes in microenvironmental conditions, as demonstrated on the structural transitions of cowpea chlorotic mottle virus capsid (CCMV) and RNA of the satellite tobacco mosaic virus (STMV).

The Amber ff99 Force Field Predicts Relative Free Energy Changes for RNA Helix Formation

The ability of the Amber ff99 force field to predict relative free energies of RNA helix formation was investigated. The test systems were three hexaloop RNA hairpins with identical loops and varying stems. The potential of mean force of stretching the hairpins from the native state to an extended conformation was calculated with umbrella sampling. Because the hairpins have identical loop sequence, the differences in free energy changes are only from the stem composition. The Amber ff99 force field was able to correctly predict the order of stabilities of the hairpins, although the magnitude of the free energy change is larger than that determined by optical melting experiments. The two measurements cannot be compared directly because the unfolded state in the optical melting experiments is a random coil, while the end state in the umbrella sampling simulations was an elongated chain. The calculations can be compared to reference data by using a thermodynamic cycle. By applying the thermodynamic cycle to the transitions between the hairpins using simulations and nearest neighbor data, agreement was found to be within the sampling error of simulations, thus demonstrating that ff99 force field is able to accurately predict relative free energies of RNA helix formation.

Hierarchical Order Parameters for Macromolecular Assembly Simulations I: Construction and Dynamical Properties of Order Parameters

Macromolecular assemblies often display a hierarchical organization of macromolecules or their sub-assemblies. To model this, we have formulated a space warping method that enables capturing overall macromolecular structure and dynamics via a set of coarse-grained order parameters (OPs). This article is the first of two describing the construction and computational implementation of an additional class of OPs that has built into them the hierarchical architecture of macromolecular assemblies. To accomplish this, first, the system is divided into subsystems, each of which is described via a representative set of OPs. Then, a global set of variables is constructed from these subsystem-centered OPs to capture overall system organization. Dynamical properties of the resulting OPs are compared to those of our previous nonhierarchical ones, and implied conceptual and computational advantages are discussed for a 100ns, 2 million atom solvated Human Papillomavirus-like particle simulation. In the second article, the hierarchical OPs are shown to enable a multiscale analysis that starts with the N-atom Liouville equation and yields rigorous Langevin equations of stochastic OP dynamics. The latter is demonstrated via a force-field based simulation algorithm that probes key structural transition pathways, simultaneously accounting for all-atom details and overall structure.

Identifying Metastable States of Folding Proteins

Recent molecular dynamics simulations of biopolymers have shown that in many cases the global features of the free energy landscape can be characterized in terms of the metastable conformational states of the system. To identify these states, a conceptionally and computationally simple approach is proposed. It consists of (i) an initial preprocessing via principal component analysis to reduce the dimensionality of the data, followed by k-means clustering to generate up to 10(4) microstates, (ii) the most probable path algorithm to identify the metastable states of the system, and (iii) boundary corrections of these states via the introduction of cluster cores in order to obtain the correct dynamics. By adopting two well-studied model problems, hepta-alanine and the villin headpiece protein, the potential and the performance of the approach are demonstrated.

Order parameters for macromolecules: application to multiscale simulation

Order parameters (OPs) characterizing the nanoscale features of macromolecules are presented. They are generated in a general fashion so that they do not need to be redesigned with each new application. They evolve on time scales much longer than 10(-14) s typical for individual atomic collisions/vibrations. The list of OPs can be automatically increased, and completeness can be determined via a correlation analysis. They serve as the basis of a multiscale analysis that starts with the N-atom Liouville equation and yields rigorous Smoluchowski/Langevin equations of stochastic OP dynamics. Such OPs and the multiscale analysis imply computational algorithms that we demonstrate in an application to ribonucleic acid structural dynamics for 50 ns.

HiRE-RNA: a high resolution coarse-grained energy model for RNA

Although RNAs play many cellular functions, the gap between their sequences and 3D structures is increasing and our knowledge of RNA thermodynamics and long time scale dynamics is still limited at an atomic level of detail. In principle, all-atom molecular dynamics (MD) and replica exchange molecular dynamics (REMD) simulations can investigate these issues, but with current computer facilities, these simulations in explicit solvent have been limited to small RNAs and to short times. To move to larger systems, we can resort to coarse-graining. In this study, we present HiRE-RNA, a generic high resolution coarse-grained model for RNA, and report MD and REMD simulations on two RNAs of 22 and 36 nucleotides. Starting from unfolded structures, the 22-mer folds within 1.8 A rmsd from the NMR structure, while the 36-mer folds within 4.6 A rmsd. Current results suggest that further optimization of the HiRE-RNA force field should open the door to a relevant model for studying large RNAs, such as riboswitches, and for predicting 3D structures from secondary structure information.