Effect of Atomic Charges on Octanol-Water Partition Coefficient Using Alchemical Free Energy Calculation.

The octanol-water partition coefficient (logPow) is an important index for measuring solubility, membrane permeability, and bioavailability in the drug discovery field. In this paper, the logPow values of 58 compounds were predicted by alchemical free energy calculation using molecular dynamics simulation. In free energy calculations, the atomic charges of the compounds are always fixed. However, they must be recalculated for each solvent. Therefore, three different sets of atomic charges were tested using quantum chemical calculations, taking into account vacuum, octanol, and water environments. The calculated atomic charges in the different environments do not necessarily influence the correlation between calculated and experimentally measured ∆Gwater values. The largest correlation coefficient values of the solvation free energy in water and octanol were 0.93 and 0.90, respectively. On the other hand, the correlation coefficient of logPow values calculated from free energies, the largest of which was 0.92, was sensitive to the combination of the solvation free energies calculated from the calculated atomic charges. These results reveal that the solvent assumed in the atomic charge calculation is an important factor determining the accuracy of predicted logPow values.

1. Introduction

The octanol–water partition coefficient (logP) is an index used in drug discovery processes to estimate the solubility, membrane permeability, and bioavailability of compounds [1,2]. Several experimental methods for measuring it have been developed [2,3,4,5,6], but they are time consuming and costly, because experimental methods use the actual compounds, and sometimes, the synthesis of a compound takes a long time. To accelerate the drug discovery process and to reduce its cost, we need a logP measurement method that is easier, faster, and more accurate than the current ones.

Computational methods for predicting the logP value have been developed using various techniques. The methods summing the contributions of a compound’s atoms and fragments (e.g., XLOGP3 [7], KOWWIN [8], etc.) rapidly calculate logP values similar to those obtained experimentally, and have therefore been used in the screening of drug-like compounds [9], making it possible to calculate the logP values for the tens of millions of compounds in databases such as PubChem [10] and ZINC [11].

The logP values of compounds can also be predicted accurately by alchemical free energy calculation techniques using molecular dynamics (MD) simulation and Monte Carlo simulation, such as thermodynamic integration (TI) [12] and free energy perturbation (FEP) [13]. These methods have been used to calculate hydration free energy of amino acids [14], nucleobases [15], and ligand/protein binding free energy [16,17,18,19,20,21,22,23,24,25,26,27,28]. They can treat water molecules either implicitly or explicitly. Thus, in the case of ligand/protein binding free energy calculation, they can be used to predict the thermodynamic behavior of a ligand in solvents. In compounds containing water molecules bridging between the protein and ligand, the calculation with explicit water molecules is efficient for calculating binding free energy. In the application of alchemical free energy calculation to logP calculation, DeBolt and Kollman successfully estimated the logP using FEP method with modified force field [29]. They were able to reproduce the structural and thermodynamic behaviors of liquid octanol. Using an implicit solvation model, Huang et al. reported a logP prediction method in which the 3D-RISM-KH theory was used to calculate the solvation free energy in water and in octanol [30]. To reproduce the experimental conditions, Bhatnagar et al. reported an adaptive biasing force method using explicit water and octanol solvents [31]. In another application, saturated octanol solvation was used for the calculation, taking the hydrate state of octanol in the experiment into account [32,33]. These methods were able to estimate logP values in agreement with the ones measured experimentally. The alchemical methods using all-atom MD simulation require time and computer resources, compared with atom- and fragment-based methods. However, these methods are able to predict the partition coefficient of several solvents, and observe the thermodynamic behavior of compounds in the solvents. Therefore, the free energy calculation using the alchemical method is expected to provide useful information on physical chemistry properties of the compounds.

