Coarse

Subscriber access provided by CARLETON UNIVERSITY

Article

Predicting partition coefficients with a simple all-atom/coarse-grained hybrid model Samuel Genheden J. Chem. Theory Comput., Just Accepted Manuscript • DOI: 10.1021/acs.jctc.5b00963 • Publication Date (Web): 20 Nov 2015 Downloaded from http://pubs.acs.org on November 29, 2015

Just Accepted “Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted online prior to technical editing, formatting for publication and author proofing. The American Chemical Society provides “Just Accepted” as a free service to the research community to expedite the dissemination of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully peer reviewed, but should not be considered the official version of record. They are accessible to all readers and citable by the Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore, the “Just Accepted” Web site may not include all articles that will be published in the journal. After a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web site and published as an ASAP article. Note that technical editing may introduce minor changes to the manuscript text and/or graphics which could affect content, and all legal disclaimers and ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or consequences arising from the use of information contained in these “Just Accepted” manuscripts.

Journal of Chemical Theory and Computation is published by the American Chemical Society. 1155 Sixteenth Street N.W., Washington, DC 20036 Published by American Chemical Society. Copyright © American Chemical Society. However, no copyright claim is made to original U.S. Government works, or works produced by employees of any Commonwealth realm Crown government in the course of their duties.

Page 1 of 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Journal of Chemical Theory and Computation

Predicting partition coefficients with a simple all-atom/coarsegrained hybrid model Samuel Genheden Department of Chemistry and Molecular Biology, Box 462, University of Gothenburg, SE-405 30 Göteborg, Sweden [email protected]

Abstract The solvation free energy is an essential quantity in force field development and in numerous applications. Here, we present the estimation of solvation free energies in polar (water, hexanol, octanol and nonanol) and in apolar (hexane, octane, nonane) media. The estimations are produced using molecular dynamics simulations employing a simple all-atom/coarse-grained hybrid model (AA/ELBA) and are therefore very efficient. More than 150 solutes were taken from the Minnesota solvation database and represent small, organic molecules. The mean absolute deviation for the different solvents ranges between 2.0 and 4.1 kJ/mol and the correlation coefficient ranges between 0.78 and 0.99, indicating that the predictions are accurate. Outliers are identified and potential avenues for improvements are discussed. Furthermore, partitions coefficients between water and the organic solvents were estimated and the percentage of the predictions that have the correct sign ranges between 74% (for octane) and 92% (for octanol and hexanol). Finally, membrane/water partition coefficients are replaced with hexane/water and octanol/water partition coefficients and the latter is found to be as accurate as the expensive membrane calculation, indicating a wider application area. Keywords: solvation free energy, partition coefficients, hybrid model, coarse-graining

Introduction The Gibbs free energy of solvation is the work needed to transfer one mol of a solute from the gas phase to solution and is a fundamental thermodynamic quantity with numerous experimental and theoretical applications. Solvation is for instance a major determinant in macromolecular complex formation and protein folding. 1 From a theoretical perspective, solvation free energies are commonly used in benchmarking molecular mechanics force fields due to the abundance of experimental data.2,3,4,5,6 In particular, the difference in free energy between two solvents with different physicochemical characteristics, i.e. the partition coefficient (log P), is a highly desirable property in for instance pharmaceutical design, separation technology or environmental sciences.7,8,9 By far the most common partition coefficient is the 1octanol/water partition coefficient (log POW). Because of its amphiphilic nature, 1-

1 Environment ACS Paragon Plus


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

octanol is a reasonable model for e.g. biological lipids and different soils. Furthermore, the hydrophobicity is usually quantified by log POW where the sign of the coefficients indicates a preference for either water or an organic phase.10 Because of its importance, several computational methods have been devised to predict its value ranging from statistical quantitative structure–property relation (QPSR) models11,12 to physics-based simulation methods. 13,14 Other solvents apart from 1– octanol are also used and one of the more commonly used is n-hexane, forming a nhexane/water partition coefficient (log PHW). N-hexane is good model of an apolar, hydrophobic phase and experimental data of log PHW has been used previously as input to QPSR11,12 and simulations methods.15 Furthermore, a biologically important partition coefficient is that between water and a membrane (log PMW). 16 The biological membrane is a heterogeneous phase with several distinct regions and is the natural place of many important processes in both biochemistry and biotechnology. Unfortunately there are relatively few experimental partition coefficients available, and therefore log PMW is usually approximated with either log POW or log PHW, 17 although this is questionable.16 The computational methods18,19 available to estimate log PMW also tends to be much more expensive than the methods used to estimate POW or log PHW. Recently, an efficient all-atom/coarse-grained (AA/CG) hybrid model was presented and applied to membrane processes. 20 For instance it was used to successfully predict partition coefficients for eleven small organic molecules in a model membrane. Because of this success to predict log PMW, it is relevant to ask if the hybrid model also can be used to predict simpler partition coefficients like POW and log PHW. In this paper, this question is investigated. It is also of general interest to expand the applicability of the AA/CG model and to see how well it performs in comparison to for instance all-atom models. Furthermore, there have not been many large-scale benchmark studies on solvation free energies in organic solvents,6 although partition coefficients are extremely important as discussed above. Unfortunately, there are no available CG models of solvents like n-hexane or 1octanol in this particular CG force field. 21 Therefore, simple, ad-hoc models were developed of hexane, hexanol, nonane and nonanol, where the latter two also can be used to model octane and octanol, respectively. This simple and computational efficient method was shown to be accurate in predicting partition coefficients. Method Solvent molecules. Coarse-grained (CG) models representing 1-hexanol, 1-nonanol, nhexane and n-nonane were created with the Elba force field.21,22 The Elba lipid force field uses a 3 to 1 mapping of the tails, i.e. three carbons are mapped to a single CG bead. Therefore, two such beads represent a hexane molecule and three such beads represent a nonane molecule. All bonded and non-bonded parameters were taken from the lipid force field, implying that the beads are uncharged Lennard-Jones sites held together with harmonic bond and angle potentials. To create a model for hexanol and nonanol, an additional bead was appended to hexane and nonane, respectively, that has the same parameters as the Elba water bead. As such, it has a point dipole on a single Lennard-Jones site. Bonded parameters were taken from the Elba lipid models. The solvent molecules are depicted in Figure 1 and all parameters are listed in Supporting Info, Table S1. Due to the uniform size of the CG beads, it is unfortunately not possible to create unique models for octane and octanol. However, octanol is an important solvent and there are plenty of experimental data available.


