Force Field Benchmark of Amino Acids. 3. Hydration with Scaled Lennard-Jones Interactions

The evolution of the Amber additive protein force field: History, current status, and future

Molecular dynamics simulations are pivotal in elucidating the intricate properties of biological molecules. Nonetheless, the reliability of their outcomes hinges on the precision of the molecular force field utilized. In this perspective, we present a comprehensive review of the developmental trajectory of the Amber additive protein force field, delving into researchers' persistent quest for higher precision force fields and the prevailing challenges. We detail the parameterization process of the Amber protein force fields, emphasizing the specific improvements and retained features in each version compared to their predecessors. Furthermore, we discuss the challenges that current force fields encounter in balancing the interactions of protein-protein, protein-water, and water-water in molecular dynamics simulations, as well as potential solutions to overcome these issues.

Likely Overstabilization of Charge-Charge Interactions in CHARMM36m(w): A Case for a99SB-disp Water

Recent years have witnessed drastic improvements in general-purpose explicit solvent protein force fields, partially driven by the need to study intrinsically disordered proteins (IDPs), and yet the state-of-the-art force fields such as CHARMM36m (c36m) and a99SB-disp still provide different performances in simulating disordered protein states, where c36m has a bias toward overcompaction for large IDPs. Here, we examine the performance of c36m and a99SB-disp in describing the stabilities of a set of 46 amino acid backbone and side chain pairs in various configurations. The free energy results show that c36m systematically predicts stronger interactions compared to a99SB-disp by an average of 0.2 kcal/mol for nonpolar pairs, 0.6 kcal/mol for polar pairs, and 0.8 kcal/mol for salt bridges. The most severe overstabilization in c36m is observed for charged pairs involving the Arg and Glu side chains by up to 2.9 kcal/mol. Importantly, the systematic overstabilization of c36m is only marginally alleviated by c36mw, an ad hoc patch to c36m that increases the dispersion interactions between TIP3P hydrogens and protein atoms. Guided by free energy decomposition, we evaluated if revising the charges alone could alleviate the severe overstabilization of salt bridges of c36m(w) vs a99SB-disp. The results suggested that the direct modification of protein-water interactions is also necessary. Toward this end, we proposed a tentative modification to c36m, referred to as c36mrb-disp, which combines modified Arg side chain charges, retuned backbone hydrogen bonding strength, and the a99SB-disp water model. The modified force field successfully reproduces the secondary structures of several intrinsically disordered peptides and proteins including (AAQAA)3, GB1p, and p53 transactivation domain, while maintaining the stability of a set of folded proteins. This work provides a set of useful systems for benchmarking and optimizing protein force fields and highlights the importance of balancing protein-protein and protein-water electrostatic interactions for accurately describing both folded and disordered proteins.

Scrutinizing the protein hydration shell from molecular dynamics simulations against consensus small-angle scattering data

Biological macromolecules in solution are surrounded by a hydration shell, whose structure differs from the structure of bulk solvent. While the importance of the hydration shell for numerous biological functions is widely acknowledged, it remains unknown how the hydration shell is regulated by macromolecular shape and surface composition, mainly because a quantitative probe of the hydration shell structure has been missing. We show that small-angle scattering in solution using X-rays (SAXS) or neutrons (SANS) provide a protein-specific probe of the protein hydration shell that enables quantitative comparison with molecular simulations. Using explicit-solvent SAXS/SANS predictions, we derived the effect of the hydration shell on the radii of gyration Rg of five proteins using 18 combinations of protein force field and water model. By comparing computed Rg values from SAXS relative to SANS in D2O with consensus SAXS/SANS data from a recent worldwide community effort, we found that several but not all force fields yield a hydration shell contrast in remarkable agreement with experiments. The hydration shell contrast captured by Rg values depends strongly on protein charge and geometric shape, thus providing a protein-specific footprint of protein-water interactions and a novel observable for scrutinizing atomistic hydration shell models against experimental data.

Identification of natural marine compounds as potential inhibitors of CDK2 using molecular docking and molecular dynamics simulation approach