In this paper, the logP values for relatively small compounds were predicted from the solvent free energies calculated by TI [34,35]. The calculation of the difference free energy, ∆G, was estimated by the Bennett’s acceptance ratio (BAR) method [36]. To analyze the contribution of the atomic charges of the compounds, three sets of atomic charges were examined: one in vacuum, one in octanol, and the other in water. Solvation free energies in water and in octanol are in good agreement with the ones measured experimentally. The correlation coefficient (R) of the solvation free energies indicated strong correlation between calculated and experimentally measured values. They did not differ significantly among the three sets of atomic charges. The logP values, however, showed slightly different behaviors. In fact, the R values of logP depended on the combination of the atomic charges calculated in the different solvents. These results reveal that the environment assumed in the calculation of the atomic charges is a very critical factor for the logP calculations. Therefore, our results provide useful information for improving the quality of the calculations. Finally, the logP values of 17 drugs calculated using the best combination of the atomic charges showed good agreement with the values measured experimentally.

2. Results and Discussion

2.1. Solvation Free Energy in Water

To calculate free energy in water, the three charge sets were used independently. In this paper, the notation of the parameter sets calculated with the atomic charges in the different solvents are simplified as P{x}, where x is either v(vacuum), o(octanol), or w(water). It should be noted that, in this work, we did not change the atomic charges in water molecules, because the force field of water molecules was established by the many studies conducted on water molecules [37].

Strong correlations between calculated and experimentally measured ∆G values were observed (Figure 1). In particular, the free energy calculated using P{w}, for which the coefficient of determination (R2) is 0.87, indicated good agreement with experimentally measured values. The free energy calculated using P{v} also showed a strong correlation with the experimental data, even though R2 (=0.84) was the smallest in the three cases. These results confirm that the calculated atomic charges in the different environments do not necessarily influence the correlation between calculated and experimentally measured ∆G values.

Figure 1

Scatter diagram of calculated and experimentally measured ∆G values in water solvent. R2 and regression line are shown.

On the other hand, the values of atomic charge were sensitive to the solvent in the quantum chemistry calculation (Figure 2). In the three solvents, somewhat large differences in atomic charge values can be found in the charged group containing heteroatoms, e.g., oxygen in acetone, nitrogen in pyridine, and oxygen in octanol. The calculated atomic charges in octanol and water were similar, and were slightly different from the ones in the vacuum state. This tendency can be understood from the graphs drawn for the atomic charges vs. dielectric constant values (Figure S1). The atomic charges dramatically increase or decrease at dielectric constants around 10, and become constant thereafter. These differences can affect the regression line in the scatter diagrams (Figure 1). The slope of the regression line for ∆G values with P(w) is close to that with P(o), whereas that with P(v) is slightly different from the others. However, the R values do not differ significantly among the three parameters. Furthermore, there were strong correlations among three.

Figure 2

Calculated atomic charges of (a) Acetone, (b) Pyridine, and (c) Octanol in vacuum, octanol, and water.

∆G values, for which the R values of P(v)–P(o), P(v)–P(w), and P(o)–P(w) pairs, were 0.9912, 0.9885, and 0.9996, respectively. In these cases, the order in which the compounds are arranged according to the ∆G, is more or less the same among the three pairs. These analyses show that the atomic charges were affected by the solvent types assumed in the calculation, but these differences did not affect the correlation of the ∆G values.

The RMSE and MAE values in the three parameter sets show opposite relations to the R values; i.e., the calculated data having strong correlation showed large errors (Table 1). The ∆G of P{v} had the smallest R and R2 values, but the errors, for which the RMSE and MAE values are 4.49 kJ/mol (=1.07 kcal/mol) and 3.65 (=0.87 kcal/mol), respectively, were the smallest in the three parameter sets. On the other hand, the calculation with P{w} shows the largest R value, but the RMSE and MAE values are the largest. In these calculations, the difference of RMSE and MAE values, however, are 1.2 kJ/mol and 0.6 kJ/mol between P{v} and P{w}. These may not be significant differences among the three parameters. Therefore, for all the selections of atomic charges, the BAR method can estimate the ∆G values accurately enough.

