The Determination of Free Energy of Hydration of Water Ions from First Principles.

We model the autoionization of water by determining the free energy of hydration of the major intermediate species of water ions. We represent the smallest ions─the hydroxide ion OH-, the hydronium ion H3O+, and the Zundel ion H5O2+─by bonded models and the more extended ionic structures by strong nonbonded interactions (e.g., the Eigen H9O4+ = H3O+ + 3(H2O) and the Stoyanov H13O6+ = H5O2+ + 4(H2O)). Our models are faithful to the precise QM energies and their components to within 1% or less. Using the calculated free energies and atomization energies, we compute the pKa of pure water from first principles as a consistency check and arrive at a value within 1.3 log units of the experimental one. From these calculations, we conclude that the hydronium ion, and its hydrated state, the Eigen cation, are the dominant species in the water autoionization process.

Combining Force Fields and Neural Networks for an Accurate Representation of Bonded Interactions.

We present a formalism of a neural network encoding bonded interactions in molecules. This intramolecular encoding is consistent with the models of intermolecular interactions previously designed by this group. Variants of the encoding fed into a corresponding neural network may be used to economically improve the representation of torsional degrees of freedom in any force field. We test the accuracy of the reproduction of the ab initio potential energy surface on a set of conformations of two dipeptides, methyl-capped ALA and ASP, in several scenarios. The encoding, either alone or in conjunction with an analytical potential, improves agreement with ab initio energies that are on par with those of other neural network-based potentials. Using the encoding and neural nets in tandem with an analytical model places the agreements firmly within "chemical accuracy" of ±0.5 kcal/mol.

Neural Network Corrections to Intermolecular Interaction Terms of a Molecular Force Field Capture Nuclear Quantum Effects in Calculations of Liquid Thermodynamic Properties.

We incorporate nuclear quantum effects (NQE) in condensed matter simulations by introducing short-range neural network (NN) corrections to the ab initio fitted molecular force field ARROW. Force field NN corrections are fitted to average interaction energies and forces of molecular dimers, which are simulated using the Path Integral Molecular Dynamics (PIMD) technique with restrained centroid positions. The NN-corrected force field allows reproduction of the NQE for computed liquid water and methane properties such as density, radial distribution function (RDF), heat of evaporation (HVAP), and solvation free energy. Accounting for NQE through molecular force field corrections circumvents the need for explicit computationally expensive PIMD simulations in accurate calculations of the properties of chemical and biological systems. The accuracy and locality of pairwise NN NQE corrections indicate that this approach could be applicable to complex heterogeneous systems, such as proteins.

Combining Force Fields and Neural Networks for an Accurate Representation of Chemically Diverse Molecular Interactions.

A key goal of molecular modeling is the accurate reproduction of the true quantum mechanical potential energy of arbitrary molecular ensembles with a tractable classical approximation. The challenges are that analytical expressions found in general purpose force fields struggle to faithfully represent the intermolecular quantum potential energy surface at close distances and in strong interaction regimes; that the more accurate neural network approximations do not capture crucial physics concepts, e.g., nonadditive inductive contributions and application of electric fields; and that the ultra-accurate narrowly targeted models have difficulty generalizing to the entire chemical space. We therefore designed a hybrid wide-coverage intermolecular interaction model consisting of an analytically polarizable force field combined with a short-range neural network correction for the total intermolecular interaction energy. Here, we describe the methodology and apply the model to accurately determine the properties of water, the free energy of solvation of neutral and charged molecules, and the binding free energy of ligands to proteins. The correction is subtyped for distinct chemical species to match the underlying force field, to segment and reduce the amount of quantum training data, and to increase accuracy and computational speed. For the systems considered, the hybrid ab initio parametrized Hamiltonian reproduces the two-body dimer quantum mechanics (QM) energies to within 0.03 kcal/mol and the nonadditive many-molecule contributions to within 2%. Simulations of molecular systems using this interaction model run at speeds of several nanoseconds per day.

Protein-Ligand Binding Free-Energy Calculations with ARROW─A Purely First-Principles Parameterized Polarizable Force Field.