Page 2 of 16

Page 3 of 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Therefore, estimates for solutes in octane and octanol were produced with the same solvent model, i.e. with identical representation and force field parameters, as for nonane and nonanol, respectively.

Figure 1 – Models of solvent molecules. Each sphere represents a Lennard-Jones particle and the arrows indicate point dipoles. Solvent simulations. For each solvent molecule a box with 512 solvent molecules (or 1024 in the case of water) was created using Packmol (version 15.217).23 The size of the box was set according to the experimental density.24 This system was equilibrated for 100 ps in the NVT ensemble and for 100 ps in the NPT ensemble, followed by a 10 ns production simulation in the NPT ensemble. All simulations were performed with the Lammps software25 with a 10 fs timestep. The temperature was fixed at 298 K using Langevin dynamics26 with a 10 ps time constant, and the pressure was fixed at 1 atm using a weak-coupling scheme27 with a 1 ps time constant. Snapshots were saved for analysis every 2 ps. The density was computed from the instantaneous volume and the enthalpy of vaporization was calculated using established protocols.28,29 A modified version of Lammps as well as equilibrated boxes can be downloaded from http://github.com/sgenheden. Solute molecules. Solute molecules were taken from the Minnesota solvation database (version 2012).30,31 For each of the solvents hexanol, octanol, nonanol, hexane, octane and nonane, all solutes that also had a determined solvation free energy in water was extracted. Furthermore, all solutes with more than 10 heavy atoms were excluded to limit the calculations to cases were sampling issues are minimised. Furthermore, hydrogen, water and tetramethylsilane were removed leaving the total number of unique solute molecules at 168. The Cartesian coordinates provided in the database were used to parameterise the solutes with the general Amber force field (GAFF)32 and AM1-BCC charges 33 using AmberTools (version 15). 34 The in-house python scripts used to setup the solutes are available from http://github.com/sgenheden. Six additional solutes (dibutylether, N,N-dimethylacetamide, N,Ndimethylformamide, 1-chlorobutane, dimethylsulfoxide, diisopropylether) were parameterised and simulated in octanol in order to make comparisons with the recent data of Zhang et al.6 Furthermore, six molecules with experimental membrane/water partition coefficients20 (glycerol, benzylalcohol, p-chlorocresol, 2,4,5-trichloroaniline, hexachlorobenzene) were parameterised and simulated. The Cartesian coordinates for these twelve additional molecules were generated from SMILES strings and parameterised with the same protocol as described above. Solvation of solutes. The solutes were inserted in the pre-equilibrated solvent boxes. To make this automatic, a spherical void with a ≈1.0 nm diameter was created at the centre of the solvent box by growing a Lennard-Jones particle from a non-interacting dummy to a fully interacting site during a 1 ns NPT simulation with Lammps.25 Each of the solutes was then inserted at the centre of this void (after removal of the 3 Environment ACS Paragon Plus


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Page 4 of 16

Lennard-Jones particle), avoiding any overlap between solvent and solute atoms. This has the same effect as overlaying the solutes with pre-equilibrated boxes and removing solvent atoms that overlap with the solute. The solute–solvent system was then subjected to 100 steps of steepest descent minimisation and equilibration for 1.5 ns in the NPT ensemble. A recently presented multiple timestep algorithm20 was used in these simulations to propagate the CG–CG forces with a 6 fs timestep and all other forces with a 2 fs timestep. Covalent bonds to hydrogen atoms within the solutes were constrained with the SHAKE algorithm.35 The solvent and solutes were coupled to two different Langevin thermostats26 keeping the temperature fixed at 298 K. Solvation free energy estimates. The Gibbs free energy of solvation was estimated using thermodynamic integration (TI)36

−∆Gsolv =

∫

1 0

∂U ∂λ

dλ

1)

λ

