Performance of 3D-RISM-KH in Predicting Hydration Free Energy: Effect of Solute Parameters

Solvation entropy, enthalpy and free energy prediction using a multi-task deep learning functional in 1D-RISM

Simultaneous calculation of entropies, enthalpies and free energies has been a long-standing challenge in computational chemistry, partly because of the difficulty in obtaining estimates of all three properties from a single consistent simulation methodology. This has been particularly true for methods from the Integral Equation Theory of Molecular Liquids such as the Reference Interaction Site Model which have traditionally given large errors in solvation thermodynamics. Recently, we presented pyRISM-CNN, a combination of the 1 Dimensional Reference Interaction Site Model (1D-RISM) solver, pyRISM, with a deep learning based free energy functional, as a method of predicting solvation free energy (SFE). With this approach, a 40-fold improvement in prediction accuracy was delivered for a multi-solvent, multi-temperature dataset when compared to the standard 1D-RISM theory [Fowles et al., Digital Discovery, 2023, 2, 177-188]. Here, we report three further developments to the pyRISM-CNN methodology. Firstly, solvation free energies have been introduced for organic molecular ions in methanol or water solvent systems at 298 K, with errors below 4 kcal mol^-1 obtained without the need for corrections or additional descriptors. Secondly, the number of solvents in the training data has been expanded from carbon tetrachloride, water and chloroform to now also include methanol. For neutral solutes, prediction errors nearing or below 1 kcal mol^-1 are obtained for each organic solvent system at 298 K and water solvent systems at 273-373 K. Lastly, pyRISM-CNN was successfully applied to the simultaneous prediction of solvation enthalpy, entropy and free energy through a multi-task learning approach, with errors of 1.04, 0.98 and 0.47 kcal mol^-1, respectively, for water solvent systems at 298 K.

Identifying Systematic Force Field Errors Using a 3D-RISM Element Counting Correction

Hydration free energies of small molecules are commonly used as benchmarks for solvation models. However, errors in predicting hydration free energies are partially due to the force fields used and not just the solvation model. To address this, we have used the 3D reference interaction site model (3D-RISM) of molecular solvation and existing benchmark explicit solvent calculations with a simple element count correction (ECC) to identify problems with the non-bond parameters in the general AMBER force field (GAFF). 3D-RISM was used to calculate hydration free energies of all 642 molecules in the FreeSolv database, and a partial molar volume correction (PMVC), ECC, and their combination (PMVECC) were applied to the results. The PMVECC produced a mean unsigned error of 1.01±0.04kcal/mol and root mean squared error of 1.44±0.07kcal/mol, better than the benchmark explicit solvent calculations from FreeSolv, and required less than 15 s of computing time per molecule on a single CPU core. Importantly, parameters for PMVECC showed systematic errors for molecules containing Cl, Br, I, and P. Applying ECC to the explicit solvent hydration free energies found the same systematic errors. The results strongly suggest that some small adjustments to the Lennard-Jones parameters for GAFF will lead to improved hydration free energy calculations for all solvent models.

Explainable Solvation Free Energy Prediction Combining Graph Neural Networks with Chemical Intuition

The prediction of a molecule's solvation Gibbs free (ΔG_solv) energy in a given solvent is an important task which has traditionally been carried out via quantum chemical continuum methods or force field-based molecular simulations. Machine learning (ML) and graph neural networks in particular have emerged as powerful techniques for elucidating structure-property relationships. This work presents a graph neural network (GNN) for the prediction of ΔG_solv which, in addition to encoding typical atom and bond-level features, incorporates chemically intuitive, solvation-relevant parameters into the featurization process: semiempirical partial atomic charges and solvent dielectric constant. Solute-solvent interactions are included via an interaction map layer which can be visualized to examine solubility-enhancing or -decreasing interactions learnt by the model. On a test set of small organic molecules, our GNN predicts ΔG_solv in water and cyclohexane with an accuracy comparable to polarizable and ab initio generated force field methods [mean absolute error (MAE) = 0.4 and 0.2 kcal mol^-1, respectively], without the need for any molecular simulation. For the FreeSolv data set of hydration free energies, the test MAE is 0.7 kcal mol^-1. Interpretability and applicability of the model is highlighted through several examples including rationalizing the increased solubility of modified diaminoanthraquinones in organic solvents. The clear explanations afforded by our GNN allow for easy understanding of the model's predictions, giving the experimental chemist confidence in employing ML models toward more optimized synthetic routes.

