Hydration Free Energy from Orthogonal Space Random Walk and Polarizable Force Field.

The orthogonal space random walk (OSRW) method has shown enhanced sampling efficiency in free energy calculations from previous studies. In this study, the implementation of OSRW in accordance with the polarizable AMOEBA force field in TINKER molecular modeling software package is discussed and subsequently applied to the hydration free energy calculation of 20 small organic molecules, among which 15 are positively charged and five are neutral. The calculated hydration free energies of these molecules are compared with the results obtained from the Bennett acceptance ratio method using the same force field, and overall an excellent agreement is obtained. The convergence and the efficiency of the OSRW are also discussed and compared with BAR. Combining enhanced sampling techniques such as OSRW with polarizable force fields is very promising for achieving both accuracy and efficiency in general free energy calculations.

Introduction

Water is a substantial component of living organisms, forming an environment where biological processes such as the transportation of ions, the folding of proteins, and the activation/deactivation of signaling pathways can take place. The interactions between water and physiologically relevant molecules, such as monatomic ions, small molecules, and macromolecules, are crucial to the efforts of understanding such biological processes and applications such as protein engineering and drug discovery. Therefore, accurately modeling the hydration process is arguably the first step in modeling these biological processes and developing accurate physical models and robust computational approaches. For instance, not only is the hydration free energy (HFE) a key property in predicting the solubility of organic molecules and their binding to proteins,[1−4] but hydration free energy is also an important measure in the development and evaluation of the accuracy of force fields[5−10] and sampling methods.[11−16] The hydration free energy of a molecule can be calculated by using explicit solvent models, e.g., TIP3P water[17] and AMOEBA water,[18] in combination with alchemical approaches, such as thermodynamic integration (TI; see review by Kollman[1]), Bennett acceptance ratio (BAR),[19] or the orthogonal space random walk (OSRW).[4,20−23] Once the force field is well-defined, the accuracy and precision of the alchemical results are somewhat more predictable.

Although the importance of including explicit polarization in molecular modeling have been demonstrated in previous studies,[24−26] the routine application of polarizable force fields, such as AMOEBA,[18,27−30] to obtain accurate thermodynamic properties is still hindered by the high computational cost of traditional alchemical approaches. Thus, enhanced sampling methods such as the OSRW method are more appealing in such simulations. Unlike BAR or TI, which requires a number of arbitrary, fixed order parameter λ to connect the two end states, the OSRW method[4,20−23] described in later sections utilizes dynamic order parameters, λ, and dU/dλ, coupled with the metadynamics approach[31] to sample the two dimensions. In this way, the alchemical perturbation between the two end states can be performed in a single molecular dynamics simulation and with improved efficiency. Previously, we have demonstrated that OSRW allows efficient sampling of configurational spaces of molecular crystals.[4]

In this paper, the OSRW method is implemented with the polarizable multipole based AMOEBA force field in TINKER and applied to compute the hydration free energy of several small organic molecules. The hydration free energy results from OSRW are compared with those computed from the conventional BAR method, which has been utilized to compute free energy of hydration and binding in combination with AMOEBA in previous studies.[5,30,32−39] The results from the two approaches are in excellent agreement (RMSD = 0.49 kcal/mol), with the OSRW method showing a significant advantage in computational efficiency.

Methods

Theoretical Background of OSRW

Because free energy is a path-independent property, a common approach to calculate the free energy difference is to define a mixed potential, so that the potential functions of the two end states of interest can be connected analytically. Such a mixed potential is usually defined aswhere U is the Hamiltonian of the mixed potential, r is the coordinate, and the scaling parameter λ = 0 and 1 corresponds to the two end states, U0 and U1, respectively. The free energy change from one state to the other can thus be given aswhere G is the free energy of each state, ΔG is the change in free energy, and ⟨ ⟩λ is the ensemble average of each λ state.

Such construction of the approach, however, relies on the assumption that sufficient conformational sampling can be done as the system adjusts to the new intermediate states. For complex systems, such transition usually requires a significant amount of simulation time, especially when there are larger changes in structure. This is often known as the “Hamiltonian lagging” problem, which exists for methods where λ is a continuous and dynamic variable.[40] For methods that perform simulation at discrete λ “windows,” a large number of intermediate steps are required, and long simulations at each step are needed to ensure sufficient equilibration and sampling at each step.