were U is the potential of the system, and λ is a coupling variable. At λ = 0, the solute interacts fully with the solvent molecules, and at λ = 1, the solute is in a decoupled state and does not interact at all with the solvent molecules. Because soft-core potentials37,38 are currently not implemented in Lammps, U was scaled with a fourthpower function f (λ ) = (1− λ )4 .39,40 Twenty-five equally spaced values of λ from 0 to 0.96 was used and λ = 1 was estimated by linear extrapolation, and the integration in Equation 1 was carried out using the trapezoid rule. One 60 ns long simulation was carried out and the value of λ was changed step-wise, every 2.4 ns. The initial 600 ps at each value of λ were discarded as equilibration and the last 1.8 ns were used to estimate the ensemble average in Equation 1. The sampling frequency was 0.6 ps. The setup of the simulations was otherwise identical to the equilibration simulations discussed above. Two independent repeats of the TI simulations were performed by both randomly rotating the solutes in the solvent box prior to equilibration and by assigning different random seeds to the velocity generation. Error analysis. The agreement between the TI estimates and the experimental data was quantified using mean absolute deviation (MAD), Pearson’s correlation coefficient (R) and Kendall’s τ.41 Systematic errors caused by specific functional groups were analysed using the approach suggested by Mobley et al. 42 Therefore, the BEDROC (Boltzmannenhanced discrimination of receiver-operating characteristic) metric43 was computed for different functional group as outlined previously.42,6 The Checkmol program (version 0.5)44 was used to identify the functional groups, and the BEDROC analysis was carried out with the CROC python package (version 1.1). 45 The uncertainty estimate of the BEDROC metric was estimated by 500 bootstrap iterations. A Student’s t-test was also performed on the absolute errors for the different chemical groups compared to the entire population of absolute errors. Table 1 – Density (ρ) and heat of vaporization (∆Hvap) for the modeled solvents. Solvent Hexane Nonane

ρ [g/l] Model Experiment24 747.5 660.6 767.4 719.2

∆Hvap [kJ/mol] Model Experiment46 33.2 31.6 52.6 46.4


Page 5 of 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Hexanol Nonanol Water

779.3 798.4 1003.1

813.6 828.0 997.0

63.8 79.1 39.2

61.6 76.9 44.0

Results and discussion Pure liquids. The density (ρ) and heat of vaporization (∆Hvap) for the five model solvents are listed in Table 1. The liquid properties of water are reproduced well as discussed previously.47 For the apolar liquids, both ρ and ∆Hvap are overestimated. The relative error in ρ is 14 and 7% for hexane and nonane, respectively, and the relative error in ∆Hvap is 5 and 13%, respectively. These errors are acceptable considering the simplicity of the model, viz. the solvents are modeled with either two or three Lennard-Jones beads. The relative errors in density are larger than has been observed previously for atomistic solvents,29 but the relative errors in ∆Hvap are on a par with atomistic solvents. Surprisingly, the liquid properties of hexanol and nonanol are better reproduced than for the corresponding apolar solvents. The densities are underestimated slightly, but the relative errors are only 4 and 3% for hexanol and nonanol, respectively. Furthermore, the relative errors in ∆Hvap are only around 4%, or 2 kJ/mol. Therefore, it must be concluded that these simple and relatively ad-hoc models of polar solvents are accurate enough to model pure liquids. Because of the relative success of reproducing liquid properties with the simple models, not further optimization was attempted. Table 2 – Statistics for the estimates of solvation free energy (∆Gsolv) and partition coefficient (log P) ∆Gsolv (scaled) ∆Gsolv log P N MAD MAD Solvent solutes [kJ/mol] R τ [kJ/mol] R τ MAD R Accuracy1 Hexane 54 6.1 0.78 0.61 2.98 0.78 0.60 0.86 0.84 80% Octane 35 6.9 0.86 0.68 2.15 0.86 0.66 0.78 0.80 74% Nonane 25 7.3 0.89 0.73 2.09 0.88 0.73 0.62 0.81 80% Hexanol 13 5.1 0.98 0.84 2.98 0.97 0.92 0.79 0.61 92% Octanol 166 5.4 0.91 0.71 3.64 0.91 0.72 0.66 0.75 92% Nonanol 9 3.4 0.99 0.78 3.22 0.99 0.83 0.84 0.69 78% Water 168 4.1 0.93 0.75 1 Defined as the percentage of solutes for which the estimate of log P has the correct sign

Solvation free energies. The calculated solvation free energies are plotted against the experimental data in Figure S1 and statistics of the performance are shown in Table 2. Full results can be found in Supporting Information. Although, there are 168 unique solutes, there are not experimental solvation free energies for all of them in all solvents. The number of solutes for each solvent is listed in Table 2 as well, ranging from 9 for nonanol to 168 for water. Furthermore, as mentioned in the Methods section, it is not possible to construct a unique model for octane and octanol, so those estimates were produce using the same models as for nonane and nonanol, respectively. Therefore, the different results for octane and nonane as well as for octanol and nonanol is due to the different set of solutes and the experimental solvation free energy. For the apolar solvents, the mean absolute deviations (MAD) compared to experiment ranges between 6 and 7 kJ/mol, which indicate significant errors. 5 Environment ACS Paragon Plus


1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

However, the correlation coefficients (R) range between 0.78 and 0.89, and are statistical significant. Kendall’s τ also indicates a positive and significant ranking of the solutes. Therefore it seems that there is a systematic error in the estimates leading to too negative estimates for the solutes, which is also clear from the correlation plots in Figure S1. For the estimates of solutes in hexanol, octanol and nonanol, the MAD ranges only between 3 and 5 kJ/mol, which is significantly lower than for the corresponding apolar solvents. The correlation is also stronger, with correlation coefficients above 0.90. Also, τ is more positive for the polar solvents. This could be an effect of that the pure liquid properties are better reproduced for the hexanol and nonanol than for hexane and nonane. Still, it is clear from the correlation plots in Figure S1 that many of the estimates are too negative. The experimental solvation free energy of solutes in polar liquids is usually for the hydrated liquid, i.e. with a significant mole fraction of water. However, we have computed the solvation free energies in anhydrous liquids, a common approximation in the literature.6,13 Typically, the difference between the two solvation free energies are within the uncertainty of the calculations.13,48 Finally, the estimates in water give an MAD of 4.1 kJ/mol and an R of 0.93, which is on a par with previous benchmarking studies in the literature. For instance the MAD and R of the more than 600 solutes in the FreeSolv database49 is 4.8 kJ/mol and 0.94, respectively.

