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.

Mechanical Allosteric Couplings of Redox-Induced Conformational Changes in Respiratory Complex I.

We apply linear response theory to calculate mechanical allosteric couplings in respiratory complex I between the iron sulfur cluster N2, located in the catalytic cavity, and the membrane part of the enzyme, separated from it by more than 50 Å. According to our hypothesis, the redox reaction of ubiquinone in the catalytic cavity of the enzyme generates an unbalanced charge that via repulsion of the charged redox center N2 produces local mechanical stress that transmits into the membrane part of the enzyme where it induces proton pumping. Using coarse-grained simulations of the enzyme, we calculated mechanistic allosteric couplings that reveal the pathways of the mechanical transmission of the stress along the enzyme. The results shed light on the recent experimental studies where a stabilization of the enzyme with an introduced disulfide bridge resulted in the abolishing of proton pumping. Simulation of the disulfide bond action indicates a dramatic change of the mechanistic coupling pathways in line with experimental findings.

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.

Quinone binding in respiratory complex I: Going through the eye of a needle. The squeeze-in mechanism of passing the narrow entrance of the quinone site.

At the joint between the membrane and hydrophilic arms of the enzyme, the structure of the respiratory complex I reveals a tunnel-like Q-chamber for ubiquinone binding and reduction. The narrow entrance of the quinone chamber located in ND1 subunit forms a bottleneck (eye of a needle) which in all resolved structures was shown to be too small for a bulky quinone to pass through, and it was suggested that a conformational change is required to open the channel. The closed bottleneck appears to be a well-established feature of all structures reported so-far, both for the so-called open and closed states of the enzyme, with no indication of a stable open state of the bottleneck. We propose a squeeze-in mechanism of the bottleneck passage, where dynamic thermal conformational fluctuations allow quinone to get in and out. Here, using molecular dynamics simulations of the bacterial enzyme, we have identified collective conformational changes that open the quinone chamber bottleneck. The model predicts a significant reduction-due to a need for a rare opening of the bottleneck-of the effective bi-molecular rate constant, in line with the available kinetic data. We discuss possible reasons for such a tight control of the quinone passage into the binding chamber and mechanistic consequences for the quinone two-electron reduction.

Non-Isothermal Decomposition as Efficient and Simple Synthesis Method of NiO/C Nanoparticles for Asymmetric Supercapacitors.

A series of NiO/C nanocomposites with NiO concentrations ranging from 10 to 90 wt% was synthesized using a simple and efficient two-step method based on non-isothermal decomposition of Nickel(II) bis(acetylacetonate). X-ray diffraction (XRD) measurements of these NiO/C nanocomposites demonstrate the presence of ß-NiO. NiO/C nanocomposites are composed of spherical particles distributed over the carbon support surface. The average diameter of nickel oxide spheres increases with the NiO content and are estimated as 36, 50 and 205 nm for nanocomposites with 10, 50 and 80 wt% NiO concentrations, respectively. In turn, each NiO sphere contains several nickel oxide nanoparticles, whose average sizes are 7-8 nm. According to the tests performed using a three-electrode cell, specific capacitance (SC) of NiO/C nanocomposites increases from 200 to 400 F/g as the NiO content achieves a maximum of 60 wt% concentration, after which the SC decreases. The study of the NiO/C composite showing the highest SC in three- and two-electrode cells reveals that its SC remains almost unchanged while increasing the current density, and the sample demonstrates excellent cycling stability properties. Finally, NiO/C (60% NiO) composites are shown to be promising materials for charging quartz clocks with a power rating of 1.5 V (30 min).

Mode Coupling at around M-Point in PZT.

The question of the microscopic origin of the M-superstructure and additional satellite peaks in the Zr-rich lead zirconate-titanate is discussed for nearly 50 years. Clear contradiction between the selection rules of the critical scattering and the superstructure was found preventing unambiguous attributing of the observed superstructure either to the rotation of the oxygen octahedra or to the antiparallel displacements of the lead cations. Detailed analysis of the satellite pattern explained it as the result of the incommensurate phase transition rather than antiphase domains. Critical dynamics is the key point for the formulated problems. Recently, the oxygen tilt soft mode in the PbZr0.976Ti0.024O3 (PZT2.4) was found. But this does not resolve the extinction rules contradiction. The results of the inelastic X-ray scattering study of the phonon spectra of PZT2.4 around M-point are reported. Strong coupling between the lead and oxygen modes resulting in mode anticrossing and creation of the wide flat part in the lowest phonon dispersion curves is identified. This flat part corresponds to the mixture of the displacements of the lead and oxygen ions and can be an explanation of the extinction rules contradiction. Moreover, a flat dispersion surface is a typical prerequisite for the incommensurate phase transition.