Table 1

Summary of free energy in water calculations.

Parameter	R	R2	RMSE a (kJ/mol)	MAE b (kJ/mol)
Pwater{v}	0.92	0.84	4.47	3.63
Pwater{o}	0.93	0.87	5.05	3.74
Pwater{w}	0.93	0.87	5.69	4.28

a Root mean square error; b Mean average error.

2.2. Solvation Free Energy in Octanol

In calculating the free energy in octanol, nine sets of the atomic charges were used: all combinations of the three sets of atomic charges of the compound (3 sets) and octanol (3 sets). The atomic charges of octanol molecule in octanol solvation, the octanol environment, should be considered. However, in this work, to understand the effect of the atomic charges in octanol on the logP calculation, we took into account the atomic charges in the octanol molecule for vacuum, octanol, and water state.

The compound and octanol parameter sets calculated with the atomic charges in octanol is defined as P{x, y}, in which x and y (for compound and octanol, respectively) are either v(vacuum), o(octanol), or w(water).

Tendencies similar to those seen in the results obtained with the ∆G calculation could be seen in the results obtained with the ∆G calculation. It can be seen in Figure 3 that all the calculated ∆G values showed strong correlations with the ones measured experimentally. The largest R value was 0.9 (P{v, o} and P{v, w}), though these values are slightly smaller than that in ∆G (0.93). The calculated ∆G values distributed toward the small range of the horizontal axis, thus, all the calculations may underestimate the ∆G value. The tendency of the relationship between the R2 values and error values (RMSE and MAE) is similar to that seen in the ∆G results; i.e., the larger the R2 values are, the larger the RMSE and MAE values are. The largest MAE is 13.26 kJ/mol (=3.16 kcal/mol) in the ∆G with P{w, w}. This value is a little bit larger than that in ∆G. This result reveals that the calculation of ∆G values accumulates more errors than that of ∆G. However, these ∆G values still showed strong correlation with the experimental data.

Figure 3

Scatter diagram of calculated and experimentally measured ∆G values in octanol solvent. R2 and regression line are shown.

In Figure 3 and Table 2, the correlations between the calculated and experimentally measured ∆G values are shown to be independent of the atomic charges. The values of R and R2 are close in all the ∆G calculations. The R, R2, RMSE, and MAE values of the ∆G with P{o, o} are comparable with the other ones. This means that the atomic charges did not affect the correlation of ∆G calculations. The ∆G calculated with either P{v, o} or P{v, w} had the largest R and R2 values. Those values indicated a strong correlation with the experimental data. On the other hand, the ∆G calculated with P{w, v} had the smallest R and R2 values. However, these values still indicated strong correlation. The differences of R and R2 values between P{v, o} (or P{v, w}) and P{w, v}, are 0.02 and 0.04, respectively. These are not necessarily significant differences. Therefore, it is possible to conclude that the calculated atomic charges in the different environments did not affect the R and R2 values in ∆G calculations.

Table 2

Summary of free energy in octanol calculations.

Parameter	R	R2	RMSE a (kJ/mol)	MAE b (kJ/mol)
Poctanol{v, v}	0.89	0.80	11.18	10.36
Poctanol{v, o}	0.90	0.81	11.23	10.49
Poctanol{v, w}	0.90	0.81	11.26	10.56
Poctanol{o, v}	0.88	0.78	13.59	12.53
Poctanol{o, o}	0.89	0.78	13.70	12.77
Poctanol{o, w}	0.89	0.79	13.65	12.72
Poctanol{w, v}	0.88	0.77	14.18	13.05
Poctanol{w, o}	0.89	0.79	14.15	13.18
Poctanol{w, w}	0.89	0.79	14.23	13.26

a Root mean square error; b Mean average error.

2.3. Effect of Water and Octanol Solvents on Free Energy