Alternatively, by combining the ideas of the dynamic λ method[40,41] and the metadynamics method,[31] Yang and co-workers proposed an efficient free energy sampling approach, which they referred to as the orthogonal space random walk (OSRW).[4,20−23] In this approach, a random walk is performed in two dimensions, λ and its orthogonal generalized force Fλ = ∂U/∂λ. The use of dynamic λ itself unnecessarily improves the computational efficiency; however, directly biasing along the ∂U/∂λ dimension can potentially accelerate the free energy calculation since the integral of ∂U/∂λ is exactly the free energy. By repetitively adding a Gaussian-like repulsive potential to λ and Fλ spaces, the low energy wells can be “flooded” to overcome the energy barriers. The potential of the system in the OSRW can be written aswhere g is the biasing potential that can be defined recursively aswhere h and w are the height and width of the Gaussian, which can be adjusted to balance the accuracy and efficiency of the method, and t is the index of states. The free energy along the reaction coordinates can thus be estimated as –g(λ,Fλ). To move from an initial state to a target state, λ, the free energy change can thus be estimated asBy adaptively adding a negative G(λ) to the system Hamiltonian as shown in eq 3, the “flooding” of the free energy surface can be accelerated along the λ space.

AMOEBA Based Alchemical Scheme in TINKER

A hybrid Hamiltonian based on eq 1 is implemented to calculate the alchemical free energy using the AMOEBA force field in TINKER. A dual topology approach is used to keep the intramolecular energies of the mutating systems (e.g., solute molecules in solution) throughout the simulation, and the mutating systems A (initial) and B (final) never interact with each other. In case of absolute hydration free energy calculations performed in this study, the two topologies are solute-in-water and water without solute, respectively. In this approach, the mixed potential can be written aswhere the subscript dt indicates the Hamiltonian is for dual topology, and the bonded term is independent of λ. A softcore van der Waals (vdW) potential[4,32] has been adapted to replace the original buffered 14–7 potential[42] in AMOEBA in this implementation:where i and j are the indices of the atoms, ε is the well-depth of the potential, α is an adjustable constant, ρ = r/r*, and r* is the equilibrium distance between two atoms. This equation prevents the numerical instability of the system when λ is small, and it reduces to the original buffered 14–7 potential when λ = 1. The real space electrostatic potential[18] is written aswhere G is the permanent multipole moment, l is the order of the multipole moment, B is the screening function, and f is the modified distance defined as

This definition, similar to the softcore potential for the vdW term, can prevent the numerical instability of the system when λ is small and atoms are very close to each other. The final mixed potential for the real space electrostatic potential is thuswhere the superscripts tot and mut respectively indicate the energy of the whole system and the energy of the part undergoing alchemical transformation. The reciprocal space part of the permanent multipole interaction energy is mixed linearly:

The polarization energy has a slightly different combination. To eliminate the potential problem[4] of the self-consistent field calculation at unphysically small atomic distances, the polarization is switched off until λ is equal to or greater than 0.75, which is to referred as λpolstart. An artificial scaling factor, given as λpol,A = (1 – λpolstart – λ)/(1 – λpolstart) and λpol,B = (λ – λpolstart)/(1 – λpolstart), is used to smoothly switch the potential across 1 – λpolstart and λpolstart, respectively. The final mixed polarization potential can be written as

Various derivatives, ∂U/∂λ, ∂2U/∂λ2, and ∂2U/∂r∂λ have been derived in the same way as reported previously.[4]

Computational Details

In this study, the hydration free energy of 20 small molecule solutes,[43] among which 15 are positively charged compounds and five are neutral compounds, were calculated in an explicit solvent model[18] using the polarizable AMOEBA force field.[27] Parameters of the compounds were obtained using POLTYPE.[38] Both orthogonal space random walk and Bennett acceptance ratio methods were used to estimate the free energy using the same simulation conditions, such as box size, simulation ensemble, boundary conditions, vdW, and electrostatic cutoffs. Thus, the comparison between OSRW and BAR methods would not be affected by potential artifacts in the calculated hydration free energy due to different boundary conditions or treatment of electrostatic interactions suggested in previous studies.[44,45] All molecular dynamics simulations were conducted with a RESPA integrator,[46] Bussi thermostat,[47] and 2 fs time step using the TINKER software package if not otherwise stated. A cutoff of 12 Å was applied to vdW’s interaction with α = 0.07 (eq 7), while a cutoff of 7.0 Å was applied in the real space softcore electrostatic calculations with α = 2.0 (eqs 8 and 9). Self-consistent induced dipole moments were converged to below 0.00001 D per atom.

Before the alchemical simulations, all the solutes were first soaked in a 30 Å cubic water box followed by a 600 ps relaxation using NPT molecular dynamics simulation at 298 K with a 2 fs time step. The resulting boxes were used in the subsequent NVT simulations with the density fixed at the average from the NPT simulations.