Figure 2 – Correlation between experimental and calculated solvation free energy. The solvent is indicated in the upper-left corner of each subplot. In order to improve the absolute estimates of the predictions in the organic solvents, the van der Waals interaction between the solvent and the solute was reduced. This is


Page 6 of 16

Page 7 of 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


similar to what was done previously to improve the barrier height of the potential of mean force for solutes in membranes.20 This re-parameterization can be seen as a compensation factor for the simplified physical description of CG–AA interactions. The parameterization was performed on the solutes in octane to limit the risk of overfitting. Reduction with 50% and 25% produced too positive estimates of the free energy of solvation and large MADs (see Table S2), but a reduction with 10%, i.e. a scaling factor of 0.9 for the Lennard-Jones ε-parameter gave improved estimates compared to experiments. Therefore, this scaling factor was applied in estimations of solvation free energy in all organic solvents. It should be emphasized that only the interaction between the solute atoms and the apolar CG beads was modified and that the interaction with the dipolar beads was left unchanged. The re-calculated solvation free energies are plotted against the experimental data in Figure 2 and statistics of the performance are shown in Table 2. It is clear that R and τ are largely unaffected by the modified parameter. The largest change is τ for estimates in hexanol that changes from 0.84 to 0.92, and τ for estimates in nonanol that changes from 0.78 to 0.83. However, the MAD is reduced significantly for all solvents. As expected, the estimates in apolar solvents are most affected with improvements in MAD between 3 and 5 kJ/mol, whereas the improvements in MAD for the polar solvents are smaller. Still, this simple modification leads to estimates that at least on average are within chemical accuracy. Therefore, this modification will be used throughout the rest of this report. For the solutes in octanol, a direct comparison with the recent atomistic study of Zhang et al. can be made.6 They computed the solvation free energy of 29 solutes and 23 of them are also included in the Minnesota solvation database.30 Therefore, ∆Gsolv was calculated for the six molecules not in the database using the same protocol as for the other solutes. The full results are listed in Table S3. For these 29 solutes, the hybrid model gives a MAD of 3.0 kJ/mol compared to 2.6 with a fully atomistic study. This difference is statistical significant, although it is not clear if the difference has any practical effect. However, R for the hybrid model is 0.85, comparable to the 0.90 for the atomistic study, indicating that hybrid model is equally capable of ranking the solutes, although it has a larger absolute error. Table 3 – BEDROC analysis for identified chemical groups.1 Group alcohol aldehyde alkane alkene alkyl bromide alkyl chloride amine aromatic compound carbonitrile carboxylic acid carboxylic acid ester ether halogen derivative

Hexane Uniform Observed 0.44 0.20 ±0.09

Octanol Uniform Observed 0.43 0.34 ±0.07 0.42 0.40 ±0.09 0.44 0.31 ±0.07 0.43 0.61 ±0.09 0.43 0.59 ±0.07 0.42 0.65 ±0.09 0.43 0.40 ±0.08

Uniform 0.43 0.42 0.44 0.43 0.43 0.42 0.43

Water Observed 0.38 ±0.07 0.67 ±0.13 0.30 ±0.05 0.57 ±0.07 0.27 ±0.06 0.41 ±0.08 0.79 ±0.07

0.48

0.83

±0.06

0.46 0.42 0.42

0.58 0.43 0.51

±0.06 ±0.15 ±0.05

0.46 0.42 0.42

0.47 0.41 0.59

±0.06 ±0.12 ±0.06

0.44

0.46

±0.08

0.43 0.43

0.38 0.29

±0.06 ±0.08

0.43 0.43

0.22 0.61

±0.07 ±0.07

0.43

0.66

±0.08

0.43

0.29

±0.08



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

heterocyclic compound 0.43 0.51 ±0.10 0.43 0.53 ±0.09 ketone 0.44 0.13 ±0.06 0.43 0.08 ±0.03 0.43 0.53 ±0.08 nitro compound 0.42 0.12 ±0.04 0.42 0.38 ±0.11 phenol 0.43 0.98 ±0.04 0.42 0.20 ±0.06 0.42 0.91 ±0.05 1 Both the expected BEDROC value from a uniform distribution and the observed one are shown. Observed BEDROC values that are greater than and significantly different from the uniform ones are indicated in bold.

Error analysis. The error distribution of the calculated solvation free energies are shown in Figure S2 and it can be seen that there a couple of outliers for each solvent except for hexanol and nonanol. In Table S4 the predictions with the largest discrepancy compared to experiment are listed. For the apolar solvents, many of the worst predictions are for aromatic compounds, and for the polar solvents, nitrogenand halogen-containing compounds, are clearly identified. In addition, for water one can identify several ethers and aldehydes. Table 4 – p-value of a t-test of the unsigned deviation and mean signed error (MSE) in kJ/mol for identified chemical groups1 Hexane p-value MSE 1.3 0.03