To investigate the effects of solvents on free energy calculation, the dielectric constant of water and octanol solvents were evaluated in the simulation. The dielectric constant was calculated as [38,39] where V, k and T are volume, Boltzmann constant, and temperature, respectively. M is the total dipole moment. To calculate the dielectric constant of water and octanol solvents, additional 5 ns MD simulations were performed to generate the products, after 10,000 steps optimization, 50 ps heat-up, and 2 ns equilibration. The dielectric constants were calculated for the trajectories of 5 ns MD simulations.

Overestimated dielectric constants of water solvents were obtained by simulation (Table 3). This value is very similar to that reported by Vega et al. [40]. In Figure 1, the largest slope of the regulations line is found in the calculation with P(v), of which absolute atomic charge values was smallest in the three parameters. Therefore, the pair of the water solvent with which the dielectric constant was underestimated, and the compound with P(v) provided the closest absolute values of ∆G in our calculation.

Table 3

Properties of solvents in simulation.

Solvents	Charge	ε		<μ> (Debye)	Flu c (Debye2)	<V> (Å3)	vdW (kJ/mol)	Elec. (kJ/mol)
		Exp.	Sim.
Water		78.49 (25 °C) a	92.70	2.35	20,914.2	22,290.7	14,371.9	−33,641.5
Octanol	Vacuum	9.34 (32.1 °C) b	4.21	2.08	836.1	26,547.7	−3860.9	−2778.5
	Octanol		3.63	2.24	678.0	26,320.5	−3768.3	−3687.8
	Water		3.42	2.27	625.2	26,347.2	−3714.9	−3879.3

a Uematsu and Frank [41]; b Smyth and Stoops [42,43]; c Flu = − 2.

On the other hand, the dielectric constants of octanol solvent were also underestimated. All the values in the different charges were smaller than the experimentally measured ones. These values were slightly smaller than that in the papers reported by DeBolt and Kollman [29]. They obtained dielectric constants near 5.1 for the liquid octanol at 40 °C. The dielectric constant in octanol solvent gradually decreased by the atomic charges in the order of vacuum, octanol, and water, i.e., these values decreased in proportion to the absolute values of atomic charges. In Table 3, <μ> of octanol with atomic charge in water was larger than that of vacuum because of increasing absolute atomic charges. This means that, owing to increasing <μ>, the values of electrostatic interaction increased. On the other hand, the van der Waals interaction drops with increasing absolute atomic charges. Since the van der Waals parameters were the same among the different simulations, the electrostatic interaction is thought to dominate the simulation. However, in Figure 3, the slopes of the regulation line were dependent on the atomic charges of compounds, rather than those of octanol solvents. Moreover, the dielectric constant with atomic charges in water solvent further differed from the experimentally measured values than that within vacuum. In the simulation of liquid octanol, of which the dielectric constant of the octanol is relatively smaller than that of water solvent, a strong contribution of the electrostatic interaction leads to smaller dielectric constant values that appear as simulation errors. Therefore, in order to improve the calculated free energy values, the balance between electrostatic and van der Waals interactions may need to be investigated in further studies.

2.4. LogPow for Tested Compounds

The logP values were calculated using the Formula (7) (See Section 3). Here, the notation for the logP value calculated from ∆G using P{x}, and ∆G using P{y, z} was simplified to P{x, {y, z}}, where x, y, and z are either v(vacuum), o(octanol), or w(water).

The logP values were sensitive to the atomic charges calculated in the solvent. The correlation of the logP with the experimental data is strongly dependent on the pair of solvation free energies ∆G and ∆G (Figure 4 and Figure S2). For instance, in Figure 4, completely different correlations of logP values can be seen in the scatter diagrams with the largest and smallest R2, which are P{w, {v, o}} and P{v, {w, o}}, respectively. In the case of P{w, {v, o}}, the logP was calculated from ∆G with P{w}, and ∆G with P{v, o}. This combination is expected to result in logP values closely correlated with the experimental ones.

Figure 4