The multifunctional enzyme cyclin-dependent kinase 2 (CDK2) protein is essential for cell proliferation, transcription and modulation of the cell cycle. There is a dysfunction that is connected to various diseases, such as cancer, making it an important treatment target in oncology and beyond. The goal of this study is to identify novel CDK2 ATP binding site inhibitors using in silico drug designing. To find competitive inhibitors for the ATP site, molecular docking, molecular dynamics (MD) simulation and free-binding energy calculations were used. Natural compounds retrieved from marine sources (fungi and algae) were docked against protein, and the best-binding compounds were further evaluated using MD simulations. LIG1, LIG2 and LIG3 (ΔGPB = -19.98, -15.82 and -12.98 kcal/mol, respectively) were placed in the top positions based on their overall binding energy calculated using MMPBSA approach. Stability of the complexes was confirmed by carefully analyzing the rmsd and rmsf patterns retrieved from the MD trajectories. Several residues and areas (Leu124, Val123, Phe80, Leu83, Glu81, Arg 126, Asn132, Leu134, Gln131, Lys88 and Glu195) appear to be critical for inhibitor retention across the active pocket, according to RMSD and RMSF. The dynamics of the ligands inside the active pocket were mapped using principle component analysis. It has been observed that LIG1-3 appear to be the best possible inhibitors due to their high binding energies, interaction pattern and retention inside the active pocket.Communicated by Ramaswamy H. Sarma.

Assessment of an Electrostatic Energy-Based Charge Model for Modeling the Electrostatic Interactions in Water Solvent

The protein force field based on the restrained electrostatic potential (RESP) charges has limitations in accurately describing hydrogen bonding interactions in proteins. To address this issue, we propose an alternative approach called the electrostatic energy-based charges (EEC) model, which shows improved performance in describing electrostatic interactions (EIs) of hydrogen bonds in proteins. In this study, we further investigate the performance of the EEC model in modeling EIs in water solvent. Our findings demonstrate that the fixed EEC model can effectively reproduce the quantum mechanics/molecular mechanics (QM/MM)-calculated EIs between a water molecule and various water solvent environments. However, to achieve the same level of computational accuracy, the electrostatic potential (ESP) charge model needs to fluctuate according to the electrostatic environment. Our analysis indicates that the requirement for charge adjustments depends on the specific mathematical and physical representation of EIs as a function of the environment for deriving charges. By comparing with widely used empirical water models calibrated to reproduce experimental properties, we confirm that the performance of the EEC model in reproducing QM/MM EIs is similar to that of general purpose TIP4P-like water models such as TIP4P-Ew and TIP4P/2005. When comparing the computed 10,000 distinct EI values within the range of -40 to 0 kcal/mol with the QM/MM results calculated at the MP2/aug-cc-pVQZ/TIP3P level, we noticed that the mean unsigned error (MUE) for the EEC model is merely 0.487 kcal/mol, which is remarkably similar to the MUE values of the TIP4P-Ew (0.63 kcal/mol) and TIP4P/2005 (0.579 kcal/mol) models. However, both the RESP method and the TIP3P model exhibit a tendency to overestimate the EIs, as evidenced by their higher MUE values of 1.761 and 1.293 kcal/mol, respectively. EEC-based molecular dynamics simulations have demonstrated that, when combined with appropriate van der Waals parameters, the EEC model can closely reproduce oxygen-oxygen radial distribution function and density of water, showing a remarkable similarity to the well-established TIP4P-like empirical water models. Our results demonstrate that the EEC model has the potential to build force fields with comparable accuracy to more sophisticated empirical TIP4P-like water models.

Solvation thermodynamics from cavity shapes of amino acids

According to common physical chemistry wisdom, the solvent cavities hosting a solute are tightly sewn around it, practically coinciding with its van der Waals surface. Solvation entropy is primarily determined by the surface and the volume of the cavity while enthalpy is determined by the solute-solvent interaction. In this work, we challenge this picture, demonstrating by molecular dynamics simulations that the cavities surrounding the 20 amino acids deviate significantly from the molecular surface. Strikingly, the shape of the cavity alone can be used to predict the solvation free energy, entropy, enthalpy, and hydrophobicity. Solute-solvent interactions involving the different chemical moieties of the amino acid, determine indirectly the cavity shape, and the properties of the branches but do not have to be taken explicitly into account in the prediction model.