Protein-ligand binding free-energy calculations using molecular dynamics (MD) simulations have emerged as a powerful tool for in silico drug design. Here, we present results obtained with the ARROW force field (FF)─a multipolar polarizable and physics-based model with all parameters fitted entirely to high-level ab initio quantum mechanical (QM) calculations. ARROW has already proven its ability to determine solvation free energy of arbitrary neutral compounds with unprecedented accuracy. The ARROW FF parameterization is now extended to include coverage of all amino acids including charged groups, allowing molecular simulations of a series of protein-ligand systems and prediction of their relative binding free energies. We ensure adequate sampling by applying a novel technique that is based on coupling the Hamiltonian Replica exchange (HREX) with a conformation reservoir generated via potential softening and nonequilibrium MD. ARROW provides predictions with near chemical accuracy (mean absolute error of ∼0.5 kcal/mol) for two of the three protein systems studied here (MCL1 and Thrombin). The third protein system (CDK2) reveals the difficulty in accurately describing dimer interaction energies involving polar and charged species. Overall, for all of the three protein systems studied here, ARROW FF predicts relative binding free energies of ligands with a similar accuracy level as leading nonpolarizable force fields.

Accurate determination of solvation free energies of neutral organic compounds from first principles.

The main goal of molecular simulation is to accurately predict experimental observables of molecular systems. Another long-standing goal is to devise models for arbitrary neutral organic molecules with little or no reliance on experimental data. While separately these goals have been met to various degrees, for an arbitrary system of molecules they have not been achieved simultaneously. For biophysical ensembles that exist at room temperature and pressure, and where the entropic contributions are on par with interaction strengths, it is the free energies that are both most important and most difficult to predict. We compute the free energies of solvation for a diverse set of neutral organic compounds using a polarizable force field fitted entirely to ab initio calculations. The mean absolute errors (MAE) of hydration, cyclohexane solvation, and corresponding partition coefficients are 0.2 kcal/mol, 0.3 kcal/mol and 0.22 log units, i.e. within chemical accuracy. The model (ARROW FF) is multipolar, polarizable, and its accompanying simulation stack includes nuclear quantum effects (NQE). The simulation tools' computational efficiency is on a par with current state-of-the-art packages. The construction of a wide-coverage molecular modelling toolset from first principles, together with its excellent predictive ability in the liquid phase is a major advance in biomolecular simulation.

AMPA Receptor Noncompetitive Inhibitors Occupy a Promiscuous Binding Site.

Noncompetitive inhibitors of AMPA receptors have attracted interest in recent years as antiepileptic drugs. However, their development is hindered by a lack of detailed understanding of the protein-inhibitor interaction mechanisms. Recently, structures of AMPA receptor complexes with the structurally dissimilar, noncompetitive, small-molecule inhibitors pyridone perampanel (PMP), GYKI 53655 (GYKI), and CP 465022 (CP) were resolved, revealing that all three share a common binding site. However, due to the low resolution of the ligands, their exact binding modes and protein-ligand interactions remain ambiguous and insufficiently detailed. We carried out molecular dynamics (MD) simulations on X-ray-resolved and docked AMPA receptor complexes, including thermodynamic integration (TI) to compute ligand binding constants, in order to investigate the inhibitor binding modes in detail and identify key protein-ligand interaction mechanisms. Our analysis and simulations show that the ligand binding pocket at the interface of the receptor's transmembrane domain exhibits features also found in the binding pockets of the multidrug-resistance proteins. The inhibitors bind to such promiscuous pockets by forming multiple weak contacts, while the large, flexible pocket undergoes adjustments to accommodate structurally different ligands in different orientations. TI was able to identify a specific more favorable binding mode for GYKI, while PMP, which has a symmetric ring structure, produced several comparable poses indicating that it may bind in several orientations.

A Hierarchical Approach to Predict Conformation-Dependent Histidine Protonation States in Stable and Flexible Proteins.

