Parameters of Monovalent Ions in the AMBER-99 Forcefield: Assessment of Inaccuracies and Proposed Improvements

Can We Ever Develop an Ideal RNA Force Field? Lessons Learned from Simulations of the UUCG RNA Tetraloop and Other Systems

Molecular dynamics (MD) simulations are an important and well-established tool for investigating RNA structural dynamics, but their accuracy relies heavily on the quality of the employed force field (ff). In this work, we present a comprehensive evaluation of widely used pair-additive and polarizable RNA ffs using the challenging UUCG tetraloop (TL) benchmark system. Extensive standard MD simulations, initiated from the NMR structure of the 14-mer UUCG TL, revealed that most ffs did not maintain the native state, instead favoring alternative loop conformations. Notably, three very recent variants of pair-additive ffs, OL3CP-gHBfix21, DES-Amber, and OL3R2.7, successfully preserved the native structure over a 10 × 20 μs time scale. To further assess these ffs, we performed enhanced sampling folding simulations of the shorter 8-mer UUCG TL, starting from the single-stranded conformation. Estimated folding free energies (ΔG°fold) varied significantly among these three ffs, with values of 0.0 ± 0.6, 2.4 ± 0.8, and 7.4 ± 0.2 kcal/mol for OL3CP-gHBfix21, DES-Amber, and OL3R2.7, respectively. The ΔG°fold value predicted by the OL3CP-gHBfix21 ff was closest to experimental estimates, ranging from -1.6 to -0.7 kcal/mol. In contrast, the higher ΔG°fold values obtained using DES-Amber and OL3R2.7 were unexpected, suggesting that key interactions are inaccurately described in the folded, unfolded, or misfolded ensembles. These discrepancies led us to further test DES-Amber and OL3R2.7 ffs on additional RNA and DNA systems, where further performance issues were observed. Our results emphasize the complexity of accurately modeling RNA dynamics and suggest that creating an RNA ff capable of reliably performing across a wide range of RNA systems remains extremely challenging. In conclusion, our study provides valuable insights into the capabilities of current RNA ffs and highlights key areas for future ff development.

An Ultra-Stable, High-Energy and Wide-Temperature-Range Aqueous Alkaline Sodium-Ion Battery with the Microporous C4N/rGO Anode

Common anode materials in aqueous alkaline electrolytes, such as cadmium, metal hydrides and zinc, usually suffer from remarkable biotoxicity, high cost, and serious side reactions. To overcome these problems, we develop a conjugated porous polymer (CPP) in-situ grown on reduced graphene oxide (rGO) and Ketjen black (KB), noted as C4N/rGO and C4N/KB respectively, as the alternative anodes. The results show that C4N/rGO electrode delivers a low redox potential (-0.905 V vs. Ag/AgCl), high specific capacity (268.8 mAh g-1 at 0.2 A g-1), ultra-stable and fast sodium ion storage behavior (216 mAh g-1 at 20 A g-1) in 2 M NaOH electrolyte. The assembled C4N/rGO//Ni(OH)2 full battery can cycle stably more than 38,000 cycles. Furthermore, by adding a small amount of antifreeze additive dimethyl sulfoxide (DMSO) to adjust the hydrogen bonding network, the low-temperature performance of the electrolyte (0.1 DMSO/2 M NaOH) is significantly improved while hydrogen evolution is inhibited. Consequently, the C4N/rGO//Ni(OH)2 full cell exhibits an energy density of 147.3 Wh Kg-1 and ultra-high cycling stability over a wide temperature range from -70 to 45 °C. This work provides an ultra-stable high-capacity CPP-based anode and antifreeze electrolyte for aqueous alkaline batteries and will facilitate their practical applications under extreme conditions.

Rational Design of Targeted Gold Nanoclusters with High Affinity to Integrin αvβ3 for Combination Cancer Therapy

The unique attributes of targeted nano-drug delivery systems (TNDDSs) over conventional cancer therapies in suppressing off-target effects make them one of the most promising options for cancer treatment. There is evidence that the density of surface-conjugated ligands is a crucial factor in achieving the desired therapeutic efficacy of TNDDSs, but this is hardly manageable in conventional nanomaterials. In this context, ligand-protected gold nanoclusters (AuNCs) are excellent candidates for developing new TNDDSs with a unique control on their surface functionalities, thus helping to achieve enhanced delivery performance. Here, we study the interactions and binding free energies between ten different functionalized Au144(SR)60 (SR = thiolate ligand) nanoclusters and integrin αvβ3 using molecular dynamics simulations and the umbrella sampling method to obtain the optimal formulations. The AuNCs were functionalized with anticancer drugs (5-fluorouracil or signaling pathways inhibitors, such as capivasertib, linifanib, tanespimycin, and taselisib) and integrin-targeting peptides (RGD4C or QS13), and we identified the optimal mixed ligand layer to enhance their binding affinity to the cancer cell receptor. The results showed that changing the proportions of the same type of ligands on the surface of AuNCs led to differences of up to 38 kcal/mol in computed binding free energies. RGD4C as the targeting peptide resulted in greater affinity for αvβ3, and in most formulations studied, a higher amount of drug than peptide was needed. Polar and charged residues, such as Ser123, Asp150, Tyr178, Arg214, and Asp251 were found to play a significant role in AuNC binding. Our simulations also revealed that Mn2+ cations are crucial for stabilizing the αvβ3-AuNC complex. These findings demonstrate the potential of carefully designing the surface composition of TNDDSs to optimize their target affinity and specificity.

A Natural Casein-Based Separator with Brick-and-Mortar Structure for Stable, High-Rate Proton Batteries

Rechargeable aqueous proton batteries with small organic molecule anodes are currently considered promising candidates for large-scale energy storage due to their low cost, stable safety, and environmental friendliness. However, the practical application is limited by the poor cycling stability caused by the shuttling of soluble organic molecules between electrodes. Herein, a cell separator is modified by a GO-casein-Cu2+ layer with a brick-and-mortar structure to inhibit the shuttling of small organic molecules. Experimental and calculation results indicate that, attributed to the synergistic effect of physical blocking of casein molecular chains and electrostatic and coordination interactions of Cu2+, bulk dissolution and shuttling of multiple small molecules can be inhibited simultaneously, while H+ transfer across the separators is not almost affected. With the protection of the GO-casein-Cu2+ separator, soluble small molecules, such as diquinoxalino[2,3-a:2',3'-c]phenazine,2,3,8,9,14,15-hexacyano (6CN-DQPZ) exhibit a high reversible capacity of 262.6 mA h g-1 and amazing stability (capacity retention of 92.9% after 1000 cycles at 1 A g-1). In addition, this strategy is also proved available to other active conjugated small molecules, such as indanthrone (IDT), providing a general green sustainable strategy for advancing the use of small organic molecule electrodes in proton cells.

