Comparison of two simulation methods to compute solvation free energies and partition coefficients

Predicting Solvation Free Energies from the Minnesota Solvation Database Using Classical Density Functional Theory Based on the PC-SAFT Equation of State

We critically assess the capabilities of classical density functional theory (DFT) based on the perturbed-chain statistical associating fluid theory (PC-SAFT) equation of state to predict the solvation free energies of small molecules in various hydrocarbon solvents. We compare DFT results with experimental data from the Minnesota solvation database and utilize statistical methods to analyze the accuracy of our approach, as well as its weaknesses. The mean absolute error of the solvation free energies is 3.7 kJ mol-1 for n-alkane solvents, ranging from pentane to hexadecane, with 473 solute-solvent systems. For solvents consisting of cyclic hydrocarbons (cyclohexane, benzene, toluene, and ethylbenzene) with 245 solute-solvent systems, we report a slightly larger mean absolute error of 4.2 kJ mol-1. We identify three possible sources of errors: (i) the neglect of solute-solvent and solvent-solvent Coulomb interactions, which limits the applicability of PC-SAFT DFT to nonpolar and weakly polar molecules; (ii) the solute's Lennard-Jones parameters supplied by the general AMBER force field, which are not parametrized toward solvation free energies; and (iii) the application of the Lorentz-Berthelot combining rules to the dispersive interactions between a segment of the PC-SAFT solvent and a Lennard-Jones interaction site of the solute. The approach is more accurate than standard implementations of phenomenological models in common chemistry software packages, which exhibit mean absolute errors larger than 9.12 kJ mol-1, even though newer phenomenological models achieve a mean absolute error of about 2 kJ mol-1. PC-SAFT DFT is more computationally efficient than state of the art explicit molecular simulations in combination with free energy perturbation methods. It is predictive with respect to solvation free energies, i.e., the input for the model is the (element-specific) molecular force field, the solute configuration from molecular dynamics simulations, and the (substance-specific) PC-SAFT parameters. The PC-SAFT parametrization uses pure-component data and does not require experimental solvation free energies. The PC-SAFT equation of state, without applying a DFT formalism, can also be used to calculate solvation free energies, provided that the PC-SAFT parameters for the solute are available. A large number of substances was recently parametrized by members of our group (Esper, T.; Bauer, G.; Rehner, P.; Gross, J. Ind. Eng. Chem. Res. 2023, 62), which enables a comparison to the DFT approach for 103 substances. Accurate results are obtained from the PC-SAFT equation of state with an MAE below 2.51 kJ mol-1. The DFT approach does not require PC-SAFT parameters for the solute and can be applied to all solutes that can be represented by the molecular force field.

Predicting octanol/water partition coefficients using molecular simulation for the SAMPL7 challenge: comparing the use of neat and water saturated 1-octanol

Blind predictions of octanol/water partition coefficients at 298.15 K for 22 drug-like compounds were made for the SAMPL7 challenge. The octanol/water partition coefficients were predicted using solvation free energies computed using molecular dynamics simulations, wherein we considered the use of both pure and water-saturated 1-octanol to model the octanol-rich phase. Water and 1-octanol were modeled using TIP4P and TrAPPE-UA, respectively, which have been shown to well reproduce the experimental mutual solubility, and the solutes were modeled using GAFF. After the close of the SAMPL7 challenge, we additionally made predictions using TIP4P/2005 water. We found that the predictions were sensitive to the choice of water force field. However, the effect of water in the octanol-rich phase was found to be even more significant and non-negligible. The effect of inclusion of water was additionally sensitive to the chemical structure of the solute.

From adaptive resolution to molecular dynamics of open systems

ABSTRACT

We provide an overview of the Adaptive Resolution Simulation method (AdResS) based on discussing its basic principles and presenting its current numerical and theoretical developments. Examples of applications to systems of interest to soft matter, chemical physics, and condensed matter illustrate the method's advantages and limitations in its practical use and thus settle the challenge for further future numerical and theoretical developments.

Performing solvation free energy calculations in LAMMPS using the decoupling approach