Solution acidity measured by pH is an important environmental factor that affects protein structure. It influences the protonation state of protein residues, which in turn may be coupled to protein conformational changes, unfolding, and ligand binding. It remains difficult to compute and measure the p Ka of individual residues, as well as to relate them to pH-dependent protein transitions. This paper presents a hierarchical approach to compute the p Ka of individual protonatable residues, specifically histidines, coupled with underlying structural changes of a protein. A fast and efficient free energy perturbation (FEP) algorithm has also been developed utilizing a fast implementation of standard molecular dynamics (MD) algorithms. Specifically, a CUDA version of the AMBER MD engine is used in this paper. Eight histidine p Ka's are computed in a diverse set of pH stable proteins to demonstrate the proposed approach's utility and assess the predictive quality of the AMBER FF99SB force field. A reference molecule is carefully selected and tested for convergence. A hierarchical approach is used to model p Ka's of the six histidine residues of the diphtheria toxin translocation domain (DTT), which exhibits a diverse ensemble of individual conformations and pH-dependent unfolding. The hierarchical approach consists of first sampling equilibrium conformational ensembles of a protein with protonated and neutral histidine residues via long equilibrium MD simulations (Flores-Canales, J. C.; et al. bioRxiv, 2019, 572040). A clustering method is then used to identify sampled protein conformations, and p Ka's of histidines in each protein conformation are computed. Finally, an ensemble averaging formalism is developed to compute weighted average histidine p Ka's. These can be compared with an apparent experimentally measured p Ka of the DTT protein and thus allows us to propose a mechanism of pH-dependent unfolding of the DTT protein.

Ion Permeation Mechanism in Epithelial Calcium Channel TRVP6.

Calcium is the most abundant metal in the human body that plays vital roles as a cellular electrolyte as well as the smallest and most frequently used signaling molecule. Calcium uptake in epithelial tissues is mediated by tetrameric calcium-selective transient receptor potential (TRP) channels TRPV6 that are implicated in a variety of human diseases, including numerous forms of cancer. We used TRPV6 crystal structures as templates for molecular dynamics simulations to identify ion binding sites and to study the permeation mechanism of calcium and other ions through TRPV6 channels. We found that at low Ca2+ concentrations, a single calcium ion binds at the selectivity filter narrow constriction formed by aspartates D541 and allows Na+ permeation. In the presence of ions, no water binds to or crosses the pore constriction. At high Ca2+ concentrations, calcium permeates the pore according to the knock-off mechanism that includes formation of a short-lived transition state with three calcium ions bound near D541. For Ba2+, the transition state lives longer and the knock-off permeation occurs slower. Gd3+ binds at D541 tightly, blocks the channel and prevents Na+ from permeating the pore. Our results provide structural foundations for understanding permeation and block in tetrameric calcium-selective ion channels.

Molecular mechanisms of bio-catalysis of heme extraction from hemoglobin.

Red blood cell hemolysis in sickle cell disease (SCD) releases free hemoglobin. Extracellular hemoglobin and its degradation products, free heme and iron, are highly toxic due to oxidative stress induction and decrease in nitric oxide availability. We propose an approach that helps to eliminate extracellular hemoglobin toxicity in SCD by employing a bacterial protein system that evolved to extract heme from extracellular hemoglobin. NEAr heme Transporter (NEAT) domains from iron-regulated surface determinant proteins from Staphylococcus aureus specifically bind free heme as well as facilitate its extraction from hemoglobin. We demonstrate that a purified NEAT domain fused with human haptoglobin ß-chain is able to remove heme from hemoglobin and reduce heme content and peroxidase activity of hemoglobin. We further use molecular dynamics (MD) simulations to resolve molecular pathway of heme transfer from hemoglobin to NEAT, and to elucidate molecular mechanism of such heme transferring process. Our study is the first of its kind, in which simulations are employed to characterize the process of heme leaving hemoglobin and subsequent rebinding with a NEAT domain. Our MD results highlight important amino acid residues that facilitate heme transfer and will guide further studies for the selection of best NEAT candidate to attenuate free hemoglobin toxicity.