In OSRW, a 2D grid along the λ and Fλ axes was constructed to store the history of the (λ, Fλ) states visited. Each point on the grid represents a bin with finite dimensions. Our implementation has a λ width of 0.005 and an Fλ width of 2.0. For the λ axis, λ ranges from 0 to 1, and for mathematical convenience the first and last λ bins are half size (the total number of the λ bins is 201) and centered at 0 and 1, respectively. For the Fλ axis, Fλ has no clear range, so the range must be dynamically updated if the calculated Fλ falls outside the initial specified range. Also for mathematical convenience, there is always an Fλ bin centered at zero. Throughout the MD simulation, each bin centered at (λ, Fλ,) in the grid represents the number of times a particular state with λ – Δλ/2 < λ < λ + Δλ/2 and Fλ – ΔFλ/2 < Fλ < Fλ + ΔFλ/2 was visited, where Δλ and ΔFλ represent the λ width and Fλ width, respectively. For each count, a 2D biasing Gaussian centered at (λ, Fλ) is added to the potential. Our implementation uses a Gaussian height of 0.005 with variances of w12 = (2Δλ)2 and w22 = (2ΔFλ)2. The Gaussians are cut off after five bins from the central bin. However, simply with the current implementation, the random walk may end up being stuck at the λ end points since λ is between 0 and 1. Thus, mirror conditions are enforced where if a Gaussian has a contributing value outside the λ range, the contribution is mirrored onto a bin with a valid range. For the first and last λ bins, the contributions are automatically doubled since they are centered at 0 and 1, respectively. The statistical error was estimated from three repeat simulations of 4 ns for each compound.

In BAR calculations, a three-step perturbation approach was applied using the AMOEBA force field.[5,30,32] To make the solute disappear in the solvent, the electrostatic and polarization interactions were perturbed in 11 windows, scaled by λ = 1.0, 0.9, 0.8, 0.7, 0.6, 0.5, 0.4, 0.3, 0.2, 0.1, 0.0. After the electrostatic interactions were scaled to zero, then the vdW’s contribution was perturbed using 14 windows with λ = 1.0, 0.9, 0.8, 0.75, 0.7, 0.65, 0.6, 0.55, 0.5, 0.4, 0.3, 0.2, 0.1, 0.0. Then, 1000 ps NVT simulations at 298 K were performed at each window. The recharging of each solute in the gas phase was modeled using 11 windows (with an interval of 0.1 for λ) of 1000 ps molecular dynamics simulation at 298 K, a stochastic integrator, and a 0.1 fs time step. Data collected from 50 to 1000 ps range was then processed using the BAR equations. The errors for BAR results were computed as a sum of errors from the individual alchemical perturbation steps.

Results and Discussion

Hydration Free Energy

The hydration free energy was calculated for the 20 compounds using both OSRW and BAR methods. In general, a good agreement between the two methods was obtained (Table 1). The energy values from the two methods were plotted against each other (Figure 1), and the r2 correlation coefficient was calculated to be 0.98 for the charged set and 0.99 for the neutral set. The unsigned average difference in the calculated hydration free energy between the two methods is 0.39 kcal/mol, and the root-mean-square difference is 0.49 kcal/mol (Table 1). For illustration, the hydration free energy over time is plotted for the OSRW (Figure 2) and BAR (Figure 3) methods for compounds 1, 11, and 17.

Table 1

Calculated Hydration Free Energy (kcal/mol) of 20 Compounds Using BAR and OSRW

–14.06 ± 0.22 kcal/mol after simulations extended from 1 to 5 ns at each of the 25 steps.

–12.97 ± 0.39 kcal/mol after simulations extended from 4 to 8 ns, averaged over seven independent simulations.

Data collected from 0.5 ns simulation in each step; otherwise 1 ns simulations were performed.

Figure 1

Comparison of the calculated hydration free energy from BAR and OSRW methods. An excellent agreement between the values from the two methods is obtained with an R2 coefficient of 0.98 ± 0.02 and 0.99 ± 0.01 for the 15 charged (upper panel) and the five neutral molecules (lower panel), respectively.

Figure 2

Plots of hydration free energy (kcal/mol) for compounds 1 (top), 11 (middle), and 17 (bottom) calculated by OSRW as a function of time.

Figure 3

Plots of hydration free energy of compounds 1 (top), 11 (middle), and 17 (bottom) calculated by BAR as a function of time.