Predicting Solvent Effects on SN2 Reaction Rates: Comparison of QM/MM, Implicit, and MM Explicit Solvent Models

Solvents are one of the key variables in the optimization of a synthesis yield or properties of a synthesis product. In this paper, contemporary solvent models are applied to predict the rates of S_N2 reactions in a range of aqueous and non-aqueous solvents. High-level CCSD(T)/CBS//M06-2X/6-31+G(d) gas phase energies were combined with solvation free energies from SMD, SM12, and ADF-COSMO-RS continuum solvent models, as well as molecular mechanics (MM) explicit solvent models with different atomic charge schemes to predict the rate constants of three S_N2 reactions in eight protic and aprotic solvents. It is revealed that the prediction of rate constants in organic solvents is not necessarily less challenging than in water and popular solvent models struggle to predict their rate constants to within 3 log units of experimental values. Among the continuum solvent models, the ADF-COSMO-RS model performed the best in predicting absolute rate contants while the SM12 model was best at predicting relative rate constants with an average accuracy of about 1.5 and 0.8 log units, respectively. The use of computationally more demanding MM explicit solvent models did not translate to improvements in absolute rate constants but was quite effective at predicting relative rate constants due to systematic error cancellation. Free energy barriers obtained from umbrella sampling with explicit solvent QM/MM simulations led to excellent agreement with experimental values, provided that a validated level of theory is used to treat the QM region.

Ab initio machine learning of phase space averages

Equilibrium structures determine material properties and biochemical functions. We here propose to machine learn phase space averages, conventionally obtained by ab initio or force-field-based molecular dynamics (MD) or Monte Carlo (MC) simulations. In analogy to ab initio MD, our ab initio machine learning (AIML) model does not require bond topologies and, therefore, enables a general machine learning pathway to obtain ensemble properties throughout the chemical compound space. We demonstrate AIML for predicting Boltzmann averaged structures after training on hundreds of MD trajectories. The AIML output is subsequently used to train machine learning models of free energies of solvation using experimental data and to reach competitive prediction errors (mean absolute error ∼ 0.8 kcal/mol) for out-of-sample molecules-within milliseconds. As such, AIML effectively bypasses the need for MD or MC-based phase space sampling, enabling exploration campaigns of Boltzmann averages throughout the chemical compound space at a much accelerated pace. We contextualize our findings by comparison to state-of-the-art methods resulting in a Pareto plot for the free energy of solvation predictions in terms of accuracy and time.

Implementation and Optimization of the Embedded Cluster Reference Interaction Site Model with Atomic Charges

In this work, we implemented the embedded cluster reference interaction site model (EC-RISM) originally developed by Kloss, Heil, and Kast (J. Phys. Chem. B2008, 112, 4337–4343). This method combines quantum mechanical calculations with the 3D reference interaction site model (3D-RISM). Numerous options, such as buffer, grid space, basis set, charge model, water model, closure relation, and so forth, were investigated to find the best settings. Additionally, the small point charges, which are derived from the solvent distribution from the 3D-RISM solution to represent the solvent in the QM calculation, were neglected to reduce the overhead without the loss of accuracy. On the MNSOL[a], MNSOL, and FreeSolv databases, our implemented and optimized method provides solvation free energies in water with 5.70, 6.32, and 6.44 kJ/mol root-mean-square deviations, respectively, but with different settings, 5.22, 6.08, and 6.63 kJ/mol can also be achieved. Only solvent models containing fitting parameters, like COSMO-RS and EC-RISM with universal correction and directly used electrostatic potential, perform better than our EC-RISM implementation with atomic charges.