Simple and Effective Conformational Sampling Strategy for Intrinsically Disordered Proteins Using the UNRES Web Server

Intrinsically disordered proteins (IDPs) contain more charged amino acids than folded proteins, resulting in a lack of hydrophobic core(s) and a tendency to adopt rapidly interconverting structures rather than well-defined structures. The structural heterogeneity of IDPs, encoded by the amino acid sequence, is closely related to their unique roles in biological pathways, which require them to interact with different binding partners. Recently, Robustelli and co-workers have demonstrated that a balanced all-atom force field can be used to sample heterogeneous structures of disordered proteins ( Proc. Natl. Acad. Sci. U.S.A. 2018, 115, E4758-E4766). However, such a solution requires extensive computational resources, such as Anton supercomputers. Here, we propose a simple and effective solution to sample the conformational space of IDPs using a publicly available web server, namely, the UNited-RESidue (UNRES) web server. Our proposed solution requires no investment in computational resources and no prior knowledge of UNRES. UNRES Replica Exchange Molecular Dynamics (REMD) simulations were carried out on a set of eight disordered proteins at temperatures spanning from 270 to 430 K. Utilizing the latest UNRES force field designed for structured proteins, with proper selections of temperatures, we were able to produce comparable results to all-atom force fields as reported in work done by Robustelli and co-workers. In addition, NMR observables and the radius of gyration calculated from UNRES ensembles were directly compared with the experimental data to further evaluate the accuracy of the UNRES model at all temperatures. Our results suggest that carrying out the UNRES simulations at optimal temperatures using the UNRES web server can be a good alternative to sample heterogeneous structures of IDPs.

An Imbalance in the Force: The Need for Standardized Benchmarks for Molecular Simulation

Force fields (FFs) for molecular simulation have been under development for more than half a century. As with any predictive model, rigorous testing and comparisons of models critically depends on the availability of standardized data sets and benchmarks. While such benchmarks are rather common in the fields of quantum chemistry, this is not the case for empirical FFs. That is, few benchmarks are reused to evaluate FFs, and development teams rather use their own training and test sets. Here we present an overview of currently available tests and benchmarks for computational chemistry, focusing on organic compounds, including halogens and common ions, as FFs for these are the most common ones. We argue that many of the benchmark data sets from quantum chemistry can in fact be reused for evaluating FFs, but new gas phase data is still needed for compounds containing phosphorus and sulfur in different valence states. In addition, more nonequilibrium interaction energies and forces, as well as molecular properties such as electrostatic potentials around compounds, would be beneficial. For the condensed phases there is a large body of experimental data available, and tools to utilize these data in an automated fashion are under development. If FF developers, as well as researchers in artificial intelligence, would adopt a number of these data sets, it would become easier to compare the relative strengths and weaknesses of different models and to, eventually, restore the balance in the force.

Molecular dynamics studies of CED-4/CED-9/EGL-1 ternary complex reveal CED-4 release mechanism in the linear apoptotic pathway of Caenorhabditis elegans

Many steps in programmed cell death are evolutionarily conserved across different species. The Caenorhabditis elegans proteins CED-9, CED-4 and EGL-1 involved in apoptosis are respectively homologous to anti-apoptotic Bcl-2 proteins, Apaf-1 and the "BH3-only" pro-apototic proteins in mammals. In the linear apoptotic pathway of C. elegans, EGL-1 binding to CED-9 leads to the release of CED-4 from CED-9/CED-4 complex. The molecular events leading to this process are not clearly elucidated. While the structures of CED-9 apo, CED-9/EGL-1 and CED-9/CED-4 complexes are known, the CED-9/CED-4/EGL-1 ternary complex structure is not yet determined. In this work, we modeled this ternary complex and performed molecular dynamics simulations of six different systems involving CED-9. CED-9 displays differential dynamics depending upon whether it is bound to CED-4 and/or EGL-1. CED-4 exists as an asymmetric dimer (CED4a and CED4b) in CED-9/CED-4 complex. CED-4a exhibits higher conformational flexibility when simulated without CED-4b. Principal Component Analysis revealed that the direction of CED-4a's winged-helix domain motion differs in the ternary complex. Upon EGL-1 binding, majority of non-covalent interactions involving CARD domain in the CED-4a-CED-9 interface have weakened and only half of the contacts found in the crystal structure between α/β domain of CED4a and CED-9 are found to be stable. Additional stable contacts in the ternary complex and differential dynamics indicate that winged-helix domain may play a key role in CED-4a's dissociation from CED-9. This study has provided a molecular level understanding of potential intermediate states that are likely to occur when CED-4a is released from CED-9.