Scatter diagram of calculated and experimentally measured logP values. The diagrams with the largest and smallest R2 are shown, and the diagram of P{w, {w, w}} are shown for reference. R2 and regression line are shown.

On the other hand, in the case of P{v, {w, o}}, the logP values showed very weak correlation with the experimental data. The calculated logP values were based on ∆G with P{v} and ∆G with P{w,o}. The free energies in water and octanol both showed strong correlations with the experimental data, and the RMSE and MAE values were not as significantly different from the others, as mentioned above. The slope of the regression line for ∆G with P{v}, however, was over 1, whereas that of the regression line for ∆G with P{w,o} was under 1. These differences raised the large RMSE and MAE values slightly (Table 4). As a result, the logP values showed very weak correlation. These results reveal that the combination of the ∆G and ∆G is very important for the correct logP. That is, the selection of the solvent in the calculation of the atomic charges is very important for the logP calculation. Thus, the inappropriate selection of the solvents for the calculation of atomic charges results in the weak correlation between the logP and the experimental data. This tendency can be observed in the logP with P{w,{w, w}}. The ∆G with P{w}, and ∆G with P{w, w}, showed the strongest correlation in our test data. The values, however, showed moderate correlation (R = 0.77) not stronger than that of the ∆G with P{w, {v, o}}.

Table 4

Summary of logP calculations.

Parameters	R	R2	RMSE 1 (kJ/mol)	MAE 2 (kJ/mol)
PlogP{v, {v, v}}	0.77	0.59	2.16	2.04
PlogP{v, {v, o}}	0.80	0.64	2.17	2.07
PlogP{v, {v, w}}	0.82	0.68	2.17	2.07
PlogP{v, {o, v}}	0.55	0.31	2.61	2.42
PlogP{v, {o, o}}	0.60	0.36	2.63	2.46
PlogP{v, {o, w}}	0.62	0.38	2.62	2.45
PlogP{v, {w, v}}	0.49	0.24	2.72	2.51
PlogP{v, {w, o}}	0.54	0.29	2.72	2.53
PlogP{v, {w, w}}	0.54	0.29	2.74	2.55
PlogP{o, {v, v}}	0.92	0.84	1.40	1.31
PlogP{o, {v, o}}	0.92	0.85	1.42	1.34
PlogP{o, {v, w}}	0.92	0.85	1.43	1.35
PlogP{o, {o, v}}	0.85	0.73	1.79	1.69
PlogP{o, {o, o}}	0.87	0.76	1.82	1.73
PlogP{o, {o, w}}	0.88	0.77	1.81	1.72
PlogP{o, {w, v}}	0.82	0.67	1.90	1.78
PlogP{o, {w, o}}	0.84	0.71	1.90	1.81
PlogP{o, {w, w}}	0.84	0.71	1.92	1.82
PlogP{w, {v, v}}	0.92	0.85	1.26	1.17
PlogP{w, {v, o}}	0.92	0.86	1.29	1.20
PlogP{w, {v, w}}	0.92	0.85	1.31	1.21
PlogP{w, {o, v}}	0.88	0.78	1.63	1.55
PlogP{w, {o, o}}	0.90	0.80	1.66	1.59
PlogP{w, {o, w}}	0.90	0.81	1.66	1.58
PlogP{w, {w, v}}	0.86	0.75	1.73	1.64
PlogP{w, {w, o}}	0.88	0.77	1.75	1.66
PlogP{w, {w, w}}	0.88	0.77	1.76	1.68

1 Root mean square error; 2 Mean average error.