Octanol p-value MSE 0.40 -2.1 0.71 3.4 0.08 -0.4 0.08 5.3 0.07 4.8 0.08 5.4 0.64 2.0

Water p-value MSE 0.90 -0.2 0.18 -2.9 0.2 0.02 5.3 0.04 2.1 0.03 0.90 3.6 4.9 0.00

Group N N N alcohol 7 16 16 aldehyde 5 5 alkane 19 19 alkene 10 10 alkyl bromide 10 10 alkyl chloride 6 6 amine 12 12 aromatic compound 20 4.6 42 0.08 4.2 42 0.89 -1.7 0.03 carbonitrile 5 0.82 4.1 5 0.88 3.5 carboxylic acid 5 0.09 -4.1 5 -5.6 0.02 carboxylic acid ester 8 0.92 -2.9 9 0.23 -2.8 11 -2.1 0.01 ether 15 0.15 1.5 15 5.2 0.04 halogen derivative 16 2.3 16 0.35 -0.7 0.03 heterocyclic compound 14 0.38 3.1 14 0.36 4.7 ketone 6 0.7 10 -0.1 10 0.30 4.6 0.00 0.00 nitro compound 7 -1.1 7 0.58 -0.2 0.00 phenol 5 0.05 8.0 7 0.5 7 -9.1 0.01 0.00 1 p-values smaller than 0.05 are indicated in bold. N is the number of solutes for each chemical group.

In order to proceed with a more formal analysis, the BEDROC metric was computed for the seventeen chemical groups identified in the solutes. Only groups that were represented at least five times were included and therefore, it was only possible to perform the analysis on all seventeen groups for octanol and water. In Table 3, the BEDROC metric is listed for these solvents as well as hexane. The p-values for a ttest of the absolute deviation for the individual chemical groups compared to the total population of deviations were also calculated. These p-values are listed in Table 4


Page 8 of 16

Page 9 of 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


together with the mean signed error (MSE). For hexane, the aromatic group shows an observed BEDROC value of 0.83, which is greater than the value expected from a uniform distribution and is in line with the qualitative observations in Table S4. Furthermore, the t-test of the absolute deviations shows that aromatic compounds are significantly different from the full population. Furthermore, phenols are also identified with a BEDROC value of 0.98, also in line with the qualitative observations. The group MSE of 8 kJ/mol is also the largest observed for hexane. For solutes in octanol, alkyl chlorides, alkyl bromides, aromatic compounds and ethers all gives an observed BEDROC value that is greater than that expected from a uniform distribution. This is in line with the compounds identified in Table S4. However, it is interesting to note that the BEDROC values are smaller than the outliers identified for hexane. Some of the identified groups also has a MSE different from the overall MSE as indicated by p-values lower than 0.05 and large MSE. For water, alkenes, amines, carboxylic acids, ethers and phenol compounds gives large observed BEDROC values, significant p-values and sizable MSE. Thus it is clear that different types of solvents are more sensitive to different kinds of solutes and it is difficult to identify clear starting points for improvements. It should also be mentioned that the analysis is dependent on the modification introduced to the solute–solvent interactions, as different groups were deemed outliers when the original parameterisation was used (not shown). Furthermore, it is clear that the groups identified here as outliers was not identified in the recent all-atom study of Zhang et al.,6 which also indicates that the main error stems from solute–solvent interactions rather than from problems with the solute force field. Partition coefficients. From the solvation free energy in water and an organic solvent S, the partition coefficient can be expressed as50

log P =

∆Gsolv (water) − ∆Gsolv (S) 2.3RT

(2)

were R is the gas-constant and T the absolute temperature. The calculated partition coefficients are plotted against the experimental coefficients in Figure 3 and statistics of the performance are listed in Table 2. For the apolar solvents the MAD ranges between 0.6 and 0.9 log units, and R ranges between 0.80 and 0.84. This is on a par with the correlations observed for the solvation free energies. If as a measure of success, we consider the sign of the partition coefficient is correctly predicted, the accuracy is between 74 and 80%, which is fairly accurate. It should be noted that this calculation has not excluded cases where the estimate is within the experimental uncertainty, which is on the order of 0.2 log units given an uncertainty in the free energies of 0.8 kJ/mol. For the polar solvents the MAD is similar to the MAD for apolar solvents and ranges between 0.61 and 0.75. The R is only 0.66 to 0.84, which is slightly worse than the correlation observed for the solvation free energy. Still, the accuracy of the estimates is between 78 and 92%, indicating that the hybrid model can correctly predict partitioning.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Figure 3 – Correlation between experimental and calculated solvent/water partition coefficients. The solvent is indicated in the upper-left corner of each subplot. Because of the relative success in predicting simple partition coefficients it is of interest if they can be used as proxies for other partition coefficients that are more expensive to compute. One example is the calculation of membrane/water partition coefficients (log PMW) that are interesting in for instance drug design. The AA/Elba model was used recently to estimate log PMW for eleven small molecules in a DMPC (1,2-dimyristoyl-sn-glycero-3-phosphocholine) membrane. The results were encouraging and on a par with all-atom simulations.20 Here, octanol/water (log POW) and hexane/water (log PHW) partition coefficients were computed for these eleven molecules because these two coefficients have previously been used as proxies for log PMW. About half of the molecules were already included in the Minnesota solvation database. The correlation between log PMW and either log PHW or log POW is illustrated in Figure 4. It is clear that hexane is a worse proxy for the membrane than octanol, indicating that the head group influences the partitioning significantly. It is also clear that the log POW estimate are systematically shifted which is an effect of the re-parameterized CG–AA interactions. The membrane simulations were performed without any modification to the interaction between the tail beads and the solute.20 In fact, log POW estimates without the modified interaction correlate better with log PMW (not shown). Because log POW correlates better with log PMW than log PHW, it is not surprising that the correlation to the experimental partition coefficients is much better for log POW than for log PHW as can be seen in Figure 4. In fact, the MAD is only 1.4 log unit, which is only slightly worse than the MAD for log PMW (see Table 5). Also R is 0.75, which is very similar to R for the calculated log PMW, whereas R for log PHW is only 0.62. Finally, the accuracy is 55 and 100% is for log PHW and log POW, respectively, indicating that the octanol/water partition coefficient is an excellent proxy for these molecules. However, a larger validation set is necessary to consolidate this.