Effect of K+ Force Fields on Ionic Conductivity and Charge Dynamics of KOH in Ethylene Glycol

Predicting ionic conductivity is crucial for developing efficient electrolytes for energy storage and conversion and other electrochemical applications. An accurate estimate of ionic conductivity requires understanding complex ion-ion and ion-solvent interactions governing the charge transport at the molecular level. Molecular simulations can provide key insights into the spatial and temporal behavior of electrolyte constituents. However, such insights depend on the ability of force fields to describe the underlying phenomena. In this work, molecular dynamics simulations were leveraged to delineate the impact of force field parameters on ionic conductivity predictions of potassium hydroxide (KOH) in ethylene glycol (EG). Four different force fields were used to represent the K+ ion. Diffusion-based Nernst-Einstein and correlation-based Einstein approaches were implemented to estimate the ionic conductivity, and the predicted values were compared with experimental measurements. The physical aspects, including ion-aggregation, charge distribution, cluster correlation, and cluster dynamics, were also examined. A force field was identified that provides reasonably accurate Einstein conductivity values and a physically coherent representation of the electrolyte at the molecular level.

Influence of ion and hydration atmospheres on RNA structure and dynamics: insights from advanced theoretical and computational methods

RNA, a highly charged biopolymer composed of negatively charged phosphate groups, defies electrostatic repulsion to adopt well-defined, compact structures. Hence, the presence of positively charged metal ions is crucial not only for RNA's charge neutralization, but they also coherently decorate the ion atmosphere of RNA to stabilize its compact fold. This feature article elucidates various modes of close RNA-ion interactions, with a special emphasis on Mg2+ as an outer-sphere and inner-sphere ion. Through examples, we highlight how inner-sphere chelated Mg2+ stabilizes RNA pseudoknots, while outer-sphere ions can also exert long-range electrostatic interactions, inducing groove narrowing, coaxial helical stacking, and RNA ring formation. In addition to investigating the RNA's ion environment, we note that the RNA's hydration environment is relatively underexplored. Our study delves into its profound interplay with the structural dynamics of RNA, employing state-of-the-art atomistic simulation techniques. Through examples, we illustrate how specific ions and water molecules are associated with RNA functions, leveraging atomistic simulations to identify preferential ion binding and hydration sites. However, understanding their impact(s) on the RNA structure remains challenging due to the involvement of large length and long time scales associated with RNA's dynamic nature. Nevertheless, our contributions and recent advances in coarse-grained simulation techniques offer insights into large-scale structural changes dynamically linked to the RNA ion atmosphere. In this connection, we also review how different cutting-edge computational simulation methods provide a microscopic lens into the influence of ions and hydration on RNA structure and dynamics, elucidating distinct ion atmospheric components and specific hydration layers and their individual and collective impacts.

AMBERff at Scale: Multimillion-Atom Simulations with AMBER Force Fields in NAMD

All-atom molecular dynamics (MD) simulations are an essential structural biology technique with increasing application to multimillion-atom systems, including viruses and cellular machinery. Classical MD simulations rely on parameter sets, such as the AMBER family of force fields (AMBERff), to accurately describe molecular motion. Here, we present an implementation of AMBERff for use in NAMD that overcomes previous limitations to enable high-performance, massively parallel simulations encompassing up to two billion atoms. Single-point potential energy comparisons and case studies on model systems demonstrate that the implementation produces results that are as accurate as running AMBERff in its native engine.

Computer-Aided Reconstruction and Application of Bacillus halodurans S7 Xylanase with Heat and Alkali Resistance

β-1,4-Endoxylanase is the most critical hydrolase for xylan degradation during lignocellulosic biomass utilization. However, its poor stability and activity in hot and alkaline environments hinder its widespread application. In this study, BhS7Xyl from Bacillus halodurans S7 was improved using a computer-aided design through isothermal compressibility (βT) perturbation engineering and by combining three thermostability prediction algorithms (ICPE-TPA). The best variant with remarkable improvement in specific activity, heat resistance (70 °C), and alkaline resistance (both pH 9.0 and 70 °C), R69F/E137M/E145L, exhibited a 4.9-fold increase by wild-type in specific activity (1368.6 U/mg), a 39.4-fold increase in temperature half-life (458.1 min), and a 57.6-fold increase in pH half-life (383.1 min). Furthermore, R69F/E137M/E145L was applied to the hydrolysis of agricultural waste (corncob and hardwood pulp) to efficiently obtain a higher yield of high-value xylooligosaccharides. Overall, the ICPE-TPA strategy has the potential to improve the functional performance of enzymes under extreme conditions for the high-value utilization of lignocellulosic biomass.

Osmotic Force Balance Evaluation of Aqueous Electrolyte Osmotic Pressures and Chemical Potentials

Concentrated aqueous salt solutions are ubiquitous in problems of biological and environmental relevance. The development of accurate force fields that capture the interactions between dissolved species in solution is crucial to simulating these systems to gain molecular insights into the underlying processes under saline conditions. The osmotic pressure is a relatively simple thermodynamic property connecting the experimental and simulation measurements of the associative properties of the ions in solution. Milner [C. Gillespie and S. T. Milner, Soft Matter, 16, 9816 (2020)] proposed a simulation approach to evaluate the osmotic pressures of salts in solution by applying a restraint potential to the ions alone in solution and determining the resulting pressure required to balance that potential, referred to here as the osmotic force balance. Here, we expand Milner's approach, demonstrating that the chemical potentials of the salts in solution as a function of concentration can be fitted to the concentration profiles determined from simulation, additionally providing an analytical expression for the osmotic pressure. This approach is used to determine the osmotic pressures of 15 alkali halide salts in water from simulations. The cross interactions between cations and anions in solution are subsequently optimized to capture their experimental osmotic pressures.

Characterizing the Folding Transition-State Ensembles in the Energy Landscape of an RNA Tetraloop

Molecular dynamics (MD) simulations have become increasingly powerful and can now describe the folding/unfolding of small biomolecules in atomic detail. However, a major challenge in MD simulations is to represent the complex energy landscape of biomolecules using a small number of reaction coordinates. In this study, we investigate the folding pathways of an RNA tetraloop, gcGCAAgc, using five classical MD simulations with a combined simulation time of approximately 120 μs. Our approach involves analyzing the tetraloop dynamics, including the folding transition state ensembles, using the energy landscape visualization method (ELViM). The ELViM is an approach that uses internal distances to compare any two conformations, allowing for a detailed description of the folding process without requiring root mean square alignment of structures. This method has previously been applied to describe the energy landscape of disordered β-amyloid peptides and other proteins. The ELViM results in a non-linear projection of the multidimensional space, providing a comprehensive representation of the tetraloop's energy landscape. Our results reveal four distinct transition-state regions and establish the paths that lead to the folded tetraloop structure. This detailed analysis of the tetraloop's folding process has important implications for understanding RNA folding, and the ELViM approach can be used to study other biomolecules.