The chemical group’s effect on the logP values of a compound is shown in Figure 5. In the logP calculation with P{w, {v, o}}, the errors of each compound were ~2.0 in all groups, except the ethers (for ethyl phenyl ether, the error was 3.16). The calculated ∆G value of ethyl phenyl ether was an overestimate. The calculated and the experimentally measured ∆G values were −10.22 kJ/mol and −17.92 kJ/mol, respectively. On the other hand, the calculated ∆G was an underestimate. The calculated and experimentally measured ∆G values were −34.91 kJ/mol and −23.65 kJ/mol, respectively. These differences caused large errors in logP values. Similar errors can be found in the logP values calculated with P{w, {w, w}}. The profile of histograms for P{w, {w, w}} is similar to that of histograms for P{w, {v, o}}. The error of ethyl phenyl ether stands out from the other compounds. The error of each compound was around 2.0, but these errors are slightly larger than that of the values calculated with P{w, {v, o}}. Therefore, the R values for P{w, {w, w}} were slightly smaller than those for P{w, {v, o}}.

Figure 5

Error of logP with (a) P, (b) P, and (c) P for the experimental values. Here, the left-to-right order of the compounds in each group is the same as the top-down order in Table S1.

On the other hand, most of the logP values calculated with P{v, {w, o}} were larger than those calculated with the other parameters. Those of the compounds in the phosphate group, which has a strong polarity at the phosphoric acid moiety, showed particularly large errors. The selected environments in the charge calculation was inconsistent, i.e., the vacuum state was selected in ∆G calculation and water environment was selected in ∆G calculations. These inadequate combinations induce large errors in the logP calculation. A similar tendency can be found in the other groups with polar moieties, such as ether and ketone groups. The selection of the environment of the calculated charges, therefore, is important in logP calculation.

To analyze the contributions of the atomic charges in different environments, the R and R2 values of logP in Table 4 were compared. In the calculated ∆G, the R and R2 values become small in the order of P{v}, P{o}, and P{w} (Figure S3a). In terms of the atomic charges of the compounds in ∆G calculations, the R and R2 values of logP become small in the order of water, octanol, and vacuum environments (Figure S3b). However, the atomic charges of the octanol in octanol solvent seem to affect the logP values very slightly (Figure S3c). To briefly summarize these results, the atomic charges of the compound in the ∆G and ∆G calculations are better in water and vacuum environments, respectively. In free energy, the atomic charges of the octanol molecule in octanol solvent should be reasonable, nonetheless, they do not influence the logP value much. It should be noted that this tendency will be independent of the functional and basis sets in the calculation of atomic charges, because the RMSE shows the same behavior (Figure S4), and the correlation efficient values of our procedure are very large with the other procedures using the different functionals and basis sets in five compounds (Table S4).

2.5. LogPow for 17 Compounds

The best combination of the parameters of the additional 17 compounds (P{w, {v, w}}) was different from those of the small compounds (P{w, {v, o}}) (Table S5). However, the R and R2 values of P{w, {v, v}}, P{w, {v, o}}, and P{w, {v, w}} were very close to each other. Therefore, we conclude that choosing one of P{w, {v, *}} (* = v, o, and w) parameter set will provide R and R2 values close to the experimental ones.

The calculated logP values of the additional 17 compounds with the best combination of the atomic charges in water and octanol (P{w, {v, w}}) showed a good agreement with the experimental data (Figure 6). The R and R2 values were 0.93 and 0.86, respectively. These values were similar to those obtained by the other methods based on the compounds’ atoms and fragments [44], though our method needs more computational time than the other methods. The RMSE and MAE were 2.58 and 2.26, respectively, which were slightly larger than those of the small compounds. These errors affect the slope of the regression line, which is somewhat smaller than that of the line obtained for the small compounds when using P{w, {v, o}} (Figure 4). These results indicate that the order of the compounds sorted by the logP values, is very similar to the order of the compounds sorted by the logP values measured experimentally. This means that the alchemical free energy calculation method with P{w, {v, w}} will be effective to predict the lipophilicity of compounds, and will be a useful tool for drug design.

Figure 6

Scatter diagram of calculated and experimentally measured logP values for 17 compounds.