Compound 17 shows the largest difference of 1.38 kcal/mol between the two methods, which reduced to 1.09 kcal/mol after we significantly extended simulations for both methods (Table 1). Compound 17 has a “complex” structure with two hydroxyl groups and one fluorine atom all connecting to adjacent carbons in the benzene ring. The interactions among these groups and with water can be rather complex. For example, from both BAR and OSRW simulations, the two hydroxyl groups were seen as very flexible with the hydrogen atoms either facing away or forming hydrogen bonds with each other, albeit with different frequencies. In addition, compound 17 is the only solute with a fluorine atom in this set. To verify the fluorine parameters, hydration free energy of fluorobenzene and 2-fluorophenol have also been calculated using OSRW and compared with experiment values. For each of the two compounds, three independent OSRW calculations were performed. The average hydration free energy from the simulations is −0.76 ± 0.33 and −5.71 ± 0.16 kcal/mol for fluorobenzene and 2-fluorophenol, respectively. This is in a good agreement with the experimental hydration free energy of −0.80 and −5.29 kcal/mol for the two compounds.[9]

We first examined the effect of van der Waals perturbation steps on the hydration free energy calculated from BAR for compound 17. In our implementation of the BAR approach, the electrostatic interaction between solute and environment was first turned off, and then the vdW interactions were scaled. During this latter stage, water and solute molecules can have significant overlap, resulting in large uncertainty in the free energy as evident by the difference between the forward and backward free energy perturbation results. We added an additional six windows in the middle of the vdW perturbation (λ = 0.775, 0.725, 0.675, 0.625, 0.575, 0.525). With these additional steps and longer simulations, the hydration free energy of compound 17 is merely increased by 0.2 kcal/mol (Supporting Information).

We have further investigated water structure near the solute sampled during the BAR simulations by plotting radial distribution functions (RDFs) between the compound 17 (O atoms in the two hydroxyl groups and F atom) and water (O). It would be best to compare the RDF at the same λ values. However, this is impossible in this study due to the limitation of our implementation. As explained in the Methods section, the BAR method decouples the vdW and electrostatic interactions separately while in the OSRW approach the two are scaled simultaneously. As a result, the same λ in the two methods actually represents different states. Therefore, we plotted the RDF for the three atom pairs using the trajectories of all the λ windows (Figure 4). Similar plots for OSRW simulations are shown in Figure 5, which naturally included all λ states visited during the simulations. The first apparent difference between the two methods is that the BAR RDFs have water peaks closer to the solute at around 1 Å. Note that, in addition to the difference in how the vdW and electrostatic solute–water interactions are scaled (sequential in BAR and simultaneous in OSRW), an equal number of coordinate structures were saved for each λ window in BAR simulations while the OSRW had an uneven distribution of λ values given its nature of importance sampling. Both factors would affect the overall RDF. According to these RDFs, there is much less structure in the water around the polar groups of the solute during OSRW simulations, which we attribute to the “flattened” energy surface by the biasing potential introduced in OSRW. Another interesting feature is that there is a consistent but faint “peak” at around 3 Å for the OSRW RDFs.

Figure 4

Plots of the radial distribution function (bin size 0.2 Å) from BAR simulations: (top) the oxygen of water and the fluorine of compound 17, (middle) oxygen of water and oxygen of compound 17 in the hydroxide group closest to the fluorine, and (bottom) oxygen of water and oxygen of compound 17 in the hydroxide group furthest from fluorine. The RDF was evaluated using the simulated structures from all perturbation windows.

Figure 5

Plots of the radial distribution function (bin size 0.2 Å) from OSRW simulations: (top) oxygen of water and fluorine of compound 17, (middle) oxygen of water and oxygen of compound 17 in the hydroxide group closest to fluorine in compound 17, and (bottom) oxygen of water and oxygen of compound 17 in the hydroxide group farthest away from fluorine in compound 17.

As seen from the time evolution of hydration free energy for compound 17 (Figure 2 bottom), the OSRW simulations actually first approached −15 kcal/mol in the beginning of the simulation but quickly increased to around −13 kcal/mol. This behavior was also observed in some other repeated OSRW simulations. The results suggest that BAR and OSRW may have sampled different phase space within the limited simulation time in this study. Thus, we have further extended the simulations for both methods. For BAR, we increased the simulation length at each window from 1 to 5 ns, while for OSRW we extended the simulation from 4 ns and 8 ns and added an additional four independent simulation runs for compound 17. The difference between the hydration free energy from the two methods dropped to 1.09 kcal/mol but was still greater than the statistical uncertainty. This indicates that, for molecules such as compound 17, substantially long simulations may be necessary to fully converge the answer.

Convergence and Efficiency of the Free Energy in OSRW and BAR