The decoupling approach to solvation free energy calculations requires scaling the interactions between the solute and the solution with all intramolecular interactions preserved. This paper reports a new procedure that makes it possible to these calculations in LAMMPS. The procedure is tested against built-in GROMACS capabilities. The model compounds chosen to test our methodology are ethanol and biphenyl. The LAMMPS and GROMACS results obtained are in good agreement with each other. This work should help perform solvation free energy calculations in LAMMPS and/or other molecular dynamics software having no built-in functions to implement the decoupling approach.

Predicting octanol/water partition coefficients for the SAMPL6 challenge using the SM12, SM8, and SMD solvation models

Blind predictions of octanol/water partition coefficients at 298 K for 11 kinase inhibitor fragment like compounds were made for the SAMPL6 challenge. We used the conventional, "untrained", free energy based approach wherein the octanol/water partition coefficient was computed directly as the difference in solvation free energy in water and 1-octanol. We additionally proposed and used two different forms of a "trained" approach. Physically, the goal of the trained approach is to relate the partition coefficient computed using pure 1-octanol to that using water-saturated 1-octanol. In the first case, we assumed the partition coefficient using water-saturated 1-octanol and pure 1-octanol are linearly correlated. In the second approach, we assume the solvation free energy in water-saturated 1-octanol can be written as a linear combination of the solvation free energy in pure water and 1-octanol. In all cases here, the solvation free energies were computed using electronic structure calculations in the SM12, SM8, and SMD universal solvent models. In the context of the present study, our results in general do not support the additional effort of the trained approach.

Steering a solute between coexisting solvation states: Revisiting nonequilibrium work relations and the calculation of free energy differences

By analogy with single-molecule pulling experiments, we present a computational framework to obtain free energy differences between complex solvation states. To illustrate our approach, we focus on the calculation of solvation free energies (SFEs). However, the method can be readily extended to cases involving more complex solutes and solvation conditions as well as to the calculation of binding free energies. The main idea is to drag the solute across the simulation box where atomistic and ideal gas representations of the solvent coexist at constant temperature and chemical potential. At finite pulling speeds, the resulting work allows one to extract SFEs via nonequilibrium relations, whereas at infinitely slow pulling speeds, this process becomes equivalent to the thermodynamic integration method. Results for small molecules well agree with literature data and pave the way to systematic studies of arbitrarily large and complex molecules.

Molecular Insights into Water Clusters Formed in Tributylphosphate-Di-(2-ethylhexyl)phosphoric Acid Extractant Systems from Experiments and Molecular Dynamics Simulations

Di-(2-ethylhexyl)phosphoric acid (D2EHPA) and tributylphosphate (TBP) are two of the most studied and researched organophosphorous extractants. D2EHPA is an acidic extractant, offering both hydrogen bond donor and acceptor sites while TBP, a neutral extractant, only offers a single acceptor site per molecule. In spite of this, it is observed that 1 M D2EHPA in dodecane is a poorer extractant for water than 1 M TBP in dodecane. The objective of present work is to look into the molecular interactions that cause such behavior. Experiments were carried out with varying molar ratios of TBP and D2EHPA in the organic dodecane phase. Total extractant concentration was kept constant at 1 M with dodecane as diluent. Water extraction was quantified by measuring the moisture content of the organic phase after equilibration. ¹H and ³¹P NMR spectra of the organic phase samples were recorded to study the change in the chemical environment upon extraction. Small angle X-ray scattering data of water saturated extractant phases were analyzed for the possibility of a reverse micellar aggregate formation. Molecular dynamics simulations could calculate free energies in quantitative agreement with experiments. Experimental and simulation studies showed that aggregation in the organic phase was promoted by the presence of water. This combined approach, of experiments and simulation, has shown that water is indispensable for the formation of ordered aggregates of extractants in nonpolar organic solvents. It is seen that, in the organic phase, around 80% of water's hydrogen bonds are with extractant molecules rather than with itself. The analysis clearly indicates that, rather than forming an aqueous core surrounded by extractant, water acts as a bridge between extractant molecules.

Quantum chemical and molecular dynamics modelling of hydroxylated polybrominated diphenyl ethers