Page 10 of 16

Page 11 of 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Figure 4 – Correlations between calculated and experimental partition coefficient (log P). a) log PHW versus calculated log PMW b) log POW versus calculated log PMW c) log PHW versus experimental log PMW and d) log POW versus experimental log PMW



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

Table 5 – Statistics for the estimates of membrane/water, hexane/water and octanol/water partition coefficient (log P) for 11 solutes with experimental membrane/water data MAD R Accuracy1 20 log PMW 1.1 0.78 91% log PHW 2.9 0.62 55% log POW 1.4 0.75 100% 1

Defined as the percentage of solutes for which the estimate of log P has the correct sign

Conclusions In this paper we report on the use of a simple AA/CG hybrid model20 to estimate solvation free energies and partition coefficients. The solvent environment is modelled using a simple CG model made up of Lennard-Jones particles and in the case of polar solvents, an additional dipolar bead. This leads to a considerable reduction in the number of particles that needs to be simulated, and thereby very efficient simulations. The solutes, not always being amenable to CG are described using a common AA force field and combined with the CG model as recently demonstrated for hydration free energies and membrane transfer free energies. The hybrid model utilises a multiple timestep algorithm20 further improving the efficiency. A single estimate of the solvation free energy takes between 10 and 15 hours on 12 cores. Considering the simplicity of the solvent environment and the computational efficiency, the method seems promising. The pure liquid properties are reproduced well but further improvements for hexane and nonane can probably be made. Coincidental, it is also the estimate of ∆Gsolv and log P in the apolar solvents that are the least satisfactory. The MAD is up to 7 kJ/mol for ∆Gsolv, although the correlation is significant and fairly good, which indicates a systematic error in the absolute estimates. Therefore, the organic solvents were made less sticky by reducing the van der Waals interaction between the solvent and the solutes. This lead to a MAD for all solvent that was at most 4 kJ/mol, i.e. reaching chemical accuracy. Partition coefficients were also accurately estimated with a MAD up to 0.9 log units and correlations between 0.61 and 0.84. This is really encouraging, considering the ad-hoc nature of the solvent model. It is basically apolar or polar beads held together with harmonic forces and no additional effort has been put into the parameterisation of the liquid force fields. Finally, an attempt was made to replace the expensive calculations of membrane/water partition coefficients with the simpler log PHW and log POW. This turned out to be an excellent approximation for log POW and the results were of the same predictive power as the calculated log PMW. This is encouraging and warrants further studies with larger test sets. However, it remains to be seen if this AA/CG hybrid model is practically useful in prediction partition coefficients. At the moment it is slightly less accurate in predicting absolute values than an all-atom simulations,6 although the ranking is on a par with the more realistic all-atom solvent model. Because it is more computational efficient than an AA simulation, it is an alternative in that perspective, although it is much slower than semi-empirical 51 or a quantitative structure–property relation models52. Still, this study show that a relative simple physical model of a solvent environment can be used predict the correct partitioning for up to 92% of the solutes and that it has a potential use in replacing expensive membrane simulations.


Page 12 of 16

Page 13 of 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


Acknowledgements SG acknowledges the Wenner–Gren foundations for financial support and C3SE at Chalmers Technical University and iSolutions at the University of Southampton for computational resources. Leif A. Eriksson and Ulf Ryde are acknowledged for fruitful discussions and encouragement. Supporting Information ELBA force field parameters for the CG liquids, additional correlation plots and statistics are provided in a Microsoft Word document. Full lists with results are provided as a Microsoft Excel spreadsheet. This information is available free of charge via the Internet at http://pubs.acs.org References 1

