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.

MoSDeF-GOMC: Python Software for the Creation of Scientific Workflows for the Monte Carlo Simulation Engine GOMC.

MoSDeF-GOMC is a python interface for the Monte Carlo software GOMC to the Molecular Simulation Design Framework (MoSDeF) ecosystem. MoSDeF-GOMC automates the process of generating initial coordinates, assigning force field parameters, and writing coordinate (PDB), connectivity (PSF), force field parameter, and simulation control files. The software lowers entry barriers for novice users while allowing advanced users to create complex workflows that encapsulate simulation setup, execution, and data analysis in a single script. All relevant simulation parameters are encoded within the workflow, ensuring reproducible simulations. MoSDeF-GOMC's capabilities are illustrated through a number of examples, including prediction of the adsorption isotherm for CO2 in IRMOF-1, free energies of hydration for neon and radon over a broad temperature range, and the vapor-liquid coexistence curve of a four-component surrogate for the jet fuel S-8. The MoSDeF-GOMC software is available on GitHub at https://github.com/GOMC-WSU/MoSDeF-GOMC.

Transferable Force Fields from Experimental Scattering Data with Machine Learning Assisted Structure Refinement.

Deriving transferable pair potentials from experimental neutron and X-ray scattering measurements has been a longstanding challenge in condensed matter physics. State-of-the-art scattering analysis techniques estimate real-space microstructure from reciprocal-space total scattering data by refining pair potentials to obtain agreement between simulated and experimental results. Prior attempts to apply these potentials with molecular simulations have revealed inaccurate predictions of thermodynamic fluid properties. In this Letter, a machine learning assisted structure-inversion method applied to neutron scattering patterns of the noble gases (Ne, Ar, Kr, and Xe) is shown to recover transferable pair potentials that accurately reproduce both microstructure and vapor-liquid equilibria from the triple to critical point. Therefore, it is concluded that a single neutron scattering measurement is sufficient to predict macroscopic thermodynamic properties over a wide range of states and provide novel insight into local atomic forces in dense monatomic systems.

py-MCMD: Python Software for Performing Hybrid Monte Carlo/Molecular Dynamics Simulations with GOMC and NAMD.

py-MCMD, an open-source Python software, provides a robust workflow layer that manages communication of relevant system information between the simulation engines NAMD and GOMC and generates coherent thermodynamic properties and trajectories for analysis. To validate the workflow and highlight its capabilities, hybrid Monte Carlo/molecular dynamics (MC/MD) simulations are performed for SPC/E water in the isobaric-isothermal (NPT) and grand canonical (GC) ensembles as well as with Gibbs ensemble Monte Carlo (GEMC). The hybrid MC/MD approach shows close agreement with reference MC simulations and has a computational efficiency that is 2 to 136 times greater than traditional Monte Carlo simulations. MC/MD simulations performed for water in a graphene slit pore illustrate significant gains in sampling efficiency when the coupled-decoupled configurational-bias MC (CD-CBMC) algorithm is used compared with simulations using a single unbiased random trial position. Simulations using CD-CBMC reach equilibrium with 25 times fewer cycles than simulations using a single unbiased random trial position, with a small increase in computational cost. In a more challenging application, hybrid grand canonical Monte Carlo/molecular dynamics (GCMC/MD) simulations are used to hydrate a buried binding pocket in bovine pancreatic trypsin inhibitor. Water occupancies produced by GCMC/MD simulations are in close agreement with crystallographically identified positions, and GCMC/MD simulations have a computational efficiency that is 5 times better than MD simulations. py-MCMD is available on GitHub at https://github.com/GOMC-WSU/py-MCMD.

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.

Self-Assembly and Biogenesis of the Cellular Membrane are Dictated by Membrane Stretch and Composition.