Molecular Mechanisms Underlying Detection Sensitivity in Nanoparticle-Assisted NMR Chemosensing

Nanoparticle-assisted nuclear magnetic resonance (NMR) chemosensing exploits monolayer-protected nanoparticles as supramolecular hosts to detect small molecules in complex mixtures via nuclear Overhauser effect experiments with detection limits down to the micromolar range. Still, the structure-sensitivity relationships at the basis of such detection limits are little understood. In this work, we integrate NMR spectroscopy and atomistic molecular dynamics simulations to examine the covariates that affect the sensitivity of different NMR chemosensing experiments [saturation transfer difference (STD), water STD, and high-power water-mediated STD]. Our results show that the intensity of the observed signals correlates with the number and duration of the spin-spin interactions between the analytes and the nanoparticles and/or between the analytes and the nanoparticles' solvation molecules. In turn, these parameters depend on the location and dynamics of each analyte inside the monolayer. This insight will eventually facilitate the tailoring of experimental and computational setups to the analyte's chemistry, making NMR chemosensing an even more effective technique in practical use.

Dissociation Rate Calculation via Constant-Force Steered Molecular Dynamics Simulation

Steered molecular dynamics (SMD) simulations are used to study molecular dissociation events by applying a harmonic force to the molecules and pulling them at a constant velocity. Instead of constant-velocity pulling, we use a constant force: the constant-force SMD (CF-SMD) simulation. The CF-SMD simulation employs a constant force to reduce the activation barrier of molecular dissociation, thereby enhancing the dissociation event. Here, we present the capability of the CF-SMD simulation to estimate the dissociation time at equilibrium. We performed all-atom CF-SMD simulations for NaCl and protein-ligand systems, producing dissociation time at various forces. We extrapolated these values to the dissociation rate without a constant force using Bell's model or the Dudko-Hummer-Szabo model. We demonstrate that the CF-SMD simulations with the models predicted the dissociation time in equilibrium. A CF-SMD simulation is a powerful tool for estimating the dissociation rate in a direct and computationally efficient manner.

Refinement of SARS-CoV-2 envelope protein structure in a native-like environment by molecular dynamics simulations

COVID-19 has become an unprecedented threat to human health. The SARS-CoV-2 envelope (E) protein plays a critical role in the viral maturation process and pathogenesis. Despite intensive investigation, its structure in physiological conditions remains mysterious: no high-resolution full-length structure is available and only an NMR structure of the transmembrane (TM) region has been determined. Here, we present a refined E protein structure, using molecular dynamics (MD) simulations to investigate its structure and dynamics in a 1-palmitoyl-2-oleoyl-sn-glycero-3-phosphocholine (POPC) bilayer system. Our initial homology model based upon the SARS-CoV E protein structure is shown to be unstable in the lipid bilayer, and the H3 helices tend to move away from the membrane center to the membrane-water interface. A more stable model was developed by replacing all H3 helices with the fully equilibrated H3 structure sampled in the MD simulations. This refined model exhibited more favorable contacts with lipids and water than the original homology model and induced local membrane curvature, decreasing local lipid order. Interestingly, the pore radius profiles showed that the channel in both homology and refined models remained in a closed state throughout the simulations. We also demonstrated the utility of this structure to develop anti-SARS-CoV-2 drugs by docking a library of FDA-approved, investigational, and experimental drugs to the refined E protein structure, identifying 20 potential channel blockers. This highlights the power of MD simulations to refine low-resolution structures of membrane proteins in a native-like membrane environment, shedding light on the structural features of the E protein and providing a platform for the development of novel antiviral treatments.

Classical Exchange Polarization: An Anisotropic Variable Polarizability Model

Polarizability, or the tendency of the electron distribution to distort under an electric field, often depends on the local chemical environment. For example, the polarizability of a chloride ion is larger in gas phase compared to a chloride ion solvated in water. This effect is due to the restriction the Pauli exclusion principle places on the allowed electron states. Because no two electrons can occupy the same state, when a highly polarizable atom comes in close contact with other atoms or molecules, the space of allowed states can dramatically decrease. This constraint suggests that an accurate molecular mechanics polarizability model should depend on the radial distance between neighboring atoms. This paper introduces a variable polarizability model within the framework of the HIPPO (Hydrogen-like Intermolecular Polarizable Potential) force field, by damping the polarizability as a function of the orbital overlap of two atoms. This effectively captures the quantum mechanical exchange polarization effects, without explicit utilization of antisymmetrized wave functions. We show that the variable polarizability model remarkably improves the two-body polarization energies and three-body energies of ion-ion and ion-water systems. Under this model, no manual tuning of atomic polarizabilities for monatomic ions is required; the gas-phase polarizability can be used because an appropriate damping function is able to correct the polarizability at short range.

Ionic Dynamics and Vibrational Spectral Diffusion of a Protic Alkylammonium Ionic Salt through Intrinsic Cationic N-H Vibrational Probe from FPMD Simulations

We employed density functional theory (DFT)-based molecular dynamics simulations to explore the structure, dynamics, and spectral properties of the protic ionic entity trimethylammonium chloride (TMACl). Structural investigations include calculating the site-site radial distribution functions (RDFs), the distribution of constituent cations and anions in three-dimensional space, and combined distribution functions of the hydrogen-bonded pair RDF versus angle, revealing the structural characteristics of the ionic solvation and the intermolecular interactions within ions. Further, we determined the instantaneous vibrational stretching frequencies of the intrinsic N-H stretch probe modes by applying the time-series wavelet method. The associated ionic dynamics within the protic ionic compound were investigated by calculating the time-evolution of the fluctuating frequencies and the frequency-time correlation functions (FFCFs). The time scale related to the local structural relaxation process and the average hydrogen bond lifetime, ion cage dynamics, and mean squared displacement were investigated. The faster decay component of the FFCFs, depicting the intermolecular motion of intact hydrogen bonds in TMACl, is 0.07 ps for the Perdew-Burke-Ernzerhof (PBE)-based simulation and 0.06 ps for the PBE-D2 representation. The slower time scale of the longer picosecond decay time component of PBE and PBE-D2 representations are 3.13 and 2.87 ps, respectively. These picosecond time scales represent more significant fluctuations of the hydrogen-bonding partners in the ionic entity and hydrogen-bond jump events accompanied by large angular jumps. The longest picosecond time scales represent structural relaxation, including large angular jumps and ion-pair dynamics. Also, ion cage lifetimes correlate with the slowest time scale of the associated dynamics of vibrational spectral diffusion despite the type of DFT functional. This study benchmarks DFT treatments of the exchange-correlation functional with and without the van der Waals (vdW) dispersion correction scheme. The inclusion of vdW interactions to the PBE functional represents a less structured state of the ionic entity and faster dynamics of the molecular motions relative to the one predicted by the PBE system. All the results illustrate the necessity of accurately describing the Coulomb interactions, vdW dispersive interactive forces, and localized hydrogen bonds required to sustain the energetic balance in this ionic salt.