Levy, Y.; Onuchic, J. N. Water mediation in protein folding and molecular recognition. Annu. Rev. Biophys Biomol. Struct., 2006, 35, 389–415. 2 Shirts, M. R.; Pitera, J. W.; Swope, W. C.; Pande, V. S. Extremely precise free energy calculations of amino acid side chain analogs: Comparison of common molecular mechanics force fields for proteins. J. Chem. Phys., 2003, 119, 5740. 3 Guthrie, J. P. A blind challenge for computational solvation free energies: introduction and overview. J. Phys. Chem. B, 2009, 113, 4501–4507. 4 Shivakumar, D.; Harder, E.; Damm, W.; Friesner, R. A.; Sherman, W. Improving the Prediction of Absolute Solvation Free Energies Using the Next Generation OPLS Force Field. J. Chem. Theory Comput., 2012, 8, 2553–2558. 5 Knight, J. L.; Yesselman, J. D.; Brooks, C. L. Assessing the quality of absolute hydration free energies among CHARMM-compatible ligand parameterization schemes. J. Comput. Chem., 2013, 34, 893–903. 6 Zhang, J.; Tuguldur, B.; van der Spoel, D. Force Field Benchmark of Organic Liquids II: Gibbs Energy of Solvation. J. Chem. Inf. Model., 2015, 55, 1192–1201. 7 Leo, A.; Hansch, C.; Elkins, D. Partition coefficients and their uses. Chem. Rev., 1971, 71, 525–616. 8 Essex, J. W.; Reynolds, C. A., & Richards, W. G. Theoretical determination of partition coefficients. J. Am. Chem. Soc., 1992, 114, 3634–3639. 9 Schulte, J.; Dürr, J.; Ritter, S.; Hauthal, W. H.; Quitzsch, K.; Maurer, G. Partition Coefficients for Environmentally Important, Multifunctional Organic Compounds in Hexane + Water. J. Chem. Eng. Data, 1998, 43, 69–73. 10 Ruelle, P. The n-octanol and n-hexane/water partition coefficient of environmentally relevant chemicals predicted from the mobile order and disorder (MOD) thermodynamics. Chemosphere, 2002, 40, 457–512. 11 Goss, K.-U.; Schwarzenbach, R. P. Linear Free Energy Relationships Used To Evaluate Equilibrium Partitioning of Organic Compounds. Environ. Sci. Technol, 2001, 35, 1–9. 12 Nguyen, T. H.; Goss, K.-U.; Ball, W. P. Polyparameter Linear Free Energy Relationships for Estimating the Equilibrium Partition of Organic Compounds between Water and the Natural Organic Matter in Soils and Sediments. Environ. Sci. Technol., 2005, 39, 913–924. 13 Bhatnagar, N.; Kamath, G.; Chelst, I.; Potoff, J. J. Direct calculation of 1-octanol-water partition coefficients from adaptive biasing force molecular dynamics simulations. J. Chem. Phys., 2012, 137, 014502 14 Duffy, E. M.; Jorgensen, W. L. Prediction of Properties from Simulations: Free Energies of Solvation in Hexadecane, Octanol, and Water. J. Am. Chem. Soc., 2000, 122, 2878–2888. 15 Garrido, N. M.; Economou, I. G.; Queimada, A. J.; Jorge, M.; Macedo, E. A. Prediction of the n -hexane/water and 1-octanol/water partition coefficients for environmentally relevant compounds using molecular simulation. AIChE J., 2012, 58, 1929–1938.



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

16

Klamt, A.; Huniar, U.; Spycher, S.; Keldenich, J. COSMOmic: a mechanistic approach to the calculation of membrane-water partition coefficients and internal distributions within membranes and micelles. J. Phys. Chem. B, 2008, 112, 12148–12157. 17 Gobas, F. A. P. C.; Mackay, D.; Shiu, W. Y.; Lahittete, J. M.; Garofalo, G. A novel method for measuring membrane-water partition coefficients of hydrophobic organic chemicals: Comparison with 1-octanol-water partitioning. J. Pharma. Sci., 1988, 77, 265–272. 18 Bemporad, D.; Luttmann, C.; Essex, J. W. Computer simulation of small molecule permeation across a lipid bilayer: dependence on bilayer properties and solute volume, size, and cross-sectional area. Biophys. J., 2004, 87, 1–13. 19 Neale, C.; Bennett, W. F. D.; Tieleman, D. P.; Pomès, R. Statistical Convergence of Equilibrium Properties in Simulations of Molecular Solutes Embedded in Lipid Bilayers. J. Chem. Theory Comput., 2011, 7, 4175–4188. 20 Genheden, S.; Essex, J. A simple and transferable all-atom/coarse-grained hybrid model to study membrane processes. J. Chem. Theory Comput., 2015, 11, 4749–4759 21 Orsi, M.; Essex, J. W. The ELBA force field for coarse-grain modeling of lipid membranes. PloS One, 2011, 6, e28637. 22 Orsi, M.; Essex, J. W. Physical properties of mixed bilayers containing lamellar and nonlamellar lipids: insights from coarse-grain molecular dynamics simulations. Faraday Discus., 2013, 161, 249-272 23 Martínez, L.; Andrade, R.; Birgin, E. G.; Martínez, J. M. PACKMOL: a package for building initial configurations for molecular dynamics simulations. J. Comput. Chem., 2009, 30, 2157–2164. 24 Lide, D. R. CRC Handbook of Chemistry and Physics 85th edition; CRC Press: Cleveland, Ohio, 2004. 25 Plimpton, S. Fast Parallel Algorithms for Short-Range Molecular Dynamics. J. Comput. Phys., 1995, 117, 1–19. 26 Hünenberger, P. H. Thermostat algorithms for molecular dynamics simulations. Ad. Polymer Sci., 2005, 173, 105–147. 27 Berendsen, H. J. C.; Postma, J. P. M.; van Gunsteren, W. F.; DiNola, A.; Haak, J. R. Molecular dynamics with coupling to an external bath. J. Chem. Phys., 1984, 81, 3684. 28 Oostenbrink, C.; Villa, A.; Mark, A. E.; van Gunsteren, W. F. A biomolecular force field based on the free enthalpy of hydration and solvation: the GROMOS force-field parameter sets 53A5 and 53A6. J. Comput. Chem., 2004, 25, 1656–1676. 29 Caleman, C.; Van Maaren, P. J.; Hong, M.; Hub, J. S.; Costa, L. T.; Van Der Spoel, D. Force field benchmark of organic liquids: Density, enthalpy of vaporization, heat capacities, surface tension, isothermal compressibility, volumetric expansion coefficient, and dielectric constant. J. Chem. Theory Comput., 2012, 8, 61–74. 30 Marenich, A. V.; Kelly, C. P.; Thompson, J. D.; Hawkins, G. D.; Chambers, C. C.; Giesen, D. J.; Winget, P.; Cramer, C. J.; Truhlar, D. G. Minnesota Solvation Database – version 2012, University of Minnesota, Minneapolis, 2012. 31 Marenich, A. V.; Cramer, C. J.; Truhlar, D. G. Universal Solvation Model Based on Solute Electron Density and on a Continuum Model of the Solvent Defined by the Bulk Dielectric Constant and Atomic Surface Tensions. J. Phys. Chem. B, 2009 113, 6378–6396. 32 Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem., 2004, 25, 1157–1174. 33 Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, efficient generation of high-quality atomic charges. AM1-BCC model: II. Parameterization and validation. J. Comput. Chem., 2002, 23, 1623–1541 34 Salomon-Ferrer, R.; Case, D. A.; Walker, R. C. An overview of the Amber biomolecular simulation package. WIREs: Comput. Mol. Sci., 2013, 3, 198–210. 35 Ryckaert, J.-P.; Ciccotti, G.; Berendsen, H. J. C. Numerical integration of the cartesian equations of motion of a system with constraints: molecular dynamics of n-alkanes.J. Comput. Phys., 1977, 23, 327–341.