To examine the computational efficiency and convergence of OSRW and BAR methods, we analyze the simulations data for compounds 1, 11, and 17 as examples. The cumulative hydration free energies of these compounds calculated by OSRW as a function of time (Figure 2) clearly demonstrate its efficiency and convergence. For the neutral compound 17, the HFE has reasonably converged in about 300 ps in a single OSRW simulation. For the two charged compounds 1 and 11, it takes about 1.5 ns to reach a similar level of convergence while the absolute value of the HFE is much large for charged compounds.

Since the free energy is dependent on the λ and ∂U/∂λ, a sufficient sampling of the system in both spaces is crucial to the accuracy of the estimated free energy. Distributions of λ and its generalized force Fλ = ∂U/∂λ against time are shown in Figure 6. It can be seen that the OSRW method samples the order parameter space multiple times during the simulation. There are frequent peaks and valleys in the ∂U/∂λ plot. By examining the free energy surface with respect to ∂U/∂λ and λ in Figure 7, it can be seen that the large positive values of ∂U/∂λ correspond to the initial state when the solute is mostly decoupled from the environment while the negative attractive ∂U/∂λ occurs toward the end state when λ = 1. There is a broad and flat region in between, which corresponds to ∂U/∂λ around zero. A significant feature of the OSRW is that it flattens the energy surface in both λ and ∂U/∂λ dimensions, while the latter is the exact quantity needed to compute free energy (G = ∫⟨∂U/∂λ⟩/dλ). Most other methods that sample only the λ space do not actively handle the energy barriers along ∂U/∂λ at a given λ. Figure 7 shows that the surface along ∂U/∂λ at a given λ can be complicated. In a typical OSRW simulation, the system quickly reaches the minimum (negative) ∂U/∂λ region and biasing Gaussian potentials are added in both λ and ∂U/∂λ dimensions until it reaches the middle flat region (∂U/∂λ ∼ zero) and then glides back and forth along the λ dimension many times. Sufficient sampling near both initial and end states is important as both contribute significantly to the free energy.

Figure 6

λ and Fλ values over simulation time from the OSRW simulations of compounds 1, 11, and 17 (data collected at every 5 ps). The spikes in Fλ correspond to barriers crossed by the system.

Figure 7

Plots of free energy surfaces for compounds 1 (top), 11 (middle), and 17 (bottom). The color represents the number of visits to a particular state.

It should be noted that our implementation involves the use of metadynamics for the calculation of the biasing potential. The choice of biasing Gaussian heights and λ particle movement has not been fully optimized in this study. In addition, other implementations alternative to metadynamics can be used and have been explored.[23]

The convergence of the free energy calculated using BAR, as seen for compounds 1, 11, and 17 in Figure 3, generally takes ∼300 ps simulations for each window to converge to a reasonable level. However, if we take the number of windows needed in BAR simulation into account, it is roughly equivalent to 8 ns of a single molecular dynamics simulation. Thus, the superior efficiency of OSRW is apparent.

The OSRW simulation is in general easier to set up and process as compared with the BAR method. In the TINKER implementation, only a single molecular dynamics simulation is required, while the BAR approach requires the maintenance of ∼25 simulations (in this study). A way to improve the sampling ability of the BAR method is to allow configuration exchanges among different simulation windows as in the multistate BAR (mBAR) method.[48]

Conclusion

We have implemented the OSRW method for the polarizable multipole based AMOEBA force field in TINKER and applied it to the hydration free energy calculations. All pairwise interactions including vdW and electrostatics in real space were treated with soft-core. The dual topology approach is applied to all nonpairwise interactions such as the reciprocal component of Ewald, polarization energy, and forces. Our implementation involves the use of metadynamics to introduce the total biasing potential while other methods used to improve the overall robustness of the OSRW method are under consideration.[23] Current implementation in TINKER can be used to evaluate hydration as well as binding free energy of host–guest systems.

The OSRW method has been shown to exhibit superior efficiency for hydration free energy calculations which can be reliably calculated by using traditional methods such as BAR. The hydration free energies for a set of 20 compounds were computed by using both BAR and OSRW. We found excellent agreement between the two, with a RMSD of 0.49 kcal/mol. A noticeable difference between the OSRW and BAR results was observed for compound 17, which has two hydroxyl groups and a fluorine atom attached to a benzene ring, even after the simulations were extended significantly. This discrepancy shows that even for a relatively small solute, the sampling of solute–solvent configurations can be challenging. In the BAR method, 1 ns at each of 25 steps or more have been performed. For OSRW, one simulation of 4 ns led to similar results within the statistical uncertainty. While further studies are required to thoroughly evaluate the different free energy methods, our results suggest that OSRW is an efficient alternative sampling method that polarizable force fields can benefit from.