Diffuse Ions Coordinate Dynamics in a Ribonucleoprotein Assembly

Proper ionic concentrations are required for the functional dynamics of RNA and ribonucleoprotein (RNP) assemblies. While experimental and computational techniques have provided many insights into the properties of chelated ions, less is known about the energetic contributions of diffuse ions to large-scale conformational rearrangements. To address this, we present a model that is designed to quantify the influence of diffuse monovalent and divalent ions on the dynamics of biomolecular assemblies. This model employs all-atom (non-H) resolution and explicit ions, where effective potentials account for hydration effects. We first show that the model accurately predicts the number of excess Mg²⁺ ions for prototypical RNA systems, at a level comparable to modern coarse-grained models. We then apply the model to a complete ribosome and show how the balance between diffuse Mg²⁺ and K⁺ ions can control the dynamics of tRNA molecules during translation. The model predicts differential effects of diffuse ions on the free-energy barrier associated with tRNA entry and the energy of tRNA binding to the ribosome. Together, this analysis reveals the direct impact of diffuse ions on the dynamics of an RNP assembly.

Combined Machine Learning and GRID-Independent Molecular Descriptor (GRIND) Models to Probe the Activity Profiles of 5-Lipoxygenase Activating Protein Inhibitors

Leukotrienes (LTs) are pro-inflammatory lipid mediators derived from arachidonic acid (AA), and their high production has been reported in multiple allergic, autoimmune, and cardiovascular disorders. The biological synthesis of leukotrienes is instigated by transfer of AA to 5-lipoxygenase (5-LO) via the 5-lipoxygenase-activating protein (FLAP). Suppression of FLAP can inhibit LT production at the earliest level, providing relief to patients requiring anti-leukotriene therapy. Over the last 3 decades, several FLAP modulators have been synthesized and pharmacologically tested, but none of them could be able to reach the market. Therefore, it is highly desirable to unveil the structural requirement of FLAP modulators. Here, in this study, supervised machine learning techniques and molecular modeling strategies are adapted to vaticinate the important 2D and 3D anti-inflammatory properties of structurally diverse FLAP inhibitors, respectively. For this purpose, multiple machine learning classification models have been developed to reveal the most relevant 2D features. Furthermore, to probe the 3D molecular basis of interaction of diverse anti-inflammatory compounds with FLAP, molecular docking studies were executed. By using the most probable binding poses from docking studies, the GRIND model was developed, which indicated the positive contribution of four hydrophobic, two hydrogen bond acceptor, and two shape-based features at certain distances from each other towards the inhibitory potency of FLAP modulators. Collectively, this study sheds light on important two-dimensional and three-dimensional structural requirements of FLAP modulators that can potentially guide the development of more potent chemotypes for the treatment of inflammatory disorders.

A Study on the Effect of the Substituent against PAK4 Inhibition Using In Silico Methods

The intrinsic inductive properties of atoms or functional groups depend on the chemical properties of either electron-withdrawing groups (EWGs) or electron-donating groups (EDGs). This study aimed to evaluate in silico methods to determine whether changes in chemical properties of the compound by single atomic substitution affect the biological activity of target proteins and whether the results depend on the properties of the functional groups. We found an imidazo[4,5-b]pyridine-based PAK4 inhibitor, compound 1, as an initial hit compound with the well-defined binding mode for PAK4. In this study, we used both experimental and in silico methods to investigate the effect of atomic substitution on biological activity to optimize the initial hit compound. In biological assays, in the case of EWG, as the size of the halogen atom became smaller and the electronegativity increased, the biological activity IC₅₀ value ranged from 5150 nM to inactive; in the case of EDG, biological activity was inactive. Furthermore, we analyzed the interactions of PAK4 with compounds, focusing on the hinge region residues, L398 and E399, and gatekeeper residues, M395 and K350, of the PAK4 protein using molecular docking studies and fragment molecular orbital (FMO) methods to determine the differences between the effect of EWG and EDG on the activity of target proteins. These results of the docking score and binding energy did not explain the differences in biological activity. However, the pair-interaction energy obtained from the results of the FMO method indicated that there was a difference in the interaction energy between the EWG and EDG in the hinge region residues, L398 and E399, as well as in M395 and K350. The two groups with different properties exhibited opposite electrostatic energy and charge transfer energy between L398 and E399. Additionally, we investigated the electron distribution of the parts interacting with the hinge region by visualizing the molecular electrostatic potential (MEP) surface of the compounds. In conclusion, we described the properties of functional groups that affect biological activity using an in silico method, FMO.

Influence of Monovalent Salts on α-Glycine Crystal Growth from Aqueous Solution: Molecular Dynamics Simulations at Constant Supersaturation Conditions

The growth of α-glycine crystals from aqueous solution is investigated at constant supersaturations by utilizing the constant chemical potential molecular dynamics method. The study considers two faces (010) and (011) that predominantly determine the α-glycine crystal morphology. The general Amber force field (GAFF) with two different charge sets derived from semi-empirical calculations using the complete neglect of differential overlap method (CNDO) and from density functional calculations using the double-numerical plus d- and p-polarization basis set (DNP) is applied to describe α-glycine. The extended simple point charge model is used to simulate water. It is observed that the GAFF/DNP set leads to a much slower integration of glycine molecules into the crystal structure than the GAFF/CNDO set. The GAFF/CNDO set, however, causes the growth even at concentrations well below the experimental solubility. For the GAFF/DNP set, the influence of potassium chloride (KCl) and sodium chloride (NaCl) on the face growth rates is investigated. The parameters recently proposed by Yagasaki et al. [J. Chem. Theory Comput. 2020, 16, 2460-2473] are used to describe salt ions, as standard GAFF parameters lead to the unexpected formation of salt clusters at a concentration lower than the experimental solubility value. According to our simulation results, both salts suppress the growth of the (011) and (010) faces. The inhibiting effect of NaCl is much stronger than that of KCl for the (011) face, while both salts have a similar inhibiting effect on the (010) face. The results are in line with the experimental observations of the impact of salt ions on the α-glycine growth rates for the (011) face reported in literature.