Solvation energetics of proteins and their aggregates analyzed by all-atom molecular dynamics simulations and the energy-representation theory of solvation

Solvation is a controlling factor for the structure and function of proteins. This article addresses the effects of solvation from an energetic perspective for the fluctuations and cosolvent-induced changes in protein structures and the equilibrium of aggregate formation for a peptide. A theoretical framework to analyze the solvation effects with an explicit solvent is introduced by adopting the energy-representation theory of solvation, and the connection of the solvation free energy to the protein structure and the aggregation tendency is quantitatively described in combination with all-atom molecular dynamics simulations. The interaction components that govern the solvation effects on the structural variations of proteins are further identified through correlation analysis, and a computational scheme to assess the shift of an aggregation equilibrium due to the addition of a cosolvent is provided.

Water Dissolved in a Variety of Polymers Studied by Molecular Dynamics Simulation and a Theory of Solutions

The performance of a polymer medium as a separation membrane is determined by the dissolution free energy ΔG and diffusion coefficient D of the permeant. In this work, ΔG and D of water are investigated with all-atom molecular dynamics simulation in a wide variety of polymer species in the amorphous state. The computed ΔG is shown to agree well with the experimental value for linear homopolymers, and the degrees of polymerization of the homopolymers do not affect ΔG when they are beyond ∼10. The copolymers of ethylene-vinylidene difluoride, ethylene-vinyl acetate, and ethylene-acrylamide are then examined by changing the repeating patterns of the constituent monomers in both the periodic and graft forms. It is found that ΔG is determined primarily by the overall compositions of the monomers and is not affected by the copolymerization topology (periodic or graft). The hydrophobicity of the copolymer is enhanced, furthermore, when the hydrophobicity and hydrophilicity of the ethylene and non-ethylene parts are well contrasted and those parts are fragmented along the polymer chain. According to the computed D, the diffusivity of water tends to be larger when the (co)polymer is more hydrophobic and ΔG is more positive. D is actually seen to vary by orders of magnitude with the polymer structures, while the effect of the polymer species on the water permeation is stronger for ΔG than for D.

Accurate prediction of hydration free energies and solvation structures using molecular density functional theory with a simple bridge functional

This paper assesses the ability of molecular density functional theory to predict efficiently and accurately the hydration free energies of molecular solutes and the surrounding microscopic water structure. A wide range of solutes were investigated, including hydrophobes, water as a solute, and the FreeSolv database containing 642 drug-like molecules having a variety of shapes and sizes. The usual second-order approximation of the theory is corrected by a third-order, angular-independent bridge functional. The overall functional is parameter-free in the sense that the only inputs are bulk water properties, independent of the solutes considered. These inputs are the direct correlation function, compressibility, liquid-gas surface tension, and excess chemical potential of the solvent. Compared to molecular simulations with the same force field and the same fixed solute geometries, the present theory is shown to describe accurately the solvation free energy and structure of both hydrophobic and hydrophilic solutes. Overall, the method yields a precision of order 0.5 k_BT for the hydration free energies of the FreeSolv database, with a computer speedup of 3 orders of magnitude. The theory remains to be improved for a better description of the H-bonding structure and the hydration free energy of charged solutes.

Biomolecular Simulations with the Three-Dimensional Reference Interaction Site Model with the Kovalenko-Hirata Closure Molecular Solvation Theory