Page 14 of 16

Page 15 of 16

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60


36

Kirkwood, J. G. Statistical mechanics of fluid mixtures, J. Chem. Phys., 1935, 3, 300-313 Beutler, T. C.; Mark, A. E.; van Schaik, R. C.; Gerber, P. R.; van Gunsteren, W. F. Avoiding singularities and numerical instabilities in free energy calculations based on molecular simulations. Chem. Phys. Lett., 1994, 222, 529–539. 38 Zacharias, M.; Straatsma, T. P.; McCammon, J. A. Separation-shifted scaling, a new scaling method for Lennard-Jones interactions in thermodynamic integration. J. Chem. Phys., 1994, 100, 9025 39 Steinbrecher, T.; Mobley, D. L.; Case, D. A. Nonlinear scaling schemes for Lennard-Jones interactions in free energy calculations. J. Chem. Phys. 2007, 127, 214108. 40 Orsi, M.; Ding, W.; Palaiokostas, M. Direct Mixing of Atomistic Solutes and CoarseGrained Water. J. Chem. Theory Comput., 2014, 10, 4684–4693. 41 Kendall, M. A New Measure of Rank Correlation. Biomet. 1938, 30, 81–89. 42 Mobley, D. L.; Bayly, C. I.; Cooper, M. D.; Shirts, M. R.; Dill, K. A. Small molecule hydration free energies in explicit solvent: An extensive test of fixed-charge atomistic simulations. J. Chem. Theory Comput., 2009, 5, 350–358. 43 Truchon, J.-F.; Bayly, C. I. Evaluating virtual screening methods: good and bad metrics for the “early recognition” problem. J. Chem. Inf. Model., 2007, 47, 488–508. 44 Haider N. Checkmol. http://merian.pch.univie.ac.at/~nhaider/cheminf/cmmm.html Accessed 08/14/2015 45 Swamidass, S. J.; Azencott, C.-A.; Daily, K.; Baldi, P. A CROC stronger than ROC: measuring, visualizing and optimizing early retrieval. Bioinformatics, 2010, 26, 1348–1356. 46 Majer, V.; Svoboda, V., Enthalpies of Vaporization of Organic Compounds: A Critical Review and Data Compilation, Blackwell Scientific Publications, Oxford, 1985, 300. 47 Orsi, M. Comparative assessment of the ELBA coarse-grained model for water. Mol. Phys., 2013, 112, 1–11. 48 Bernazzani, L.; Cabani, S.; Conti, G.; Mollica, V. Thermodynamic Study of the Partitioning of Organic Compounds between Water and Octan-1-Ol. Effects of Water as Cosolvent in the Organic Phase. J. Chem. Soc. Faraday Trans. 1995, 91, 649-655. 49 Mobley, D. L.; Guthrie, J. P. FreeSolv: a database of experimental and calculated hydration free energies, with input files. J. Comput-Aided Mol. Des., 2014, 28, 711–720. 50 Michel, J., Orsi, M., Essex, J. W. Prediction of partition coefficients by multiscale hybrid atomic-level/coarse-grain simulations. J. Phys Chem. B., 2008, 112, 657–660. 51 Klamt, A. Conductor-like Screening Model for Real Solvents: A New Approach to the Quantitative Calculation of Solvation Phenomena. J. Phys. Chem., 1995, 99, 2224–2235. 52 Katritzky, A. R.; Oliferenko, A. A.; Oliferenko, P. V.; Petrukhin, R.; Tatham, D. B.; Maran, U.; Lomaka, A.; Acree, W. E. A general treatment of solubility. 1. The QSPR correlation of solvation free energies of single solutes in series of solvents. J. Chem. Inf. Comput. Sci., 2003, 43, 1794–1805. 37



1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60

For Table of Contents Only

“Predicting partition coefficients with a simple all-atom/coarse-grained hybrid model” – Samuel Genheden


Page 16 of 16

Coarse

Recommend Documents