The logP values of the other additional 32 compounds were calculated for comparing our results with the previous results reported by Bannan et al. [45]. The logP values show strong correlation with the experimentally measured values, of which the R and R2 values were 0.85 and 0.72, respectively (Figure 7). These values are very close to those calculated using GAFF-DC that were reported by Bannan et al., of which R and R2 values are 0.85 and 0.74, respectively. These results reveal that the consideration of the atomic charge is important for improving the logP values.

Figure 7

Scatter diagrams of calculated using (a) GAFF, (b) GAFF-DC, and (c) P{w, {v, o}}, and experimentally measured logP values for 32 compounds. The values calculated by GAFF and GAFF-DC, and the experimentally measured values were taken from the reference by Bannan et al. [45].

3. Materials and Methods

3.1. Setting Parameters of Compounds

The general amber force field (GAFF) [46] was used for bonding parameters (bond length, bond angle, torsion angle, and improper torsion angle) and for nonbonding van der Waals parameters (radius and well-depths) of the compounds and 1-octanol. Atomic charges were calculated by using Gaussian09 software [47] at the B3LYP/6-31G(d) level. The polarizable continuum model (PCM) [48] was used to consider the solvent effect on the calculation of atomic charges. Antechamber software was used for restrained electrostatic potential (RESP) charge-fitting [49,50,51] from the results of quantum chemical calculation. In the preparations for the simulations, biomolecule’s atomic charges are calculated once in one solvent, i.e., in water solvent, in protein, or in vacuum state. When calculating a compound’s free energy in octanol, however, better results can be obtained using the atomic charges calculated for the compound in octanol. Therefore, in this study, the atomic charges of the compounds and 1-octanol were calculated for three states: in vacuum, in 1-octanol, and in water. In the calculation of the free energy in water, the three charge sets were used independently. To calculate the free energy in octanol, nine sets of the atomic charges were used: all combinations of the three sets of atomic charges of the compound (3 sets) and octanol (3 sets). The compound and octanol parameter sets calculated with the atomic charges in octanol is defined as P{x, y}, in which x and y (for compound and octanol, respectively) are either v(vacuum), o(octanol), or w(water). The logP values based on these free energies were calculated. The dependency of the calculated atomic charges in the different environments on the free energy and logP values were examined.

3.2. Initial Models of Compounds in Solvent Box

We constructed the model systems as follows. All the compounds were hydrated using a TIP3P [37] water box consisting of 723 water molecules, and all the compounds were placed into an octanol box consisting of 100 1-octanol molecules.

All the systems were minimized by 10,000 steps of the steepest-descent method, first with constraint on the positions of the heavy atoms, and then without constraint. After the minimization, 7 ns MD simulation for the compound was performed with Gromacs 4.5.3 software suite in a periodic boundary condition, by using particle-mesh Ewald (PME) calculations [52,53] for coulombic interactions. The temperature was set at 298.15 K by using a stochastic dynamics algorithm, and pressure was maintained at 1 bar by scaling the box size in accordance with the Berendsen algorithm [54]. The frictional coefficient γ and time step were 5 ps−1 and 2.0 fs. Snapshot structures were extracted from equilibration runs as initial structures for successive annihilation MD simulations. The snapshot after 2 ns of MD simulation was used as the initial structure for the free energy calculation.

3.3. Free Energy Calculation with BAR Method

Annihilation with the BAR method is widely used in the calculation of ligand/protein binding free energy [25,26,27,28], and is briefly explained here. We consider the changes in free energy for transformation from vacuum to solvent using the BAR method, based on the following formulas: where ΔG represents the free energy change from the compound in gas (S) to the compound in solvent (S) in Formula (3), and ΔG represents the free energy change from the compound in vacuum (S0) to the compound in gas (S) in Formula (4). In Formula (3), to calculate the interaction by considering the reverse processes, the first annihilation is performed for the compound from solvent to gas. And in Formula (4), the second annihilation is performed to remove the internal energy in the compounds. In these processes, the free energy is estimated by summing the free energy differences between neighboring partially annihilated states (here states are denoted as windows). The free-energy difference between the neighboring i-th and (i + 1)-th windows is calculated by using the BAR method [55]. In the BAR method, for U representing the potential energy of the i-th window, nonequilibrium work W (=U+1–U) was sampled in the MD trajectory. According to the maximum likelihood argument of BAR, the free-energy estimator for two equal-length simulations can be written as where ΔG is the annihilation free-energy change from the i-th window to the (i + 1)-th window, and where a pair of angle brackets indicates an ensemble average over a window. The free energy difference is then defined as follows:

In the free energy calculation using the BAR method, the annihilating simulations with the soft core potential [56,57] were independently performed for the sets of λ points composing λC and λLJ for Coulomb and van der Waals interactions, respectively. With λLJ = 0.0, the λC was increased in steps of 0.05 from 0 to 1, and with λC = 1, the λLJ was increased in steps of 0.05 from 0 to 1. Totally, 41 λ-points were considered in this paper.

For each window, 1 ns runs for data collection were performed for water and octanol solvent after 1 ns MD simulation for equilibration. Nonequilibrium work between neighboring windows was measured every 2 ps. Other conditions for each simulation were the same as those for the simulation generating initial structure, as mentioned above. To reduce artifacts due to the initial structure selection and to improve the statistical reliability, three independent sets of annihilation simulations were performed.

3.4. Calculation of LogPow

The logP values were calculated using the following formula: where R and T are the gas constant and temperature, respectively, and ΔG and ΔG are the free energy differences (calculated using the BAR method as described) between the compound in solvent and the compound in vacuum. Here, the notation for the logP value calculated from ΔG using P{x}, and ΔG using P{y, z}, was simplified to P{x, {y, z}}, where x, y, and z are either v(vacuum), o(octanol) or w(water). Equation (6) indicates that we should consider the effects of free energy differences caused by the perturbation of the atomic charge, by the change of the environments, and by the different atomic charges. In this work, these free energy differences were neglected, because they cannot be accurately reproduced using currently available methods.

3.5. Test Compounds

The 58 relatively small sized compounds examined to measure the logP are listed in Table S1, along with their experimentally measured free energies extracted from a paper by Wang et al. [58]. In each compound, logP values were calculated by the combination of the free energies of the water and octanol solvents, i.e., three P{x} and nine P{x, y} parameters. Finally, 27 logP values were calculated in total.

The logP values of an additional 17 compounds were calculated using the same combination of the atomic charges in water and octanol for free energy calculations, as mentioned above. These compounds (Table S2) were used as drugs in the work reported by Daina et al. [44]. The calculation conditions used were the same as those used for the other 58 compounds.

As additional tests, the logP values of 32 compounds tested in the previous work reported by Bannan et al. [45], were calculated using the best combination of the atomic charge calculations. In this case, the P{w} and P{v, o} showed the best correlation, and were used in additional calculations. All the calculations were performed with the same conditions as the previous work by Bannan et al. [45], except for the set of λ points, for which the values of this work were used, as mentioned above.

4. Conclusions

A logP calculation method using alchemical free energy calculation with the BAR method was proposed, and the results showed good agreement with the values measured experimentally. The free energies in water and octanol both showed strong correlations with the experimental data. However, the correlation of logP values was sensitive to the combination of the ΔG and ΔG values. The weak correlation was due to the different slopes of regression line between ΔG and ΔG. Finally, we found that the calculation with P{w, {v, o}} provides the best correlations between the calculated logP values and experimentally measured values.

We calculated the logP values for the drugs using the best parameter set in our trials, and obtained values closely correlated with the ones obtained experimentally. The R values of logP values were better than those obtained in previous studies. These results show that the solvent assumed in the charge calculation greatly influences the logP values. The calculation of the atomic charges was time consuming, but dramatically increased the accuracy of the logP values. We expect these accurate results to be useful in the field of drug discovery, and to accelerate drug design processes.