The cell plasma membrane is a highly dynamic organelle governing a wide range of cellular activities including ion transport, secretion, cell division, growth, and development. The fundamental process involved in the addition of new membranes to pre-existing plasma membranes, however, is unclear. Here, we report, using biophysical, morphological, biochemical, and molecular dynamic simulations, the selective incorporation of proteins and lipids from the cytosol into the cell plasma membrane dictated by membrane stretch and composition. Stretching of the cell membrane as a consequence of volume increase following incubation in a hypotonic solution and results in the incorporation of cytosolic proteins and lipids into the existing plasma membrane. Molecular dynamic simulations further confirm that increased membrane stretch results in the rapid insertion of lipids into the existing plasma membrane. Similarly, depletion of cholesterol from the cell plasma membrane selectively alters the incorporation of lipids into the membrane.

Molecular exchange Monte Carlo: A generalized method for identity exchanges in grand canonical Monte Carlo simulations.

A generalized identity exchange algorithm is presented for Monte Carlo simulations in the grand canonical ensemble. The algorithm, referred to as molecular exchange Monte Carlo, may be applied to multicomponent systems of arbitrary molecular topology and provides significant enhancements in the sampling of phase space over a wide range of compositions and temperatures. Three different approaches are presented for the insertion of large molecules, and the pros and cons of each method are discussed. The performance of the algorithms is highlighted through grand canonical Monte Carlo histogram-reweighting simulations performed on a number of systems, which include methane+n-alkanes, butane+perfluorobutane, water+impurity, and 2,2,4-trimethylpentane+neopentane. Relative acceptance efficiencies for molecule transfers of up to 400 times that of standard configurational-bias Monte Carlo are obtained.

Optimized Mie potentials for phase equilibria: Application to noble gases and their mixtures with n-alkanes.

Transferrable force fields, based on n-6 Mie potentials, are presented for noble gases. By tuning the repulsive exponent, ni, it is possible to simultaneously reproduce experimental saturated liquid densities and vapor pressures with high accuracy, from the normal boiling point to the critical point. Vapor-liquid coexistence curves for pure fluids are calculated using histogram reweighting Monte Carlo simulations in the grand canonical ensemble. For all noble gases, saturated liquid densities and vapor pressures are reproduced to within 1% and 4% of experiment, respectively. Radial distribution functions, extracted from NVT and NPT Monte Carlo simulations, are in similarly excellent agreement with experimental data. The transferability of the optimized force fields is assessed through calculations of binary mixture vapor-liquid equilibria. These mixtures include argon + krypton, krypton + xenon, methane + krypton, methane + xenon, krypton + ethane, and xenon + ethane. For all mixtures, excellent agreement with experiment is achieved without the introduction of any binary interaction parameters or multi-body interactions.

Biomolecular simulations with the transferable potentials for phase equilibria: extension to phospholipids.

The Transferable Potentials for Phase Equilibria (TraPPE) is extended to zwitterionic and charged lipids including phosphatidylcholine (PC), phosphatidylethanolamine (PE), phosphatidylserine (PS), and phosphatidylglycerol (PG). The performance of the force field is validated through isothermal-isobaric ensemble (NPT) molecular dynamics simulations of hydrated lipid bilayers performed with the aforementioned head groups combined with saturated and unsaturated alkyl tails containing 12-18 carbon atoms. The effects of water model and sodium ion parameters on the performance of the lipid force field are determined. The predictions of the TraPPE force field for the area per lipid, bilayer thickness, and volume per lipid are within 1-5% of experimental values. Key structural properties of the bilayer, such as order parameter splitting in the sn-2 chain and X-ray form factors, are found to be in close agreement with experimental data.

Prediction of 1-octanol-water and air-water partition coefficients for nitro-aromatic compounds from molecular dynamics simulations.