Corrections in the CHARMM36 Parametrization of Chloride Interactions with Proteins, Lipids, and Alkali Cations, and Extension to Other Halide Anions

The nonpolarizable CHARMM force field is one of the most widely used energy functions for all-atom biomolecular simulations. Chloride is the only halide ion included in the latest version, CHARMM36m, and is used widely in simulation studies, often as an electrolyte ion but also as the biological substrate of transport proteins and enzymes. Here, we find that existing parameters systematically underestimate the interaction of Cl^- with proteins and lipids. Accordingly, when examined in solution, little to no Cl^-association can be observed with most components of the protein, including backbone, polar side chains and aromatic rings. The strength of the interaction with cationic side chains and with alkali ions is also incongruent with experimental measurements, specifically osmotic coefficients of concentrated solutions. Consistent with these findings, a 4-μs trajectory of the Cl^--specific transport protein CLC-ec1 shows irreversible Cl^- dissociation from the so-called S_cen binding site, even in a 150 mM NaCl buffer. To correct for these deficiencies, we formulate a series of pair-specific Lennard-Jones parameters that override those resulting from the conventional Lorentz-Berthelot combination rules. These parameters, referred to as NBFIX, are systematically calibrated against available experimental data as well as ab initio geometry optimizations and energy evaluations, for a wide set of binary and ternary Cl^- complexes with protein and lipid analogs and alkali cations. Analogously, we also formulate parameter sets for the other three biological halide ions, namely, fluoride, bromide, and iodide. The resulting parameters are used to calculate the potential of mean force defining the interaction of each anion and each of the protein and lipid analogues in bulk water, revealing association free energies in the range of -0.3 to -3.3 kcal/mol, with the F^- complexes being the least stable. The NBFIX corrections also preserve the Cl^- occupancy of CLC-ec1 in a second 4-μs trajectory. We posit that these optimized molecular-mechanics models provide a more realistic foundation for all-atom simulation studies of processes entailing changes in hydration, recognition, or transport of halide anions.

Remodeling hydrogen bond interactions results in relaxed specificity of Caspase-3

Caspase (or cysteinyl-aspartate specific proteases) enzymes play important roles in apoptosis and inflammation, and the non-identical but overlapping specificity profiles (that is, cleavage recognition sequence) direct cells to different fates. Although all caspases prefer aspartate at the P1 position of the substrate, the caspase-6 subfamily shows preference for valine at the P4 position, while caspase-3 shows preference for aspartate. In comparison with human caspases, caspase-3a from zebrafish has relaxed specificity and demonstrates equal selection for either valine or aspartate at the P4 position. In the context of the caspase-3 conformational landscape, we show that changes in hydrogen bonding near the S3 subsite affect selection of the P4 amino acid. Swapping specificity with caspase-6 requires accessing new conformational space, where each landscape results in optimal binding of DxxD (caspase-3) or VxxD (caspase-6) substrate and simultaneously disfavors binding of the other substrate. Within the context of the caspase-3 conformational landscape, substitutions near the active site result in nearly equal activity against DxxD and VxxD by disrupting a hydrogen bonding network in the substrate binding pocket. The converse substitutions in zebrafish caspase-3a result in increased selection for P4 aspartate over valine. Overall, the data show that the shift in specificity that results in a dual function protease, as in zebrafish caspase-3a, requires fewer amino acid substitutions compared with those required to access new conformational space for swapping substrate specificity, such as between caspases-3 and -6.

Supercomputer-Based Ensemble Docking Drug Discovery Pipeline with Application to Covid-19

We present a supercomputer-driven pipeline for in silico drug discovery using enhanced sampling molecular dynamics (MD) and ensemble docking. Ensemble docking makes use of MD results by docking compound databases into representative protein binding-site conformations, thus taking into account the dynamic properties of the binding sites. We also describe preliminary results obtained for 24 systems involving eight proteins of the proteome of SARS-CoV-2. The MD involves temperature replica exchange enhanced sampling, making use of massively parallel supercomputing to quickly sample the configurational space of protein drug targets. Using the Summit supercomputer at the Oak Ridge National Laboratory, more than 1 ms of enhanced sampling MD can be generated per day. We have ensemble docked repurposing databases to 10 configurations of each of the 24 SARS-CoV-2 systems using AutoDock Vina. Comparison to experiment demonstrates remarkably high hit rates for the top scoring tranches of compounds identified by our ensemble approach. We also demonstrate that, using Autodock-GPU on Summit, it is possible to perform exhaustive docking of one billion compounds in under 24 h. Finally, we discuss preliminary results and planned improvements to the pipeline, including the use of quantum mechanical (QM), machine learning, and artificial intelligence (AI) methods to cluster MD trajectories and rescore docking poses.

Interfacial and bulk properties of concentrated solutions of ammonium nitrate

We conducted molecular dynamics (MD) simulations to calculate the density and surface tension of concentrated ammonium nitrate (AN) solutions up to the solubility limit of ammonium nitrate in water, by combining the SPC/E, SPCE/F and TIP4P/2005 water models with OPLS model for ammonium and nitrate ions. This is the first time that the properties of concentrated solutions of nitrates, especially AN, have been studied by molecular dynamics. We effectively account for the polarisation effects by the electronic continuum correction (ECC), practically realised via rescaling of the ionic charges. We found that, the full-charge force field MD simulations overestimate the experimental results, as the ions experience repulsion from the interface and prefer to remain in the subsurface layer and the bulk solution. In contrast, reducing the ionic charges results in the behaviour that fits well with the experimental data. The nitrate anions display a greater propensity for the interface than the ammonium cations. We accurately predict both the density and the rise in the surface tension of concentrated solutions of AN, recommending TIP4P/2005 for water and the scaled-charge OPLS model (OPLS/ECC) for the ions in the solutions. We observe that, the adsorption of anions to the interface accompanies their depletion in the subsurface layer, which is preferentially occupied by cations, resulting in the formation of the electric double layer. We demonstrate the ion deficiency for up to 3 Å below the surface and establish the requirement to include the polarisability effects in the OPLS model for AN. While these results confirmed the findings of the previous studies for dilute solutions, they are new in the solubility limit. Concentrated solutions exhibit a strong effect of the abundance of solute on the coordination numbers of ions and on the degree of ion pairing. Surprisingly, ion pairing decreases significantly at the interface compared with the bulk. The present study identifies OPLS/ECC, along with TIP4P/2005, to yield accurate predictions of physical properties of concentrated AN, with precision required for industrial applications, such as a formulation of emulsion and fuel-oil explosives that now predominate the civilian use of AN. An application of this model will allow one to predict the surface properties of supersaturated solutions of AN which fall outside the capability of the present laboratory experiments but are important industrially.