A series of 19 hydroxylated polybrominated diphenyl ethers (OH-PBDEs) have been studied using density functional theory (DFT) and molecular dynamics simulations with the purpose of investigating eventual correlations between their physicochemical properties and toxic action. Dissociation constants (pK_a), solvation free energies and octanol-water partition coefficients (log P) have been computed. Additionally, metadynamics simulations of OH-PBDEs passing through a lipid bilayer have been carried out for four OH-PBDE species. No correlations between computed pK_a values and toxicity data have been found. Medium correlations were found between partition coefficients and the ability of OH-PBDEs to alter membrane potential in cell cultures, which is attributed to higher uptake of molecules with larger log P parameters. It was also demonstrated that in lipid bilayers, OH-PBDE molecules differ in their orientational distributions and can adopt different conformations which can affect the uptake of these molecules and influence the pathways of their toxic action.

Approaches for calculating solvation free energies and enthalpies demonstrated with an update of the FreeSolv database

Solvation free energies can now be calculated precisely from molecular simulations, providing a valuable test of the energy functions underlying these simulations. Here, we briefly review "alchemical" approaches for calculating the solvation free energies of small, neutral organic molecules from molecular simulations, and illustrate by applying them to calculate aqueous solvation free energies (hydration free energies). These approaches use a non-physical pathway to compute free energy differences from a simulation or set of simulations and appear to be a particularly robust and general-purpose approach for this task. We also present an update (version 0.5) to our FreeSolv database of experimental and calculated hydration free energies of neutral compounds and provide input files in formats for several simulation packages. This revision to FreeSolv provides calculated values generated with a single protocol and software version, rather than the heterogeneous protocols used in the prior version of the database. We also further update the database to provide calculated enthalpies and entropies of hydration and some experimental enthalpies and entropies, as well as electrostatic and nonpolar components of solvation free energies.

Blind prediction of cyclohexane-water distribution coefficients from the SAMPL5 challenge

In the recent SAMPL5 challenge, participants submitted predictions for cyclohexane/water distribution coefficients for a set of 53 small molecules. Distribution coefficients (log D) replace the hydration free energies that were a central part of the past five SAMPL challenges. A wide variety of computational methods were represented by the 76 submissions from 18 participating groups. Here, we analyze submissions by a variety of error metrics and provide details for a number of reference calculations we performed. As in the SAMPL4 challenge, we assessed the ability of participants to evaluate not just their statistical uncertainty, but their model uncertainty-how well they can predict the magnitude of their model or force field error for specific predictions. Unfortunately, this remains an area where prediction and analysis need improvement. In SAMPL4 the top performing submissions achieved a root-mean-squared error (RMSE) around 1.5 kcal/mol. If we anticipate accuracy in log D predictions to be similar to the hydration free energy predictions in SAMPL4, the expected error here would be around 1.54 log units. Only a few submissions had an RMSE below 2.5 log units in their predicted log D values. However, distribution coefficients introduced complexities not present in past SAMPL challenges, including tautomer enumeration, that are likely to be important in predicting biomolecular properties of interest to drug discovery, therefore some decrease in accuracy would be expected. Overall, the SAMPL5 distribution coefficient challenge provided great insight into the importance of modeling a variety of physical effects. We believe these types of measurements will be a promising source of data for future blind challenges, especially in view of the relatively straightforward nature of the experiments and the level of insight provided.

Study of procaine and tetracaine in the lipid bilayer using molecular dynamics simulation

Despite available experimental results, the molecular mechanism of action of local anesthetics upon the nervous system and contribution of the cell membrane to the process are still controversial. In this work, molecular dynamics simulations were performed to investigate the effect of two clinically used local anesthetics, procaine and tetracaine, on the structure and dynamics of a fully hydrated dimyristoylphosphatidylcholine lipid bilayer. We focused on comparing the main effects of uncharged and charged drugs on various properties of the lipid membrane: mass density distribution, diffusion coefficient, order parameter, radial distribution function, hydrogen bonding, electrostatic potential, headgroup angle, and water dipole orientation. To compare the diffusive nature of anesthetic through the lipid membrane quantitatively, we investigated the hexadecane/water partition coefficient using expanded ensemble simulation. We predicted the permeability coefficient of anesthetics in the following order: uncharged tetracaine > uncharged procaine > charged tetracaine > charged procaine. We also shown that the charged forms of drugs are more potent in hydrogen bonding, disturbing the lipid headgroups, changing the orientation of water dipoles, and increasing the headgroup electrostatic potential more than uncharged drugs, while the uncharged drugs make the lipid diffusion faster and increase the tail order parameter. The results of these simulation studies suggest that the different forms of anesthetics induce different structural modifications in the lipid bilayer, which provides new insights into their molecular mechanism.