Dynamic and Electronic Polarization Corrections to the Dielectric Constant of Water.

The standard approach to calculating the dielectric constant from molecular dynamics (MD) simulations employs a variant of the Kirkwood-Fröhlich methodology. Many popular nonpolarizable models of water, such as TIPnP, give a reasonable agreement with the experimental value of 78. However, it has been argued in the literature that the dipole moments of these models are effective, being smaller than the real dipole of a liquid water molecule by about a factor of [Formula: see text], or roughly [Formula: see text]. If the total or corrected dipole moment is used in calculations, the dielectric constant comes out nearly twice as large, i.e., in the range of 160, which is twice as high as the experimental value. Here we discuss possible reasons for such a discrepancy. One approach takes into account dynamic corrections due to the dependence of the dielectric response of the medium producing the reaction field on the time scale of dipole fluctuations computed in the Kirkwood-Fröhlich method. When dynamic corrections are incorporated into the computational scheme, a much better agreement with the experimental value of the dielectric constant is found when the corrected (real) dipole moment of liquid water is used. However, a formal analysis indicates that the static properties, such as dielectric constant, should not depend on dynamics. We discuss the resulting conundrum and related issues of simulations of electrostatic interactions using periodic boundary conditions in the context of our findings.

On the importance of accounting for nuclear quantum effects in ab initio calibrated force fields in biological simulations.

In many important processes in chemistry, physics, and biology the nuclear degrees of freedom cannot be described using the laws of classical mechanics. At the same time, the vast majority of molecular simulations that employ wide-coverage force fields treat atomic motion classically. In light of the increasing desire for and accelerated development of quantum mechanics (QM)-parameterized interaction models, we reexamine whether the classical treatment is sufficient for a simple but crucial chemical species: alkanes. We show that when using an interaction model or force field in excellent agreement with the "gold standard" QM data, even very basic simulated properties of liquid alkanes, such as densities and heats of vaporization, deviate significantly from experimental values. Inclusion of nuclear quantum effects via techniques that treat nuclear degrees of freedom using the laws of classical mechanics brings the simulated properties much closer to reality.

Prediction of cyclohexane-water distribution coefficient for SAMPL5 drug-like compounds with the QMPFF3 and ARROW polarizable force fields.

We present the performance of blind predictions of water-cyclohexane distribution coefficients for 53 drug-like compounds in the SAMPL5 challenge by three methods currently in use within our group. Two of them utilize QMPFF3 and ARROW, polarizable force-fields of varying complexity, and the third uses the General Amber Force-Field (GAFF). The polarizable FF's are implemented in an in-house MD package, Arbalest. We find that when we had time to parametrize the functional groups with care (batch 0), the polarizable force-fields outperformed the non-polarizable one. Conversely, on the full set of 53 compounds, GAFF performed better than both QMPFF3 and ARROW. We also describe the torsion-restrain method we used to improve sampling of molecular conformational space and thus the overall accuracy of prediction. The SAMPL5 challenge highlighted several drawbacks of our force-fields, such as our significant systematic over-estimation of hydrophobic interactions, specifically for alkanes and aromatic rings.

Docking and migration of carbon monoxide in nitrogenase: the case for gated pockets from infrared spectroscopy and molecular dynamics.

Evidence of a CO docking site near the FeMo cofactor in nitrogenase has been obtained by Fourier transform infrared spectroscopy-monitored low-temperature photolysis. We investigated the possible migration paths for CO from this docking site using molecular dynamics calculations. The simulations support the notion of a gas channel with multiple internal pockets from the active site to the protein exterior. Travel between pockets is gated by the motion of protein residues. Implications for the mechanism of nitrogenase reactions with CO and N2 are discussed.

Polarizable molecular interactions in condensed phase and their equivalent nonpolarizable models.