Supercomputer-Based Ensemble Docking Drug Discovery Pipeline with Application to Covid-19

We present a supercomputer-driven pipeline for in-silico drug discovery using enhanced sampling molecular dynamics (MD) and ensemble docking. We also describe preliminary results obtained for 23 systems involving eight protein targets of the proteome of SARS CoV-2. THe MD performed is temperature replica-exchange enhanced sampling, making use of the massively parallel supercomputing on the SUMMIT supercomputer at Oak Ridge National Laboratory, with which more than 1ms of enhanced sampling MD can be generated per day. We have ensemble docked repurposing databases to ten configurations of each of the 23 SARS CoV-2 systems using AutoDock Vina. We also demonstrate that using Autodock-GPU on SUMMIT, it is possible to perform exhaustive docking of one billion compounds in under 24 hours. Finally, we discuss preliminary results and planned improvements to the pipeline, including the use of quantum mechanical (QM), machine learning, and AI methods to cluster MD trajectories and rescore docking poses.

Dielectric constant and density of aqueous alkali halide solutions by molecular dynamics: A force field assessment

The concentration dependence of the dielectric constant and the density of 11 aqueous alkali halide solutions (LiCl, NaCl, KCl, RbCl, CsCl, LiI, NaI, KI, CsI, KF, and CsF) is investigated by molecular simulation. Predictions using eight non-polarizable ion force fields combined with the TIP4P/ε water model are compared to experimental data. The influence of the water model and the temperature on the results for the NaCl brine are also addressed. The TIP4P/ε water model improves the accuracy of dielectric constant predictions compared to the SPC/E water model. The solution density is predicted well by most ion models. Almost all ion force fields qualitatively capture the decline of the dielectric constant with the increase of concentration for all solutions and with the increase of temperature for NaCl brine. However, the sampled dielectric constant is mostly in poor quantitative agreement with experimental data. These results are related to the microscopic solution structure, ion pairing, and ultimately the force field parameters. Ion force fields with excessive contact ion pairing and precipitation below the experimental solubility limit generally yield higher dielectric constant values. An adequate reproduction of the experimental solubility limit should therefore be a prerequisite for further investigations of the dielectric constant of aqueous electrolyte solutions by molecular simulation.

Lennard-Jones Parameters Determined to Reproduce the Solubility of NaCl and KCl in SPC/E, TIP3P, and TIP4P/2005 Water

Most classical nonpolarizable ion potential models underestimate the solubility values of NaCl and KCl in water significantly. We determine Lennard-Jones parameters of Na⁺, K⁺, and Cl^- that reproduce the solubility as well as the hydration free energy in dilute aqueous solutions for three water potential models, SPC/E, TIP3P, and TIP4P/2005. The ion-oxygen distance in the solution and the cation-anion distance in salt are also considered in the parametrization. In addition to the target properties, the hydration enthalpy, hydration entropy, self-diffusion coefficient, coordination number, lattice energy, enthalpy of solution, density, viscosity, and number of contact ion pairs are calculated for comparison with 17 frequently used or recently developed ion potential models. The overall performance of each ion model is represented by a global score using a scheme that was originally developed for comparison of water potential models. The global score is better for our models than for the other 17 models not only because of the quite good prediction for the solubility but also because of the relatively small deviation from the experimental value for many of the other properties.

Improved Cation Binding to Lipid Bilayers with Negatively Charged POPS by Effective Inclusion of Electronic Polarization

Phosphatidylserine (PS) lipids are important signaling molecules and the most common negatively charged lipids in eukaryotic membranes. The signaling can be often regulated by calcium, but its interactions with PS headgroups are not fully understood. Classical molecular dynamics (MD) simulations can potentially give detailed description of lipid-ion interactions, but the results strongly depend on the used force field. Here, we apply the electronic continuum correction (ECC) to the Amber Lipid17 parameters of 1-palmitoyl-2-oleoyl-sn-glycero-3-phospho-l-serine (POPS) lipid to improve its interactions with K⁺, Na⁺, and Ca²⁺ ions. The partial charges of the headgroup, glycerol backbone, and carbonyls of POPS, bearing a unit negative charge, were scaled with a factor of 0.75, derived for monovalent ions, and the Lennard-Jones σ parameters of the same segments were scaled with a factor of 0.89. The resulting ECC-POPS model gives more realistic interactions with Na⁺ and Ca²⁺ cations than the original Amber Lipid17 parameters when validated using headgroup order parameters and the "electrometer concept". In ECC-lipids simulations, populations of complexes of Ca²⁺ cations with more than two PS lipids are negligible, and interactions of Ca²⁺ cations with only carboxylate groups are twice more likely than with only phosphate groups, while interactions with carbonyls almost entirely involve other groups as well. Our results pave the way for more realistic MD simulations of biomolecular systems with anionic membranes, allowing signaling processes involving PS and Ca²⁺ to be elucidated.

Solvation Structure and Dynamics of Alkali Metal Halides in an Ionic Liquid from Classical Molecular Dynamics Simulations

The structure and dynamics of the solvation of several alkali metal halides in an ionic liquid (IL), 1-butyl-3-methylimidazolium trifluoromethanesulfonate ([BMIm][OTf]), were investigated by classical molecular dynamics simulations. Various properties such as density, radial distribution functions, coordination numbers, spatial distribution functions, mean-square displacements, self-diffusion coefficients, and velocity-velocity autocorrelation functions were calculated to understand the solvation environment of alkali metal halide salts in IL at various salt concentrations. We observe that the halide anions are coordinated in two different ways with [BMIm]⁺ in all of the mixtures. However, the alkali metal cations interact more with the anion of the ionic liquid as we go from fluorides to iodides. When a common anion was used for the salt and the ionic liquid, we observe significant coordination of Na⁺ with the anion of the ionic liquid in two different ways, which was not observed in the case of lithium salt. We also find that the Li⁺ and Na⁺ ions are involved in the formation of a aggregate-like, stable kinetic entity with anions in their first solvation shells. These aggregate-like entities were seen to be relatively stable, and we noted a rattling motion of salt ions inside them.

