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

Abstract

The octanol–water partition coefficient (logP_ow) is an important index for measuring solubility, membrane permeability, and bioavailability in the drug discovery field. In this paper, the logP_ow 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 ∆G_water 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 logP_ow 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 logP_ow values.

1. Introduction

The octanol–water partition coefficient (logP_ow) 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_ow measurement method that is easier, faster, and more accurate than the current ones.

Computational methods for predicting the logP_ow 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_ow 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_ow values for the tens of millions of compounds in databases such as PubChem [10] and ZINC [11].

The logP_ow 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_ow calculation, DeBolt and Kollman successfully estimated the logP_ow 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_ow 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_ow 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_ow 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_ow values, however, showed slightly different behaviors. In fact, the R values of logP_ow 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_ow calculations. Therefore, our results provide useful information for improving the quality of the calculations. Finally, the logP_ow 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_water{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_water values were observed (Figure 1). In particular, the free energy calculated using P_water{w}, for which the coefficient of determination (R²) is 0.87, indicated good agreement with experimentally measured values. The free energy calculated using P_water{v} also showed a strong correlation with the experimental data, even though R² (=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_water values.

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_water values with P_water(w) is close to that with P_water(o), whereas that with P_water(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.

∆G_water values, for which the R values of P_water(v)–P_water(o), P_water(v)–P_water(w), and P_water(o)–P_water(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_water, 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_water 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. 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.

). The ∆G_water of P_water{v} had the smallest R and R² 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_water{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_water{v} and P_water{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_water values accurately enough.

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_ow 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_octanol{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_water calculation could be seen in the results obtained with the ∆G_octanol calculation. It can be seen in Figure 3 that all the calculated ∆G_octanol values showed strong correlations with the ones measured experimentally. The largest R value was 0.9 (P_octanol{v, o} and P_octanol{v, w}), though these values are slightly smaller than that in ∆G_water (0.93). The calculated ∆G_octanol values distributed toward the small range of the horizontal axis, thus, all the calculations may underestimate the ∆G_octanol value. The tendency of the relationship between the R² values and error values (RMSE and MAE) is similar to that seen in the ∆G_water results; i.e., the larger the R² 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_octanol with P_octanol{w, w}. This value is a little bit larger than that in ∆G_water. This result reveals that the calculation of ∆G_octanol values accumulates more errors than that of ∆G_water. However, these ∆G_octanol values still showed strong correlation with the experimental data.

In Figure 3 and

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.

, the correlations between the calculated and experimentally measured ∆G_octanol values are shown to be independent of the atomic charges. The values of R and R² are close in all the ∆G_octanol calculations. The R, R², RMSE, and MAE values of the ∆G_octanol with P_octanol{o, o} are comparable with the other ones. This means that the atomic charges did not affect the correlation of ∆G_octanol calculations. The ∆G_octanol calculated with either P_octanol{v, o} or P_octanol{v, w} had the largest R and R² values. Those values indicated a strong correlation with the experimental data. On the other hand, the ∆G_octanol calculated with P_octanol{w, v} had the smallest R and R² values. However, these values still indicated strong correlation. The differences of R and R² values between P_octanol{v, o} (or P_octanol{v, w}) and P_octanol{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 R² values in ∆G_octanol calculations.

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]

$$\epsilon =1+\frac{4\pi }{3k_{b}T\langle V\rangle }\left\{\langle M^2\rangle -{\langle M\rangle }^2\right\}$$

(1)

where V, k_b 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. 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 = <M²> − <M>².

). 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_water(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_water(v) provided the closest absolute values of ∆G in our calculation.

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. LogP_ow for Tested Compounds

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

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

On the other hand, in the case of P_logP{v, {w, o}}, the logP_ow values showed very weak correlation with the experimental data. The calculated logP_ow values were based on ∆G_water with P_water{v} and ∆G_octanol with P_octanol{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_water with P_water{v}, however, was over 1, whereas that of the regression line for ∆G_octanol with P_octanol{w,o} was under 1. These differences raised the large RMSE and MAE values slightly (

Table 4. Summary of logP_ow 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

¹ Root mean square error; ² Mean average error.

). As a result, the logP_ow values showed very weak correlation. These results reveal that the combination of the ∆G_water and ∆G_octanol is very important for the correct logP_ow. That is, the selection of the solvent in the calculation of the atomic charges is very important for the logP_ow calculation. Thus, the inappropriate selection of the solvents for the calculation of atomic charges results in the weak correlation between the logP_ow and the experimental data. This tendency can be observed in the logP_ow with P_logP{w,{w, w}}. The ∆G_water with P_water{w}, and ∆G_octanol with P_octanol{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_water with P_logP{w, {v, o}}.

The chemical group’s effect on the logP_ow values of a compound is shown in Figure 5. In the logP_ow calculation with P_logP{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_water value of ethyl phenyl ether was an overestimate. The calculated and the experimentally measured ∆G_water values were −10.22 kJ/mol and −17.92 kJ/mol, respectively. On the other hand, the calculated ∆G_octanol was an underestimate. The calculated and experimentally measured ∆G_octanol values were −34.91 kJ/mol and −23.65 kJ/mol, respectively. These differences caused large errors in logP_ow values. Similar errors can be found in the logP_ow values calculated with P_logP{w, {w, w}}. The profile of histograms for P_logP{w, {w, w}} is similar to that of histograms for P_logP{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_logP{w, {v, o}}. Therefore, the R values for P_logP{w, {w, w}} were slightly smaller than those for P_logP{w, {v, o}}.

On the other hand, most of the logP_ow values calculated with P_logP{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_water calculation and water environment was selected in ∆G_octanol calculations. These inadequate combinations induce large errors in the logP_ow 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_ow calculation.

To analyze the contributions of the atomic charges in different environments, the R and R² values of logP_ow in Table 4 were compared. In the calculated ∆G_water, the R and R² values become small in the order of P_water{v}, P_water{o}, and P_water{w} (Figure S3a). In terms of the atomic charges of the compounds in ∆G_octanol calculations, the R and R² values of logP_ow 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_ow values very slightly (Figure S3c). To briefly summarize these results, the atomic charges of the compound in the ∆G_water and ∆G_octanol 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_ow 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. LogP_ow for 17 Compounds

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

The calculated logP_ow values of the additional 17 compounds with the best combination of the atomic charges in water and octanol (P_logP{w, {v, w}}) showed a good agreement with the experimental data (Figure 6). The R and R² 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_logP{w, {v, o}} (Figure 4). These results indicate that the order of the compounds sorted by the logP_ow values, is very similar to the order of the compounds sorted by the logP_ow values measured experimentally. This means that the alchemical free energy calculation method with P_logP{w, {v, w}} will be effective to predict the lipophilicity of compounds, and will be a useful tool for drug design.

The logP_ow 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_ow values show strong correlation with the experimentally measured values, of which the R and R² 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 R² values are 0.85 and 0.74, respectively. These results reveal that the consideration of the atomic charge is important for improving the logP_ow values.

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_octanol{x, y}, in which x and y (for compound and octanol, respectively) are either v(vacuum), o(octanol), or w(water). The logP_ow 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_ow 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⁻¹ 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:

$$\mathsf{\Delta }G=\mathsf{\Delta }{G}_{solvent}-\mathsf{\Delta }{G}_{gas}$$

(2)

$${\mathrm{S}}_{gas}\stackrel{\hspace{1em}\hspace{1em}\mathsf{\Delta }{G}_{solvent}\hspace{1em}\hspace{1em}}{\to }{\mathrm{S}}_{solvent}$$

(3)

$${\mathrm{S}}_0\stackrel{\hspace{1em}\hspace{1em}\mathsf{\Delta }{G}_{gas}\hspace{1em}\hspace{1em}}{\to }{\mathrm{S}}_{gas}$$

(4)

where ΔG_solvent represents the free energy change from the compound in gas (S_gas) to the compound in solvent (S_solvent) in Formula (3), and ΔG_gas represents the free energy change from the compound in vacuum (S₀) to the compound in gas (S_gas) 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_i representing the potential energy of the i-th window, nonequilibrium work W (=U_i+1–U_i) 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

$$<{\left[1+\exp \left\{\left(W-\mathsf{\Delta }{G}_i\right)/RT\right\}\right]}^{-1}>{}_i=<{\left[1+\exp \left\{-\left(W-\mathsf{\Delta }{G}_i\right)/RT\right\}\right]}^{-1}>{}_{i+1}$$

(5)

where ΔG_i 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:

$$\mathsf{\Delta }G=\sum _i\mathsf{\Delta }{G}_i$$

(6)

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 LogP_ow

The logP_ow values were calculated using the following formula:

$$\begin{array}{ll}\log P_{ow}& =-\frac{\mathsf{\Delta }{G}_{transfer}^{octanol\to water}}{2.303 RT}\\ & =-\frac{\mathsf{\Delta }{G}_{octanol}-\mathsf{\Delta }{G}_{water}}{2.303 RT}\end{array}$$

(7)

where R and T are the gas constant and temperature, respectively, and ΔG_octanol and ΔG_water 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_ow value calculated from ΔG_water using P_water{x}, and ΔG_octanol using P_octanol{y, z}, was simplified to P_logP{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_ow 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_ow values were calculated by the combination of the free energies of the water and octanol solvents, i.e., three P_water{x} and nine P_octanol{x, y} parameters. Finally, 27 logP_ow values were calculated in total.

The logP_ow 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_ow 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_water{w} and P_octanol{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_ow 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_ow values was sensitive to the combination of the ΔG_water and ΔG_octanol values. The weak correlation was due to the different slopes of regression line between ΔG_water and ΔG_octanol. Finally, we found that the calculation with P_logP{w, {v, o}} provides the best correlations between the calculated logP_ow values and experimentally measured values.

We calculated the logP_ow 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_ow values were better than those obtained in previous studies. These results show that the solvent assumed in the charge calculation greatly influences the logP_ow values. The calculation of the atomic charges was time consuming, but dramatically increased the accuracy of the logP_ow values. We expect these accurate results to be useful in the field of drug discovery, and to accelerate drug design processes.