Earlier, using phenomenological approach, we showed that in some cases polarizable models of condensed phase systems can be reduced to nonpolarizable equivalent models with scaled charges. Examples of such systems include ionic liquids, TIPnP-type models of water, protein force fields, and others, where interactions and dynamics of inherently polarizable species can be accurately described by nonpolarizable models. To describe electrostatic interactions, the effective charges of simple ionic liquids are obtained by scaling the actual charges of ions by a factor of 1/√(ε(el)), which is due to electronic polarization screening effect; the scaling factor of neutral species is more complicated. Here, using several theoretical models, we examine how exactly the scaling factors appear in theory, and how, and under what conditions, polarizable Hamiltonians are reduced to nonpolarizable ones. These models allow one to trace the origin of the scaling factors, determine their values, and obtain important insights on the nature of polarizable interactions in condensed matter systems.

Polarizable Mean-Field Model of Water for Biological Simulations with Amber and Charmm force fields.

Although a great number of computational models of water are available today, the majority of current biological simulations are done with simple models, such as TIP3P and SPC, developed almost thirty years ago and only slightly modified since then. The reason is that the non-polarizable force fields that are mostly used to describe proteins and other biological molecules are incompatible with more sophisticated modern polarizable models of water. The issue is electronic polarizability: in liquid state, in protein, and in vacuum the water molecule is polarized differently, and therefore has different properties; thus the only way to describe all these different media with the same model is to use a polarizable water model. However, to be compatible with the force field of the rest of the system, e.g. a protein, the latter should be polarizable as well. Here we describe a novel model of water that is in effect polarizable, and yet compatible with the standard non-polarizable force fields such as AMBER, CHARMM, GROMOS, OPLS, etc. Thus the model resolves the outstanding problem of incompatibility.

Accounting for electronic polarization in non-polarizable force fields.

The issues of electronic polarizability in molecular dynamics simulations are discussed. We argue that the charges of ionized groups in proteins, and charges of ions in conventional non-polarizable force fields such as CHARMM, AMBER, GROMOS, etc should be scaled by a factor about 0.7. Our model explains why a neglect of electronic solvation energy, which typically amounts to about a half of total solvation energy, in non-polarizable simulations with un-scaled charges can produce a correct result; however, the correct solvation energy of ions does not guarantee the correctness of ion-ion pair interactions in many non-polarizable simulations. The inclusion of electronic screening for charged moieties is shown to result in significant changes in protein dynamics and can give rise to new qualitative results compared with the traditional non-polarizable force field simulations. The model also explains the striking difference between the value of water dipole µâˆ¼ 3D reported in recent ab initio and experimental studies with the value µ(eff)∼ 2.3D typically used in the empirical potentials, such as TIP3P or SPC/E. It is shown that the effective dipole of water can be understood as a scaled value µ(eff) = µ/√ε(el), where ε(el) = 1.78 is the electronic (high-frequency) dielectric constant of water. This simple theoretical framework provides important insights into the nature of the effective parameters, which is crucial when the computational models of liquid water are used for simulations in different environments, such as proteins, or for interaction with solutes.

Computation of hydration free energies of organic solutes with an implicit water model.

A new approach for computing hydration free energies DeltaG(solv) of organic solutes is formulated and parameterized. The method combines a conventional PCM (polarizable continuum model) computation for the electrostatic component DeltaG(el) of DeltaG(solv) and a specially detailed algorithm for treating the complementary nonelectrostatic contributions (DeltaG(nel)). The novel features include the following: (a) two different cavities are used for treating DeltaG(el) and DeltaG(nel). For the latter case the cavity is larger and based on thermal atomic radii (i.e., slightly reduced van der Waals radii). (b) The cavitation component of DeltaG(nel) is taken to be proportional to the volume of the large cavity. (c) In the treatment of van der Waals interactions, all solute atoms are counted explicitly. The corresponding interaction energies are computed as integrals over the surface of the larger cavity; they are based on Lennard Jones (LJ) type potentials for individual solute atoms. The weighting coefficients of these LJ terms are considered as fitting parameters. Testing this method on a collection of 278 uncharged organic solutes gave satisfactory results. The average error (RMSD) between calculated and experimental free energy values varies between 0.15 and 0.5 kcal/mol for different classes of solutes. The larger deviations found for the case of oxygen compounds are probably due to a poor approximation of H-bonding in terms of LJ potentials. For the seven compounds with poorest fit to experiment, the error exceeds 1.5 kcal/mol; these outlier points were not included in the parameterization procedure. Several possible origins of these errors are discussed.