United-atom force fields, based on the Transferable Potentials for Phase Equilibria (TraPPE), are developed for twelve nitro-aromatic compounds, which include 2,4-dinitrotoluene (2,4-DNT), 2,6-dinitrotoluene (2,6-DNT), 3-nitrotoluene (3-NT), 4-nitrotoluene (4-NT), 1,3-dinitrobenzene (1,3-DNB), 1,4-dinitrobenzene (1,4-DNB), 2,4-dinitroanisole (DNAN), 1,3,5-trinitrobenzene (TNB), 2,4,6-trinitrotoluene (TNT), 2-nitroanisole (2-NAN), 4-nitroanisole (4-NAN) and n-methyl-p-nitroaniline (MNA). 1-Octanol-water and air-water partition coefficients are predicted for the optimized TraPPE-UA force field with adaptive biasing force molecular dynamics simulations, and compared to available experimental data. Log Kow values are predicted with an average absolute deviation of 0.2 log units, while Henry's law constants are predicted to with an average absolute deviation of 0.5 log units. Two additional models are presented for energetic materials with five membered rings for which no experimental data are available in the open literature: 3,5-dinitropyrazole (DNP) and 3-nitro-1,2,4-triazole-5-one (NTO). Investigation of the local microstructure around each solute reveals that 1-octanol is able to form hydrogen bonded chains around the solute, while little organized microstructure was observed around the solutes in water.

Direct calculation of 1-octanol-water partition coefficients from adaptive biasing force molecular dynamics simulations.

The 1-octanol-water partition coefficient log K(ow) of a solute is a key parameter used in the prediction of a wide variety of complex phenomena such as drug availability and bioaccumulation potential of trace contaminants. In this work, adaptive biasing force molecular dynamics simulations are used to determine absolute free energies of hydration, solvation, and 1-octanol-water partition coefficients for n-alkanes from methane to octane. Two approaches are evaluated; the direct transfer of the solute from 1-octanol to water phase, and separate transfers of the solute from the water or 1-octanol phase to vacuum, with both methods yielding statistically indistinguishable results. Calculations performed with the TIP4P and SPC∕E water models and the TraPPE united-atom force field for n-alkanes show that the choice of water model has a negligible effect on predicted free energies of transfer and partition coefficients for n-alkanes. A comparison of calculations using wet and dry octanol phases shows that the predictions for log K(ow) using wet octanol are 0.2-0.4 log units lower than for dry octanol, although this is within the statistical uncertainty of the calculation.

Computational prediction of ionic liquid 1-octanol/water partition coefficients.

Wet 1-octanol/water partition coefficients (log K(ow)) predicted for imidazolium-based ionic liquids using adaptive bias force-molecular dynamics (ABF-MD) simulations lie in excellent agreement with experimental values. These encouraging results suggest prospects for this computational tool in the a priori prediction of log K(ow) values of ionic liquids broadly with possible screening implications as well (e.g., prediction of CO(2)-philic ionic liquids).

Monte Carlo predictions of phase equilibria and structure for dimethyl ether + sulfur dioxide and dimethyl ether + carbon dioxide.

A new force field for dimethyl ether (DME) based on the Lennard-Jones (LJ) 12-6 plus point charge functional form is presented in this work. This force field reproduces experimental saturated liquid and vapor densities, vapor pressures, heats of vaporization, and critical properties to within the statistical uncertainty of the combined experimental and simulation measurements for temperatures between the normal boiling and critical point. Critical parameters and normal boiling point are predicted to within 0.1% of experiment. This force field is used in grand canonical histogram reweighting Monte Carlo simulations to predict the pressure composition diagrams for the binary mixtures DME + SO(2) at 363.15 K and DME + CO(2) at 335.15 and 308.15 K. For the DME + SO(2) mixture, simulation is able to qualitatively reproduce the minimum pressure azeotropy observed experimentally for this mixture, but quantitative errors exist, suggesting that multibody effects may be important in this system. For the DME + CO(2) mixture, simulation is able to predict the pressure-composition behavior within 1% of experimental data. Simulations in the isobaric-isothermal ensemble are used to determine the microstructure of DME + SO(2) and DME + CO(2) mixtures. The DME + SO(2) shows weak pairing between DME and SO(2) molecules, while no specific pairing or aggregation is observed for mixtures of DME + CO(2).

Development of an optimized intermolecular potential for sulfur dioxide.