Spontaneous binding of single-stranded RNAs to RRM proteins visualized by unbiased atomistic simulations with a rescaled RNA force field

Recognition of single-stranded RNA (ssRNA) by RNA recognition motif (RRM) domains is an important class of protein-RNA interactions. Many such complexes were characterized using nuclear magnetic resonance (NMR) and/or X-ray crystallography techniques, revealing ensemble-averaged pictures of the bound states. However, it is becoming widely accepted that better understanding of protein-RNA interactions would be obtained from ensemble descriptions. Indeed, earlier molecular dynamics simulations of bound states indicated visible dynamics at the RNA-RRM interfaces. Here, we report the first atomistic simulation study of spontaneous binding of short RNA sequences to RRM domains of HuR and SRSF1 proteins. Using a millisecond-scale aggregate ensemble of unbiased simulations, we were able to observe a few dozen binding events. HuR RRM3 utilizes a pre-binding state to navigate the RNA sequence to its partially disordered bound state and then to dynamically scan its different binding registers. SRSF1 RRM2 binding is more straightforward but still multiple-pathway. The present study necessitated development of a goal-specific force field modification, scaling down the intramolecular van der Waals interactions of the RNA which also improves description of the RNA-RRM bound state. Our study opens up a new avenue for large-scale atomistic investigations of binding landscapes of protein-RNA complexes, and future perspectives of such research are discussed.

Impacts of targeting different hydration free energy references on the development of ion potentials

Hydration free energy (HFE) as the most important solvation parameter is often targeted in ion model development, even though the reported values differ by dozens of kcal mol^-1 mainly due to the experimentally undetermined HFE of the proton ΔG°(H⁺). The choice of ΔG°(H⁺) obviously affects the hydration of single ions and the relative HFE between the ions with different (magnitude or sign) charges, and the impacts of targeted HFEs on the ion solvation and ion-ion interactions are largely unrevealed. Here we designed point charge models of K⁺, Mg²⁺, Al³⁺, and Cl^- ions targeting a variety of HFE references and then investigated the HFE influences on the simulations of dilute and concentrated ion solutions and of the salt ion pairs in gas, liquid, and solid phases. Targeting one more property of ion-water oxygen distances (IOD) leaves the ion-water binding distance invariant, while the binding strength increases with the decreasing (more negative) HFE of ions as a result of a decrease in ΔG°(H⁺) for the cation and an increase in ΔG°(H⁺) for the anion. The increase in ΔG°(H⁺) leads to strengthened cation-anion interactions and thus to close ion-ion contacts, low osmotic pressures, and small activity derivatives in concentrated ion solutions as well as too stable ion pairs of the salts in different phases. The ion diffusivity and water exchange rates around the ions are simply not HFE dependent but rather more complex. Targeting both the aqueous IOD and salt crystal properties of KCl was also attempted and the comparison between different models indicates the complexity and challenge in obtaining a balanced performance between different phases using classical force fields. Our results also support that a real ΔG°(H⁺) value of -259.8 kcal mol^-1 recommended by Hünenberger and Reif guides ion models to reproduce ion-water and ion-ion interactions reasonably at relatively low salt concentrations. Simulations of a metalloprotein show that a relatively more positive ΔG°(H⁺) for Mg²⁺ model is better for a reasonable description of the metal binding network.

Molecular simulations of complex carbohydrates and glycoconjugates