Molecular Docking Guided Grid-Independent Descriptor Analysis to Probe the Impact of Water Molecules on Conformational Changes of hERG Inhibitors in Drug Trapping Phenomenon

Human ether a-go-go related gene (hERG) or KV11.1 potassium channels mediate the rapid delayed rectifier current (I_Kr) in cardiac myocytes. Drug-induced inhibition of hERG channels has been implicated in the development of acquired long QT syndrome type (aLQTS) and fatal arrhythmias. Several marketed drugs have been withdrawn for this reason. Therefore, there is considerable interest in developing better tests for predicting drugs which can block the hERG channel. The drug-binding pocket in hERG channels, which lies below the selectivity filter, normally contains K⁺ ions and water molecules. In this study, we test the hypothesis that these water molecules impact drug binding to hERG. We developed 3D QSAR models based on alignment independent descriptors (GRIND) using docked ligands in open and closed conformations of hERG in the presence (solvated) and absence (non-solvated) of water molecules. The ligand–protein interaction fingerprints (PLIF) scheme was used to summarize and compare the interactions. All models delineated similar 3D hERG binding features, however, small deviations of about ~0.4 Å were observed between important hotspots of molecular interaction fields (MIFs) between solvated and non-solvated hERG models. These small changes in conformations do not affect the performance and predictive power of the model to any significant extent. The model that exhibits the best statistical values was attained with a cryo_EM structure of the hERG channel in open state without water. This model also showed the best R² of 0.58 and 0.51 for the internal and external validation test sets respectively. Our results suggest that the inclusion of water molecules during the docking process has little effect on conformations and this conformational change does not impact the predictive ability of the 3D QSAR models.

Single nucleotide polymorphisms in the human ATP7B gene modify the properties of the ATP7B protein

Single nucleotide polymorphisms (SNPs) are the largest source of sequence variation in the human genome. However, their functional significance is not well understood. We show that SNPs in the Wilson disease gene, ATP7B, that produce amino-acid substitutions K832R and R952K, modulate ATP7B properties in vitro and influence serum copper (Cu) status in vivo. The presence of R832 is associated with a lower ATP7B abundance and a diminished trafficking in response to elevated Cu. The K832R substitution alters surface exposure of amino acid residues in the actuator domain and increases its conformational flexibility. All SNP-related ATP7B variants (R832/R952, R832/K952, K832/K952, and K832/R952) have Cu-transport activity. However, the activity of ATP7B-K832/K952 is lower compared to other variants. In humans, the presence of K952 is associated with a higher fraction of exchangeable Cu in serum. Thus, SNPs may modulate the properties of ATP7B and the organism Cu status.

Trimeric coiled coils expand the range of strength, toughness and dynamics of coiled coil motifs under shear

Coiled coils are widespread protein motifs in nature, and promising building blocks for bio-inspired nanomaterials and nanoscale force sensors. Detailed structural insight into their mechanical response is required to understand their role in tissues and to design building blocks for applications. We use all-atom molecular dynamics simulations to elucidate the mechanical response of two types of coiled coils under shear: dimers and trimers. The amino acid sequences of both systems are similar, thus enabling universal (vs. system-specific) features to be identified. The trimer is mechanically more stable - it is both stronger and tougher - than the dimer, withstanding higher forces (127 pN vs. 49 pN at v = 10-3 nm ns-1) and dissipating up to five times more energy before rupture. The deformation mechanism of the trimer at all pull speeds is dominated by progressive helix unfolding. In contrast, at the lowest pull speeds, dimers deform by unfolding/refolding-assisted sliding. The additional helix in the trimer thus both determines the stability of the structure and affects the deformation mechanism, preventing helix sliding. The mechanical response of the coiled coils is not only sensitive to the oligomerization state but also to helix stability: preventing helix unfolding doubles the mechanical strength of the trimer, but decreases its toughness to half. Our results show that coiled coil trimers expand the range of coiled coil responses to an applied shear force. Altering the stability of individual helices against deformation emerges as one possible route towards fine-tuning this response, enabling the use of these motifs as nanomechanical building blocks.

Molecular Dynamics Simulations of Macromolecular Crystals

The structures of biological macromolecules would not be known to their present extent without X-ray crystallography. Most simulations of globular proteins in solution begin by surrounding the crystal structure of the monomer in a bath of water molecules, but the standard simulation employing periodic boundary conditions is already close to a crystal lattice environment. With simple protocols, the same software and molecular models can perform simulations of the crystal lattice, including all asymmetric units and solvent to fill the box. Throughout the history of molecular dynamics, studies of crystal lattices have served to investigate the quality of the underlying force fields, correlate the simulated ensembles to experimental structure factors, and extrapolate the behavior in lattices to behavior in solution. Powerful new computers are enabling molecular simulations with greater realism and statistical convergence. Meanwhile, the advent of exciting new methods in crystallography, including femtosecond free-electron lasers and image reconstruction for time-resolved crystallography on slurries of small crystals, is expanding the range of structures accessible to X-ray diffraction. We review past fusions of simulations and crystallography, then look ahead to the ways that simulations of crystal structures will enhance structural biology in the future.

Enriched Conformational Sampling of DNA and Proteins with a Hybrid Hamiltonian Derived from the Protein Data Bank

In this article, we present a method for the enhanced molecular dynamics simulation of protein and DNA systems called potential of mean force (PMF)-enriched sampling. The method uses partitions derived from the potentials of mean force, which we determined from DNA and protein structures in the Protein Data Bank (PDB). We define a partition function from a set of PDB-derived PMFs, which efficiently compensates for the error introduced by the assumption of a homogeneous partition function from the PDB datasets. The bias based on the PDB-derived partitions is added in the form of a hybrid Hamiltonian using a renormalization method, which adds the PMF-enriched gradient to the system depending on a linear weighting factor and the underlying force field. We validated the method using simulations of dialanine, the folding of TrpCage, and the conformational sampling of the Dickerson⁻Drew DNA dodecamer. Our results show the potential for the PMF-enriched simulation technique to enrich the conformational space of biomolecules along their order parameters, while we also observe a considerable speed increase in the sampling by factors ranging from 13.1 to 82. The novel method can effectively be combined with enhanced sampling or coarse-graining methods to enrich conformational sampling with a partition derived from the PDB.

Concentration-dependent and configuration-dependent interactions of monovalent ions with an RNA tetraloop