The statistical mechanics-based 3-dimensional reference interaction site model with the Kovalenko-Hirata closure (3D-RISM-KH) molecular solvation theory has proven to be an essential part of a multiscale modeling framework, covering a vast region of molecular simulation techniques. The successful application ranges from the small molecule solvation energy to the bulk phase behavior of polymers, macromolecules, etc. The 3D-RISM-KH successfully predicts and explains the molecular mechanisms of self-assembly and aggregation of proteins and peptides related to neurodegeneration, protein-ligand binding, and structure-function related solvation properties. Upon coupling the 3D-RISM-KH theory with a novel multiple time-step molecular dynamic (MD) of the solute biomolecule stabilized by the optimized isokinetic Nosé-Hoover chain thermostat driven by effective solvation forces obtained from 3D-RISM-KH and extrapolated forward by generalized solvation force extrapolation (GSFE), gigantic outer time-steps up to picoseconds to accurately calculate equilibrium properties were obtained in this new quasidynamics protocol. The multiscale OIN/GSFE/3D-RISM-KH algorithm was implemented in the Amber package and well documented for fully flexible model of alanine dipeptide, miniprotein 1L2Y, and protein G in aqueous solution, with a solvent sampling rate ~150 times faster than a standard MD simulation in explicit water. Further acceleration in computation can be achieved by modifying the extent of solvation layers considered in the calculation, as well as by modifying existing closure relations. This enhanced simulation technique has proven applications in protein-ligand binding energy calculations, ligand/solvent binding site prediction, molecular solvation energy calculations, etc. Applications of the RISM-KH theory in molecular simulation are discussed in this work.

Machine learning of free energies in chemical compound space using ensemble representations: Reaching experimental uncertainty for solvation

Free energies govern the behavior of soft and liquid matter, and improving their predictions could have a large impact on the development of drugs, electrolytes, or homogeneous catalysts. Unfortunately, it is challenging to devise an accurate description of effects governing solvation such as hydrogen-bonding, van der Waals interactions, or conformational sampling. We present a Free energy Machine Learning (FML) model applicable throughout chemical compound space and based on a representation that employs Boltzmann averages to account for an approximated sampling of configurational space. Using the FreeSolv database, FML's out-of-sample prediction errors of experimental hydration free energies decay systematically with training set size, and experimental uncertainty (0.6 kcal/mol) is reached after training on 490 molecules (80% of FreeSolv). Corresponding FML model errors are on par with state-of-the art physics based approaches. To generate the input representation for a new query compound, FML requires approximate and short molecular dynamics runs. We showcase its usefulness through analysis of solvation free energies for 116k organic molecules (all force-field compatible molecules in the QM9 database), identifying the most and least solvated systems and rediscovering quasi-linear structure-property relationships in terms of simple descriptors such as hydrogen-bond donors, number of NH or OH groups, number of oxygen atoms in hydrocarbons, and number of heavy atoms. FML's accuracy is maximal when the temperature used for the molecular dynamics simulation to generate averaged input representation samples in training is the same as for the query compounds. The sampling time for the representation converges rapidly with respect to the prediction error.

Predicting PAMPA permeability using the 3D-RISM-KH theory: are we there yet?

The parallel artificial membrane permeability assay (PAMPA), a non-cellular lab-based assay, is extensively used to measure the permeability of pharmaceutical compounds. PAMPA experiments provide a working mimic of a molecule passing through cells and PAMPA values are widely used to estimate drug absorption parameters. There is an increased interest in developing computational methods to predict PAMPA permeability values. We developed an in silico model to predict the permeability of compounds based on the PAMPA assay. We used the three-dimensional reference interaction site model (3D-RISM) theory with the Kovalenko-Hirata (KH) closure to calculate the excess chemical potentials of a large set of compounds and predicted their apparent permeability with good accuracy (mean absolute error or MAE = 0.69 units) when compared to a published experimental data set. Furthermore, our in silico PAMPA protocol performed very well in the binary prediction of 288 compounds as being permeable or impermeable (precision = 94%, accuracy = 93%). This suggests that our in silico protocol can mimic the PAMPA assay and could aid in the rapid discovery or screening of potentially therapeutic drug leads that can be delivered to a desired tissue.

Application of the Approximate 3D-Reference Interaction Site Model (RISM) Molecular Solvation Theory to Acetonitrile as Solvent