A new force field for sulfur dioxide, capable of predicting accurately the vapor-liquid equilibria, critical properties, vapor pressure, and heats of vaporization is presented. The new force field reproduces the saturated liquid densities, vapor pressures and heats of vaporization to within 0.5, 2, and 2% of experiment, respectively. The predicted critical properties and the normal boiling point are in excellent agreement with experimental results. Pair distribution functions are calculated for the S-S, S-O, and O-O interactions are in close agreement with neutron and X-ray scattering experiments. In addition to the new force field, similar calculations are performed for four SO(2) intermolecular potentials proposed by Sokolic et al. (Sokolic, F.; Guissani, Y. and Guillot, B. J. Phys. Chem. 1985, 89, 3023], which show that these models work reasonably well near the state point where they were originally parametrized, but large errors in the predicted coexistence properties are displayed at higher and lower temperatures. Comparison of the radial distribution functions show the local structure is only weakly affected by the different force field parameters.

Membrane-directed molecular assembly of the neuronal SNARE complex.

Since the discovery and implication of N-ethylmaleimide-sensitive factor (NSF)-attachment protein receptor (SNARE) proteins in membrane fusion almost two decades ago, there have been significant efforts to understand their involvement at the molecular level. In the current study, we report for the first time the molecular interaction between full-length recombinant t-SNAREs and v-SNARE present in opposing liposomes, leading to the assembly of a t-/v-SNARE ring complex. Using high-resolution electron microscopy, the electron density maps and 3D topography of the membrane-directed SNARE ring complex was determined at nanometre resolution. Similar to the t-/v-SNARE ring complex formed when 50 nm v-SNARE liposomes meet a t-SNARE-reconstituted planer membrane, SNARE rings are also formed when 50 nm diameter isolated synaptic vesicles (SVs) meet a t-SNARE-reconstituted planer lipid membrane. Furthermore, the mathematical prediction of the SNARE ring complex size with reasonable accuracy, and the possible mechanism of membrane-directed t-/v-SNARE ring complex assembly, was determined from the study. Therefore in the present study, using both lipososome-reconstituted recombinant t-/v-SNARE proteins, and native v-SNARE present in isolated SV membrane, the membrane-directed molecular assembly of the neuronal SNARE complex was determined for the first time and its size mathematically predicted. These results provide a new molecular understanding of the universal machinery and mechanism of membrane fusion in cells, having fundamental implications in human health and disease.

Ca(2+) bridging of apposed phospholipid bilayers.

In an effort to provide insight into the mechanism of Ca(2+)-induced fusion of lipid vesicles, molecular dynamics simulations in the isobaric-isothermal ensemble are used to investigate interactions of Ca(2+) with apposed lipid bilayers in close proximity. Simulations reveal the formation of a Ca(2+)-phospholipid "anhydrous complex" between apposed bilayers, whereas similar calculations performed with Na(+) display only complexation between neighboring lipids within the same bilayer. The binding of Ca(2+) to apposed phospholipids brings large regions of the bilayers into close contact (<4 Å), displacing water from phospholipid head groups in the process and creating regions of local dehydration. Dehydration of the apposed bilayers leads to ordering of the phospholipid tails, which is partially disrupted by the presence of Ca(2+)-phospholipid bridges.

Mie potentials for phase equilibria calculations: application to alkanes and perfluoroalkanes.

Transferable united-atom force fields, based on n - 6 Lennard-Jones potentials, are presented for normal alkanes and perfluorocarbons. It is shown that by varying the repulsive exponent the range of the potential can be altered, leading to improved predictions of vapor pressures while also reproducing saturated liquid densities to high accuracy. Histogram-reweighting Monte Carlo simulations in the grand canonical ensemble are used to determine the vapor liquid coexistence curves, vapor pressures, heats of vaporization, and critical points for normal alkanes methane through tetradecane, and perfluorocarbons perfluoromethane through perfluorooctane. For all molecules studied, saturated liquid densities are reproduced to within 1% of experiment. Vapor pressures for normal alkanes and perfluorocarbons were predicted to within 3% and 6% of experiment, respectively. Calculations performed for binary mixture vapor-liquid equilibria for propane + pentane show excellent agreement with experiment, while slight deviations are observed for the ethane + perfluoroethane mixture.