Complex carbohydrates (glycans) are the most abundant and versatile biopolymers in nature. The broad diversity of biochemical functions that carbohydrates cover is a direct consequence of the variety of 3D architectures they can adopt, displaying branched or linear arrangements, widely ranging in sizes, and with the highest diversity of building blocks of any other natural biopolymer. Despite this unparalleled complexity, a common denominator can be found in the glycans' inherent flexibility, which hinders experimental characterization, but that can be addressed by high-performance computing (HPC)-based molecular simulations. In this short review, I present and discuss the state-of-the-art of molecular simulations of complex carbohydrates and glycoconjugates, highlighting methodological strengths and weaknesses, important insights through emblematic case studies, and suggesting perspectives for future developments.

Atomistic Simulation of Lysozyme in Solutions Crowded by Tetraethylene Glycol: Force Field Dependence

The behavior of biomolecules in crowded environments remains largely unknown due to the accuracy of simulation models and the limited experimental data for comparison. Here we chose a small crowder of tetraethylene glycol (PEG-4) to investigate the self-crowding of PEG-4 solutions and molecular crowding effects on the structure and diffusion of lysozyme at varied concentrations from dilute water to pure PEG-4 liquid. Two Amber-like force fields of Amber14SB and a99SB-disp were examined with TIP3P (fast diffusivity and low viscosity) and a99SB-disp (slow diffusivity and high viscosity) water models, respectively. Compared to the Amber14SB protein simulations, the a99SB-disp model yields more coordinated water and less PEG-4 molecules, less intramolecular hydrogen bonds (HBs), more protein-water HBs, and less protein-PEG HBs as well as stronger interactions and more hydrophilic and less hydrophobic contacts with solvent molecules. The a99SB-disp model offers comparable protein-solvent interactions in concentrated PEG-4 solutions to that in pure water. The PEG-4 crowding leads to a slow-down in the diffusivity of water, PEG-4, and protein, and the decline in the diffusion from atomistic simulations is close to or faster than the hard sphere model that neglects attractive interactions. Despite these differences, the overall structure of lysozyme appears to be maintained well at different PEG-4 concentrations for both force fields, except a slightly large deviation at 370 K at low concentrations with the a99SB-disp model. This is mainly attributed to the strong intramolecular interactions of the protein in the Amber14SB force field and to the large viscosity of the a99SB-disp water model. The results indicate that the protein force fields and the viscosity of crowder solutions affect the simulation of biomolecules under crowding conditions.

Rational Design of Nonbonded Point Charge Models for Monovalent Ions with Lennard-Jones 12-6 Potential

Ions are of central importance in nature, and a variety of potential models was proposed to model ions in different phases for an in-depth exploration of ion-related systems. Here, we developed point charge models of 14 monovalent ions with the traditional 12-6 Lennard-Jones (LJ) potential for use in conjunction with 11 water models of TIP3P, OPC3, SPC/E, SPC/E_b, TIP3P-FB, a99SB-disp, TIP4P-Ew, OPC, TIP4P/2005, TIP4P-D, and TIP4P-FB. The designed models reproduced the real hydration free energy (HFE) of ions and the ion-oxygen distance (IOD) in the first hydration shell accurately and simultaneously, a performance similar to the previously reported 12-6-4 LJ-type ion models (12-6 LJ plus an attractive C₄ term for cations or a repulsive one for anions). This work, along with our previous work on di-, tri-, and tetravalent metal cations (J. Chem. Inf. Model. 2021, 61, 4031-4044; J. Chem. Inf. Model. 2021, 61, 4613-4629), demonstrates the feasibility of the simple 12-6 LJ potential in ion modeling. In order for the 12-6 LJ potential to reproduce both the HFE and IOD, the LJ R parameters need to be close to Shannon's ionic radii for the highly charged cations and to the Stokes's van der Waals (vdW) radii for the monovalent ions. With an additional C₄ term, the R parameters of 12-6-4 LJ ion models agree well with the Stokes's vdW radii and have a more physical meaning. It appears that the C₄ term can be merged into the 12-6 LJ potential by a rational tuning of R and the LJ well depth. Simulations of the osmotic coefficients of alkali chloride solutions and the properties of gaseous and solid alkali halides indicate the necessity of further optimizing ion-ion interactions via, for instance, targeting more properties or using a more physical (polarizable) model.