The integral equation formalism based reference interaction site model (RISM) molecular solvation theory is applied to pure liquid acetonitrile and water-acetonitrile binary mixtures of different compositions. Solvate formation of d- and f-block ions by ACN is also calculated to check applicability of the RISM theory. The generalized Amber force field (GAFF) and the universal force field (UFF) parameters were found to be suitable for applications with the RISM theory for acetonitrile solvent. The presence of local microsolvated clusters for water-acetonitrile mixtures is validated in the RISM calculations.

Multisolvent Models for Solvation Free Energy Predictions Using 3D-RISM Hydration Thermodynamic Descriptors

The potential to predict solvation free energies (SFEs) in any solvent using a machine learning (ML) model based on thermodynamic output, extracted exclusively from 3D-RISM simulations in water is investigated. The models on multiple solvents take into account both the solute and solvent description and offer the possibility to predict SFEs of any solute in any solvent with root mean squared errors less than 1 kcal/mol. Validations that involve exclusion of fractions or clusters of the solutes or solvents exemplify the model's capability to predict SFEs of novel solutes and solvents with diverse chemical profiles. In addition to being predictive, our models can identify the solute and solvent features that influence SFE predictions. Furthermore, using 3D-RISM hydration thermodynamic output to predict SFEs in any organic solvent reduces the need to run 3D-RISM simulations in all these solvents. Altogether, our multisolvent models for SFE predictions that take advantage of the solvation effects are expected to have an impact in the property prediction space.

Prediction of P-glycoprotein inhibitors with machine learning classification models and 3D-RISM-KH theory based solvation energy descriptors

Development of novel in silico methods for questing novel PgP inhibitors is crucial for the reversal of multi-drug resistance in cancer therapy. Here, we report machine learning based binary classification schemes to identify the PgP inhibitors from non-inhibitors using molecular solvation theory with excellent accuracy and precision. The excess chemical potential and partial molar volume in various solvents are calculated for PgP± (PgP inhibitors and non-inhibitors) compounds with the statistical-mechanical based three-dimensional reference interaction site model with the Kovalenko-Hirata closure approximation (3D-RISM-KH molecular theory of solvation). The statistical importance analysis of descriptors identified the 3D-RISM-KH based descriptors as top molecular descriptors for classification. Among the constructed classification models, the support vector machine predicted the test set of Pgp± compounds with highest accuracy and precision of ~ 97% for test set. The validation of models confirms the robustness of state-of-the-art molecular solvation theory based descriptors in identification of the Pgp± compounds.

Application of the 3D-RISM-KH molecular solvation theory for DMSO as solvent

The molecular solvation theory in the form of the Three-Dimensional Reference Interaction Site Model (3D-RISM) with Kovalenko-Hirata (KH) closure relation is benchmarked for use with dimethyl sulfoxide (DMSO) as solvent for (bio)-chemical simulation within the framework of integral equation formalism. Several force field parameters have been tested to correctly reproduce solvation free energy in DMSO, ion solvation in DMSO, and DMSO coordination prediction. Our findings establish a united atom (UA) type parameterization as the best model of DMSO for use in 3D-RISM-KH theory based calculations.

To Pass or Not To Pass: Predicting the Blood-Brain Barrier Permeability with the 3D-RISM-KH Molecular Solvation Theory

Predicting the ability of chemical species to cross the blood-brain barrier (BBB) is an active field of research for development and mechanistic understanding in the pharmaceutical industry. Here, we report the BBB permeability of a large data set of compounds by incorporating molecular solvation energy descriptors computed by the 3D-RISM-KH molecular solvation theory. We have been able to show, for the first time, that the computed excess chemical potential in different solvents can be successfully used to predict permeability of compounds in a binary manner (yes/no) via a minimum-descriptor-based model. Our findings successfully combine the molecular solvation theory with the machine learning approach to address one of the most daunting challenges in predictive structure-activity relationship modeling. The workflow presented in this work is simple enough to be used by nonexperts with ease.