Extension of the transferable potentials for phase equilibria force field to dimethylmethyl phosphonate, sarin, and soman.

The transferable potentials for phase equilibria force field is extended to dimethylmethylphosphonate (DMMP), sarin, and soman by introducing a new interaction site representing the phosphorus atom. Parameters for the phosphorus atom are optimized to reproduce the liquid densities at 303 and 373 K and the normal boiling point of DMMP. Calculations for sarin and soman are performed in predictive mode, without further parameter optimization. Vapor-liquid coexistence curves, critical properties, vapor pressures and heats of vaporization are predicted over a wide range of temperatures with histogram reweighting Monte Carlo simulations in the grand canonical ensemble. Excellent agreement with experiment is achieved for all compounds, with unsigned errors of less than 1% for vapor pressures and normal boiling points and under 5% for heats of vaporization and liquid densities at ambient conditions.

All-atom force field for the prediction of vapor-liquid equilibria and interfacial properties of HFA134a.

A new all-atom force field capable of accurately predicting the bulk and interfacial properties of 1,1,1,2-tetrafluoroethane (HFA134a) is reported. Parameterization of several force fields with different initial charge configurations from ab initio calculations was performed using the histogram reweighting method and Monte Carlo simulations in the grand canonical ensemble. The 12-6 Lennard-Jones well depth and diameter for the different HFA134a models were determined by fitting the simulation results to pure-component vapor-equilibrium data. Initial screening of the force fields was achieved by comparing the calculated and experimental bulk properties. The surface tension of pure HFA134a served as an additional screening property to help discriminate an optimum model. The proposed model reproduces the experimental saturated liquid and vapor densities, and the vapor pressure for HFA134a within average errors of 0.7%, 4.4%, and 3.1%, respectively. Critical density, temperature, vapor pressure, normal boiling point, and heat of vaporization at 298 K are also in good agreement with experimental data with errors of 0.2%, 0.1%, 6.2%, 0%, 2.2%, respectively. The calculated surface tension is found to be within the experimental range of 7.7-8.1 mN.m(-1). The dipole moment of the different models was found to significantly affect the prediction of the vapor pressure and surface tension. The ability of the HFA134a models in predicting the interfacial tension against water is also discussed. The results presented here are relevant in the development of technologies where the more environmentally friendly HFA134a is utilized as a substitute to the ozone depleting chlorofluorocarbon propellants.

Ca2+-dimethylphosphate complex formation: providing insight into Ca2+-mediated local dehydration and membrane fusion in cells.

Earlier studies using X-ray diffraction, light scattering, photon correlation spectroscopy, and atomic force microscopy, strongly suggest that SNARE-induced membrane fusion in cells proceeds as a result of calcium bridging opposing bilayers. The bridging of phospholipid heads groups in the opposing bilayers by calcium leads to the release of water from hydrated Ca(2+) ions as well as the loosely coordinated water at PO-lipid head groups. Local dehydration of phospholipid head groups and the calcium, bridging opposing bilayers, then leads to destabilization of the lipid bilayers and membrane fusion. This hypothesis was tested in the current study by atomistic molecular dynamic simulations in the isobaric-isothermal ensemble using hydrated dimethylphosphate anions (DMP(-)) and calcium cations. Results from the study demonstrate, formation of DMP-Ca(2+) complexes and the consequent removal of water, supporting the hypothesis. Our study further demonstrates that as a result of Ca(2+)-DMP self-assembly, the distance between anionic oxygens between the two DMP molecules is reduced to 2.92A, which is in close agreement with the 2.8A SNARE-induced apposition established between opposing bilayers, reported earlier from X-ray diffraction measurements.