Replica-Exchange Enveloping Distribution Sampling Using Generalized AMBER Force-Field Topologies: Application to Relative Hydration Free-Energy Calculations for Large Sets of Molecules

Abstract

Free-energy differences between pairs of end-states can be estimated based on molecular dynamics (MD) simulations using standard pathway-dependent methods such as thermodynamic integration (TI), free-energy perturbation, or Bennett's acceptance ratio. Replica-exchange enveloping distribution sampling (RE-EDS), on the other hand, allows for the sampling of multiple end-states in a single simulation without the specification of any pathways. In this work, we use the RE-EDS method as implemented in GROMOS together with generalized AMBER force-field (GAFF) topologies, converted to a GROMOS-compatible format with a newly developed GROMOS++ program amber2gromos, to compute relative hydration free energies for a series of benzene derivatives. The results obtained with RE-EDS are compared to the experimental data as well as calculated values from the literature. In addition, the estimated free-energy differences in water and in vacuum are compared to values from TI calculations carried out with GROMACS. The hydration free energies obtained using RE-EDS for multiple molecules are found to be in good agreement with both the experimental data and the results calculated using other free-energy methods. While all considered free-energy methods delivered accurate results, the RE-EDS calculations required the least amount of total simulation time. This work serves as a validation for the use of GAFF topologies with the GROMOS simulation package and the RE-EDS approach. Furthermore, the performance of RE-EDS for a large set of 28 end-states is assessed with promising results.

Figure 1

Thermodynamic cycle to calculate relative hydration free energies ΔΔG_hyd for three small molecules i (aniline), j (1,3-dichlorobenzene), and k (anisole). Here, ΔG_hydⁱ is the hydration free energy of molecule i, ΔG_vac^ji is the free-energy difference between molecules i and j in vacuum, ΔG_wat^ji is the free-energy difference between the two molecules in water, and ΔΔG_hyd^ji is the hydration free-energy difference between the two molecules (relative hydration free energy). The free-energy difference between two molecules can be calculated from multiple pairwise simulations (as shown on the left, e.g., with TI) or from one simulation with multiple molecules (as shown on the right, e.g., with RE-EDS).

Figure 2

Set A consists of six benzene derivatives, selected from the FreeSolv^, database. The molecule indices, the FreeSolv identifiers, the SMILES strings, and the names of the molecules can be found in Table S1 in the Supporting Information.

Figure 3

Set B consists of 28 benzene derivatives, selected from the FreeSolv^, database. Set B was further subdivided into subset Ba (B1–B14) and subset Bb (B1 along with B15–B28). The molecule indices, the FreeSolv identifiers, the SMILES strings, and the names of the molecules can be found in Table S2 in the Supporting Information.

Figure 4

Schematic illustration of the RE-EDS input file preparation. The input files (topology, perturbed topology, coordinates, and distance restraints) for the RE-EDS simulations in GROMOS were created from the frcmod and mol2 files of the FreeSolv^, database. This workflow can easily be extended to also perform the molecule parametrization (i.e., to generate mol2 and frcmod files) using antechamber and parmchk.⁻

Figure 5

Potential-energy distributions of single-molecule simulations of set A based on 5 ns vacuum simulations in GROMOS (orange bars) and GROMACS (pink lines). The topologies for the simulations are based on the AMBER topologies taken from the FreeSolv^, database. The GROMACS topologies were translated with ParmEd, and the GROMOS topologies were converted with amber2gromos. Energies were written every 100 time steps (i.e., every 200 fs), and the first 1.25 ns of the simulations were discarded as equilibration.

Figure 6

Comparison of the relative hydration free energies of set A: RE-EDS (ΔΔG_hyd^RE-EDS) versus MBAR (ΔΔG_hyd^MBAR) (left) and RE-EDS versus experiment (ΔΔG_hyd^exp) (right) as reported by the FreeSolv^, database. The gray diagonal lines correspond to perfect alignment within ±4.184 kJ mol^–1 (±1 kcal mol^–1). The RE-EDS results were averaged over five independent production runs in vacuum/water, and the errors of the ΔG values correspond to the standard deviation over the five repeats. The error estimate of the ΔΔG values was calculated via Gaussian error propagation. The numerical values are provided in Table S3 in the Supporting Information. A plot of the deviations from experiment for the different methods is shown in Figure S7 in the Supporting Information. All pairwise comparisons between the different simulation methods and the experimental results are provided in Figure S5.

Figure 7

Comparison of the relative hydration free energies of set B: RE-EDS (ΔΔG_hyd^RE-EDS) versus MBAR (ΔΔG_hyd^MBAR) (left) and RE-EDS versus experiment (ΔΔG_hyd^exp) (right) as reported by the FreeSolv^, database. The gray diagonal lines correspond to perfect alignment within ±4.184 kJ mol^–1 (±1 kcal mol^–1). The ΔΔG_hyd^ji values are colored according to end-state i (i.e., the “reference molecule” for the calculation). The RE-EDS results were averaged over five independent production runs in vacuum/water, and the errors of the ΔG values correspond to the standard deviation over the five repeats. The error estimate of the ΔΔG values was calculated via Gaussian error propagation. The numerical values are provided in Table S4 in the Supporting Information. A plot of the deviations from experiment for the different methods is shown in Figure S12 in the Supporting Information. All pairwise comparisons between the different simulation methods and the experimental results are provided in Figure S9.

Figure 8

Comparison of the relative hydration free energies with RE-EDS of subsets Ba and Bb (ΔΔG_hyd^{RE-EDS,subsets}) compared to the full set B (ΔΔG_hyd^RE-EDS) (left) and compared to experiment (ΔΔG_hyd^exp) (right) as reported by the FreeSolv^, database. The gray diagonal lines correspond to perfect alignment within ±4.184 kJ mol^–1 (±1 kcal mol^–1). The ΔΔG_hyd^ji values are colored according to end-state i (i.e., the “reference molecule” for the calculation). The RE-EDS results were averaged over five independent production runs in vacuum/water, and the errors of the ΔG values correspond to the standard deviation over the five repeats. The error estimate of the ΔΔG values was calculated via Gaussian error propagation. The numerical values are provided in Table S4 in the Supporting Information. A plot of the deviations from experiment for the different methods is shown in Figure S12 in the Supporting Information. All pairwise comparisons between the different simulation methods and the experimental results are provided in Figure S9.