Monovalent salt solutions have strongly coupled interactions with biopolymers, from large polyelectrolytes to small RNA oligomers. High salt concentrations have been known to induce transitions in the structure of RNA, producing non-canonical configurations and even driving RNA to precipitate out of solution. Using all-atom molecular dynamics simulations, we model a monovalent salt species (KCL) at high concentrations (0.1-3m) and calculate the equilibrium distributions of water and ions around a small tetraloop-forming RNA oligomer in a variety of structural arrangements: folded A-RNA (canonical) and Z-RNA (non-canonical) tetraloops and unfolded configurations. From these data, we calculate the ion preferential binding coefficients and Donnan coefficients for the RNA oligomer as a function of concentration and structure. We find that cation accumulation is highest around non-canonical Z-RNA configurations at concentrations below 0.5m, while unfolded configurations accumulate the most co-ions in all concentrations. By contrast, canonical A-RNA structures consistently show the lowest accumulations for all ion species. Water distributions vary markedly with RNA configuration but show little dependency on KCL concentration. Based on Donnan coefficient calculations, the net charge of the solution at the surface of the RNA decreases linearly as a function of salt concentration and becomes net-neutral near 2.5-3m KCL for folded configurations, while unfolded configurations still show a positive solution charge. Our findings show that all-atom molecular dynamics can describe the equilibrium distributions of monovalent salt in the presence of small RNA oligomers at KCL concentrations where ion correlation effects become important. Furthermore, these results provide valuable insights into the distributions of water and ions near the RNA oligomer surface as a function of structural configuration.

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.

Modifications to a common phosphorylation network provide individualized control in caspases

Caspase-3 activation and function have been well-defined during programmed cell death, but caspase activity, at low levels, is also required for developmental processes such as lymphoid proliferation and erythroid differentiation. Post-translational modification of caspase-3 is one method used by cells to fine-tune activity below the threshold required for apoptosis, but the allosteric mechanism that reduces activity is unknown. Phosphorylation of caspase-3 at a conserved allosteric site by p38-MAPK (mitogen-activated protein kinase) promotes survival in human neutrophils, and the modification of the loop is thought to be a key regulator in many developmental processes. We utilized phylogenetic, structural, and biophysical studies to define the interaction networks that facilitate the allosteric mechanism in caspase-3. We show that, within the modified loop, Ser¹⁵⁰ evolved with the apoptotic caspases, whereas Thr¹⁵² is a more recent evolutionary event in mammalian caspase-3. Substitutions at Ser¹⁵⁰ result in a pH-dependent decrease in dimer stability, and localized changes in the modified loop propagate to the active site of the same protomer through a connecting surface helix. Likewise, a cluster of hydrophobic amino acids connects the conserved loop to the active site of the second protomer. The presence of Thr¹⁵² in the conserved loop introduces a "kill switch" in mammalian caspase-3, whereas the more ancient Ser¹⁵⁰ reduces without abolishing enzyme activity. These data reveal how evolutionary changes in a conserved allosteric site result in a common pathway for lowering activity during development or a more recent cluster-specific switch to abolish activity.

RNA force field with accuracy comparable to state-of-the-art protein force fields

The complex and often highly dynamic 3D structures of RNA molecules are central to their diverse cellular functions. Molecular dynamics (MD) simulations have played a major role in characterizing the structure and dynamics of proteins, but the physical models (“force fields”) used for simulating nucleic acids are substantially less accurate overall than those used in protein simulations, creating a major challenge for MD studies of RNA. Here, we report an RNA force field capable of describing the structural and thermodynamic properties of RNA molecules with accuracy comparable to state-of-the-art protein force fields. This force field should facilitate the use of MD simulation as a tool for the study of biologically significant RNA molecules and protein–RNA complexes.

Molecular dynamics (MD) simulation has become a powerful tool for characterizing at an atomic level of detail the conformational changes undergone by proteins. The application of such simulations to RNA structures, however, has proven more challenging, due in large part to the fact that the physical models (“force fields”) available for MD simulations of RNA molecules are substantially less accurate in many respects than those currently available for proteins. Here, we introduce an extensive revision of a widely used RNA force field in which the parameters have been modified, based on quantum mechanical calculations and existing experimental information, to more accurately reflect the fundamental forces that stabilize RNA structures. We evaluate these revised parameters through long-timescale MD simulations of a set of RNA molecules that covers a wide range of structural complexity, including single-stranded RNAs, RNA duplexes, RNA hairpins, and riboswitches. The structural and thermodynamic properties measured in these simulations exhibited dramatically improved agreement with experimentally determined values. Based on the comparisons we performed, this RNA force field appears to achieve a level of accuracy comparable to that of state-of-the-art protein force fields, thus significantly advancing the utility of MD simulation as a tool for elucidating the structural dynamics and function of RNA molecules and RNA-containing biological assemblies.

Molecular Dynamics Studies on Liquid/Vapor Interface Properties and Structures of 1-Ethyl-3-methylimidazolium Dimethylphosphate-Water

1-Ethyl-3-methylimidazolium dimethylphosphate ([Emim][Dmp])-water binary solution is one of the promising new working-pairs for absorption heat pump and absorption chillers, which are widely used to recover industrial waste heat. In the absorption process, the mass and heat transfer at the interface greatly depend on interface microscopic structure. Therefore, in order to understand the absorption process, it is very important to study the interface microscopic structure. The liquid-vapor interface properties, as well as the orientation of [Emim]⁺, [Dmp]^-, water at the interface and its aqueous solution with different water mole fraction, were studied using classical all-atom force field by molecular dynamic simulations. The simulated bulk mass density fitted by hyperbolic tangent function for each system was in good agreement with the experiment data, with the relative deviation between simulated and experimental value within 2%. The simulated results indicate that anion is always distributed at the outmost layer of the interface, followed by cation and water molecule. In [Emim][Dmp], the tilt angle of imidazolium rings to the surface normal is in the range of 0° < θ < 12°; for most cation, their ethyl and methyl tilted toward gas phase and bulk, respectively, but for a few cation, their ethyl and the methyl take the opposite orientation. For anion, one methyl prefers to turn toward gas phase and another methyl (PC vector from P atom to C atom) lie nearly parallel to the surface, while one PO vector (from P atom to O atom) turns toward liquid bulk and another PO vector is nearly parallel to the surface. In aqueous solution of [Emim][Dmp], the tilt angle of the imidazolium ring to the surface normal becomes larger (0° < θ < 37°) at the interface, but almost all ethyl intend to tilt toward gas phase and the methyl tilt toward liquid bulk compared with pure [Emim][Dmp]. Two methyl in anion prefer to turn toward gas phase and its two PO vectors toward liquid bulk. This orientation indicates that pure [Emim][Dmp] absorb water in gas phase more easily than [Emim][Dmp]+H₂O system does. Water molecules are distributed in the inner layer of the interface with two OH vectors (from O atom to H atom) tilting toward external surface.

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.