Calculating Partition Coefficients of Small Molecules in Octanol/Water and Cyclohexane/Water

Partition coefficients describe how a solute is distributed between two immiscible solvents. They are used in drug design as a measure of a solute's hydrophobicity and a proxy for its membrane permeability. We calculate partition coefficients from transfer free energies using molecular dynamics simulations in explicit solvent. Setup is done by our new Solvation Toolkit which automates the process of creating input files for any combination of solutes and solvents for many popular molecular dynamics software packages. We calculate partition coefficients between octanol/water and cyclohexane/water with the Generalized AMBER Force Field (GAFF) and the Dielectric Corrected GAFF (GAFF-DC). With similar methods in the past we found a root-mean-squared error (RMSE) of 6.3 kJ/mol in hydration free energies which would correspond to an error of around 1.6 log units in partition coefficients if solvation free energies in both solvents were estimated with comparable accuracy. Here we find an overall RMSE of about 1.2 log units with both force fields. Results from GAFF and GAFF-DC seem to exhibit systematic biases in opposite directions for calculated cyclohexane/water partition coefficients.

Predictions of the physicochemical properties of amino acid side chain analogs using molecular simulation

A candidate drug compound is released for clinical trails (in vivo activity) only if its physicochemical properties meet desirable bioavailability and partitioning criteria.

Practical Aspects of Free-Energy Calculations: A Review

Free-energy calculations in the framework of classical molecular dynamics simulations are nowadays used in a wide range of research areas including solvation thermodynamics, molecular recognition, and protein folding. The basic components of a free-energy calculation, that is, a suitable model Hamiltonian, a sampling protocol, and an estimator for the free energy, are independent of the specific application. However, the attention that one has to pay to these components depends considerably on the specific application. Here, we review six different areas of application and discuss the relative importance of the three main components to provide the reader with an organigram and to make nonexperts aware of the many pitfalls present in free energy calculations.

Molecular dynamics simulations of water/mucus partition coefficients for feeding stimulants in fish and the implications for olfaction

The odorant partition coefficient is a physicochemical property that has been shown to dramatically influence odorant deposition patterns in the mammalian nose, leading to a chromatographic separation of odorants along the sensory epithelium. It is unknown whether a similar phenomenon occurs in fish. Here we utilize molecular dynamics simulations, based on a simplified molecular model of olfactory mucus, to calculate water/mucus partition coefficients for amino acid odorants (alanine, glycine, cysteine, and valine) that are known to elicit feeding behavior in fish. Both fresh water and salt water environments are considered. In fresh water, all four amino acids prefer the olfactory mucus phase to water, and the partition coefficient is shown to correlate with amino acid hydrophobicity. In salt water, a reversal in odorant partitioning is found, where each of the feeding stimulants (except glycine) prefer the water phase to olfactory mucus. This is due to the interactions between the salt ions and the odorant molecules (in the water phase), and between the salt and simplified mucin (in the olfactory mucus phase). Thus, slightly different odorant deposition patterns may occur in the fish olfactory organ in fresh and salt water environments. However, in both underwater environments we found that the variation of the water/mucus odorant partition coefficient is approximately one order of magnitude, in stark contrast to air/mucus odorant partition coefficients that can span up to six orders of magnitude. We therefore anticipate relatively similar deposition patterns for most amino acid odorants in the fish olfactory chamber. Thus, in contrast to terrestrial species, living in an underwater environment may preclude appreciable chromatographic odorant separation that may be used for spatial coding of odor identity across the olfactory epithelium. This is consistent with the reported lack of spatial organization of olfactory receptor neurons in the fish olfactory epithelium.