Accelerated adaptive integration method.

Conformational changes that occur upon ligand binding may be too slow to observe on the time scales routinely accessible using molecular dynamics simulations. The adaptive integration method (AIM) leverages the notion that when a ligand is either fully coupled or decoupled, according to λ, barrier heights may change, making some conformational transitions more accessible at certain λ values. AIM adaptively changes the value of λ in a single simulation so that conformations sampled at one value of λ seed the conformational space sampled at another λ value. Adapting the value of λ throughout a simulation, however, does not resolve issues in sampling when barriers remain high regardless of the λ value. In this work, we introduce a new method, called Accelerated AIM (AcclAIM), in which the potential energy function is flattened at intermediate values of λ, promoting the exploration of conformational space as the ligand is decoupled from its receptor. We show, with both a simple model system (Bromocyclohexane) and the more complex biomolecule Thrombin, that AcclAIM is a promising approach to overcome high barriers in the calculation of free energies, without the need for any statistical reweighting or additional processors.

Introduction

The calculation of free energies is critical in understanding many biological phenomena, for example in the molecular recognition between a small molecule and its receptor. Computational approaches for obtaining such free energies are of particular interest in drug design and in understanding fundamental enzymatic mechanisms.

Absolute binding free energies can be computed at low computational cost, for example by docking scoring functions, but generally lack accuracy and therefore cannot be directly compared to experimental values.[1] More accurate and computationally expensive approaches redefine the Hamiltonian of a molecular dynamics (MD) simulation as a function of both the atomic coordinates and a ligand coupling parameter λ, and, in this way, compute binding free energies using estimators like free energy perturbation (FEP), thermodynamic integration (TI), or the multistate Bennett acceptance ratio (MBAR).[2−4] Even these computationally costly efforts, however, can yield grossly inaccurate free energies due to inadequate sampling of the phase space of the system.[5,6]

A clever way to avoid inadequate sampling is to introduce external forces that guide sampling during the free energy calculation. This has been done, for example, by coupling TI with umbrella sampling, or in the Confine-and-Release method, in which binding free energy calculations are performed independently for discrete conformational minima known to interchange slowly.[7,8]

One way to facilitate conformational transitions without knowledge of the relevant conformational minima a priori is to simply add a bias that effectively raises energy minima of the Hamiltonian.[9−11] With the added bias known, the equilibrium ensemble average can be recovered by reweighting. In our group, we use a variation of this approach called accelerated MD (aMD), which has been used as a means to more rapidly explore conformational space of biomolecules and improve free energy calculations.[12−15] Modifying the underlying Hamiltonian in complex systems, however, usually results in a great deal of statistical noise upon reweighting, making it impractical to recover equilibrium ensemble averages.[16,17]

Numerous other free energy methods utilize a replica-exchange framework in which multiple replicas, each with different values of λ or otherwise modified Hamiltonians, are simulated in parallel and exchange their configurations periodically, to improve convergence of free energy calculations.[18−22] In such methods, the number of replicas scales as the square root of the total number of degrees of freedom, again hindering the applicability of such approaches to complex biomolecular systems, although recent efforts have shown that specifically tempering the solute degrees of freedom can ameliorate this problem.[23,24]

Similar to a replica-exchange framework, the adaptive integration method (AIM) enhances mixing in λ space by adaptively changing λ throughout a single simulation.[25] By introducing an additional biasing factor that improves the likelihood of transitions between λ values, AIM promotes coverage of λ space and allows conformations explored at one value of λ to inform the conformational sampling at other values of λ. The AIM approach, however, does not prevent the system from being trapped by high barriers that may persist regardless of the value of λ.

Here we show that by introducing a modified Hamiltonian only at λ values between the λ = 0 and λ = 1 end points, we are able to accelerate the exploration of phase space as well as of λ space, without the need for additional CPU resources or any statistical reweighting. We evaluate this approach, which we call Accelerated AIM (AcclAIM), using three alchemical transformations: (i) a symmetric self-transformation of Bromocyclohexane (BRC), (ii) an asymmetric transformation of charge removal from BRC, and (iii) a calculation of the relative binding free energy of two ligands that bind to Thrombin.

Theory and Methods

Adaptive Integration Method

The Adaptive Integration Method (AIM) can help overcome some of the convergence problems encountered with traditional sampling protocols in free energy calculations. The sampling protocols we discuss here could in principle be combined with any of a number of robust free energy estimators; in this work, we specifically use thermodynamic integration (TI)[25] to obtain free energies.

With TI, the free energy is calculated using[26,27]where ΔF is the change in free energy, U is the potential energy, and λ is a parameter which links the initial state, where λ = 0, to the final state, where λ = 1.

The sampling protocol commonly used in free energy calculations consists of running multiple MD simulations at fixed λ values and then numerical integration is applied to eq 1. Two commonly used methods of numerical integration are trapezoidal and cubic spline interpolation, where the latter has been shown to improve results.[28]

Running multiple simulations at fixed λ values can be inefficient because even though some simulations may cover extensive regions of phase space at particular values of λ, other simulations may remain trapped in their local minima.[25]

AIM improves sampling by changing λ throughout a single simulation, which allows conformations sampled at one λ to seed sampling at other values of λ.[25] This is similar to λ-dynamics,[29] where λ is propagated as the coordinate of a fictitious particle, changing its value during a single simulation. However, instead of propagating λ as the coordinate of this particle, AIM uses a periodic Monte Carlo (MC) step to change in λ during the simulation. For a trial move from λ to λ, the acceptance probability is given bywhere ΔU = U(λ) – U(λ) is the difference in the potential energy, and ΔF = ΔF is a biasing factor which can help the simulation sample λ values more efficiently by overcoming barriers in U(λ).[25] The efficiency of moving between λ values thus depends on the difference ΔE between the potential energy and the biasing factor given by

We explore two biasing factors, denoted as MC1 and MC2. MC1 is the biasing factor originally used with AIM,[25] which uses the free energy as calculated with TIThe brackets indicate averages over samples at a given λ value.[25] In practice, discrete values are used for λ and ∂U/∂λ values are stored for each λ value so that numerical integration can be used to evaluate eq 4.

MC2 is a modified biasing factor which approximates the free energy between λ and λ + 1 usingwhere N is the number of attempted moves between the adjacent λ values. This approximation assumes that the system essentially stays in equilibrium at adjacent λ values.[30] If this assumption does not hold for a particular system, then at long sampling times the MC2 biasing factor will not sample from a flat distribution, which may reduce the efficiency of moving between λ values.

The motivation for defining an alternative biasing factor was to assess whether considering the difference in the potential energy could improve the rate at which the system accepts transitions between different λ values. Importantly, both eqs 4 and 5 follow the relationwhich is necessary for the MC step to satisfy balance. Note that while the biasing factors given by eqs 4 and 5 change at each step, they do converge given enough sampling time. Although the addition of this biasing term does not satisfy detailed balance, sampling will converge to the Boltzmann distribution asymptotically.[25,31,32]

Accelerated Adaptive Integration Method

Dynamically changing the value of λ during the simulation augments the overlap of phase space among all values of λ. Such an approach effectively resolves any inadequacy in sampling, provided that all of the relevant conformational minima are accessible in at least some of the values of λ visited. As illustrated in Figure 1a, however, the AIM method can not resolve issues in sampling in which high barriers separate conformational minima regardless of the value of λ.

Figure 1

Example potential energy surface showing the difference between (a) the Adaptive Integration Method (AIM) and (b) the Accelerated Adaptive Integration Method (AcclAIM). When using AIM, barriers in the potential energy surface that are present at all λ values can hinder conformational sampling. When using AcclAIM, barriers in the potential energy are lowered at intermediate λ values by the scaling factor γ, which allows the system to explore more conformations.

To facilitate sampling when such large barriers are present at all λ values, the potential energy can be modified according towhere α is an adjustable parameter that controls the maximum level of scaling and γ is the scaling factor. Then ∂U/∂λ is modified accordingly

Importantly, at λ equal to zero or one, the scaling factor γ is equal to one, resulting in the original unbiased potential. At intermediate values of λ, however, the scaling factor is between α and 1, decreasing the potential energy and the size of barriers, as illustrated in Figure 1b.

The idea of modifying the potential energy to improve conformational sampling has been explored previously with Hamiltonian replica exchange (HREX)[23,33] and the method of hidden restraints.[34] However, this is the first time that this idea has been combined with AIM, alleviating the requirements for multiple replicas (needed with HREX) and predefined target states (needed with the hidden restraints method). This specific construction of a λ-dependent biasing approach results in no bias applied at the end points, λ = 0 and λ = 1, such that statistical reweighting is avoided entirely, in contrast to methods like accelerated MD.

Simulation Details

The pmemd module in Amber 12 has been modified to include alchemical transformation calculations.[27] We have added AIM and AcclAIM to this previously modified version of pmemd. Similar parameters were used for the model systems and any model-specific parameters are detailed below.

System Parametrization

Small molecules were parametrized using the generalized Amber force field (GAFF),[35] and the protein was modeled using the Amber 99SB-ILDN force field.[36,37] All simulations included water molecules, which were modeled using TIP3P.[38] Atomic charges for the small molecules were calculated using RESP[39] to fit the electrostatic potentials determined using Gaussian03[40] at the Hartree–Fock/6-31G* level of theory. Missing hydrogens were added using LEaP.[41] A truncated octahedral box was used with a minimum distance of 10 Å to the solute. Particle mesh Ewald (PME)[42] with a 1 Å grid was used to calculate long-range electrostatics. An 8 Å cutoff for short-range nonbonded interactions was used. A 2 fs time step was used, with the SHAKE[43] and SETTLE[44] algorithms used to constrain bonds to hydrogen atoms.

Equilibration Protocol

Equilibrated systems for all the simulations were prepared using a similar protocol. First, a 20000-step steepest descent minimization was run. Then the systems were heated to 300 K over 500 ps using a Langevin thermostat[45,46] with a 2 ps–1 collision frequency. The density was adjusted using a Berendsen barostat,[47] with the pressure set to 1 bar and a 2 ps coupling constant. Finally, the system was equilibrated for 500 ps at constant temperature and pressure.

Distribution of λ Values

The following protocol was used to determine the spacing of λ values for the AIM and AcclAIM simulations. First, λ was set to 0 and the equilibrated structures were used as a starting point for both AIM and AcclAIM simulations. A short 110 ps simulation was run starting with 11 equally spaced λ values. Every 50 steps, the biasing factor ΔF (either MC1 or MC2, where relevant) was recorded, along with the difference ΔE (eq 3) between the potential energy and biasing factor, which is related to the probability of accepting an attempted MC move byEvery 5 ps, λ was changed to its adjacent value. This continued until λ = 1 at which point λ was decreased until λ = 0. The average ΔE value describes the difficulty in moving from one λ to the next. If ΔE for any pair of adjacent λ windows exceeded the cutoff of 2 (kcal/mol), N additional values were added linearly, where N is given byThe number of λ values resulting from this method are summarized in Tables S4 and S5 in the Supporting Information.

Production Simulations

Long 55 ns simulations were run with the new set of λ windows, starting from the equilibrated structure. For the BRC transformations, changes in λ were attempted every 50 steps, while for the Thrombin ligand transformation changes in λ were attempted every 500 steps. The acceptance probability is given by eq 2. Energies and coordinates were saved at the same frequency for later analysis.

Acceleration Parameters

The use of scaling in AcclAIM is generalizable to any combination of the components of the potential energy and can be limited to designated regions of the system. Some components of the potential energy that are of particular relevance in binding free energy calculations include the dihedrals (D), the 1–4 van der Waals and 1–4 electrostatic interactions (4), which are the interactions between atoms connected by three bonds, as well as the short-range nonbonded van der Waals and electrostatic interactions (N), which are between atoms not bonded to each other but within the short-range cutoff. Scaling the bond or angle terms, for instance, would not lead to improved conformational sampling, but instead would distort the geometry of the molecules.

In both cases of BRC and Thrombin examined here, the D, 4, and N components of the potential energy were scaled, as these terms are relevant in the problematic chair-flipping and ring-flipping conformational changes in these systems. We did, however, verify that various other combinations of scaled D, 4, or N terms would also be adequate for a simple model system like BRC (Supporting Information Figures S3–S6).

In the BRC transformations, all solute atoms were scaled, whereas only a select group of relevant atoms in the Thrombin ligand transformation were scaled. Scaling the entire potential energy surface of a biomolecule like Thrombin would reduce the overlap in the potential energy between adjacent λ values, decreasing the exchange probability and thus hindering convergence of these simulations.

For the BRC transformations, the scaling term α was set to zero for the D and 4 terms and was set to 0.1 for the N terms. For the Thrombin ligand transformation, α was set to 0.168 16 for all terms, which was used in previous work that examined this transformation using a HREX method.[24] Although we did test both more and less aggressive acceleration parameters (α = 0.05 or 0.25, respectively) for the more complicated case of Thrombin ligand binding (Supporting Information Figures S7 and S8), we found no dependence of the free energies on the choice of scaling parameter in this range. We note that in all cases we avoided scaling the nonbonded terms to zero, to prevent issues that occur when linearly scaling the potential for nonbonded terms. These issues are similar to those that led to the development of softcore potentials.[48−50] In principle, softcore scaling could be used instead of the linear scaling that we used, but this scheme was simpler for the initial implementation. A suggested protocol for choosing an appropriate value for α is included in the Results and Discussion section.

To specifically address the effect of introducing the scaling function, we compared results from AcclAIM simulations to those obtained from AIM simulations, in which no component of the potential energy is scaled.

Biasing Factor

In all of the production simulations, we used the alternative biasing factor, MC2, given by eq 5 at the MC steps. We separately evaluated the impact of the biasing factor by running AIM simulations of the BRC transformations with MC1 and MC2 (Supporting Information Figures S3–S6). We also examined how the biasing factor affects the number of full round-trips (i.e., λ = 0 → λ = 1 → λ = 0) over time in AIM and AcclAIM solvation free energy simulations.

Alchemical Transformations

Bromocyclohexane Self Transformation

The first alchemical transformation we considered was the self-transformation of Bromocyclohexane (BRC). In this transformation, both λ = 0 and λ = 1 describe BRC, but they differ in the position of the Br substituent, such that the equatorial or axial conformations, separated by a high barrier, were coupled at λ = 0 and λ = 1, respectively (see Figure 2a). The free energy of the transformation is 0 (kcal/mol), since this is a transformation to the same molecule, but this result depends on the accurate sampling of both conformations at both end points. The AcclAIM simulations scaled the parts of the potential involving the BRC molecule.

Figure 2

Alchemical transformations studied with AcclAIM. (a) Bromocyclohexane self-transformation, where the carbon atoms are offset in the λ = 0 and λ = 1 states. The transparent hydrogen and bromine atoms interact fully at λ = 0 and the solid hydrogen and bromine atoms interact fully at λ = 1 such that λ = 0 corresponds to the equatorial chair conformation and λ = 1 corresponds to the axial chair conformation. (b) Removal of charges from BRC. Two sets of simulations were run, starting from either the equatorial (top) or axial (bottom) conformations. (c) Thrombin (green ribbon) bound to CDA and CDB. CDB has an additional methyl group shown in red. The P1 pyridine ring, which is highlighted in yellow, can favor either the In conformation (shown) or the Out conformation depending on its substitution. For this model, the relative binding free energy between CDA and CDB was calculated. Two sets of simulations were run, starting from either the In or Out conformations. Images were rendered using Tachyon[51] in visual molecular dynamics (VMD).[52]

A modified softcore potential was used for the perturbed H and Br atoms. As implemented in Amber, softcore transformations decouple the energy of the perturbed part of the molecule.[27,48−50] However, the goal of this calculation is to determine how the free energy relates to conformational sampling of BRC. In order to approximate the change in the free energy due to the conformational change in BRC from λ = 0 to λ = 1, the dihedral and 1–4 terms were scaled linearly and were not decoupled from the rest of the system. This modification was tested by comparing the resulting free energy to the result from an umbrella sampling (US) calculation as described below.

Preparation of the equilibrated structures for US was similar to that described above, except a restraining potential was enabled. Production simulations were run for 5 ns at constant temperature and volume. Energies and coordinates were recorded every 50 steps. The reaction coordinate ϵ for this system is given by[53]where d1 is the dihedral between C1, C2, C3, and C4; d2 is the dihedral between C2, C3, C4, and C5; and so on. This type of restraining potential was added to pmemd so that US calculations could be carried out. Thirty-three windows were used, equally spaced between −80° and 80°. The force constant for this calculation was 131 (kcal/mol·rad2).

Bromocyclohexane Charge Removal

To assess AcclAIM in an asymmetric transformation that results in a nonzero change in free energy, we considered the free energy associated with the removal of charges from BRC. To probe convergence, we compared the results from independent simulations that were started from either axial or equatorial conformations (see Figure 2b), as simulations starting from different initial states should converge to the same free energy value in the limit of infinite sampling. As in the previous transformation, the AcclAIM simulations only scaled the parts of the potential involving BRC. Additional US calculations were run using the same protocol as described above, except that all of the charges on BRC were removed, to obtain appropriate reference values.

Thrombin Ligand Transformation

The final test transformation considered was the relative binding free energy to Thrombin for the ligands 2-(6-chloro-3-[2,2-difluoro-2(2-pyridinyl)ethyl]amino-2-oxo-1(2H)-pyrazinyl)-N-[(2-fluoro-6-pyridinyl)methyl]acetamide and 2-(6-chloro-3-[2,2-difluoro-2(2-pyridinyl)ethyl]amino-2-oxo-1(2H)-pyrazinyl)-N-[(2-fluoro-3-methyl-6-pyridinyl)methyl]acetamide, abbreviated as CDA or CDB, respectively. This transformation has been studied previously using replica exchange with solute tempering (REST) an HREX method.[24] The only difference between CDA and CDB is the addition of a methyl group on the P1 pyridine ring (See Figure 2c). This causes CDB to favor the conformation where the P1 pyridine ring is pointing in that is, where the methyl group is pointing in toward the protein. As with the previous study, AcclAIM was applied to the P1 pyridine ring and the conformational sampling of this ring was examined. As with the BRC charge transformation, separate simulations were started with the P1 pyridine ring in the In and Out states to assess convergence. Also, as with the previous study, the F atom in the P1 pyridine ring was set to be softcore to reduce the effective size of the ring and enhance sampling.[24] This modification was applied to all simulations, including the AIM simulations.

Initial coordinates came from the crystal structure for Thrombin bound to CDA (PDB 1MU6) and the structure for CDB came from PDB 1MU8.[54] The loop comprising of residues 146 to 150 near the ligand binding site was modeled from a more complete structure (PDB 4MLF).[55] The solvation free energy (SFE) was calculated by running AIM and AcclAIM simulations for the ligand transformation in solution. The relative binding free energy was calculated using

Methods of Analysis and Error

We performed production simulations using the US, AIM and AcclAIM simulations four times each. The potential of mean force (pmf) was calculated from the US simulations using the weighted histogram analysis method (WHAM)[56] as implemented in WHAM version 2.0.7 by Alan Grossfield (http://membrane.urmc.rochester.edu/content/wham/). We averaged the pmf across the four trials. The error for each trial was determined using a MC bootstrap analysis,[56] and the error was propagated according to

Free energies for AIM and AcclAIM were calculated using cubic spline interpolation with subsampling as implemented in pyMBAR (https://simtk.org/home/pymbar).[28,57] For the calculation of free energy as a function of simulation time, the free energy was calculated over all ∂U/∂λ values collected up to that point using pyMBAR. The free energy was averaged across the four simulations and the error was propagated according to eq 13.

The probability of finding the model system in each of its major conformations at λ = 0 or 1 was also calculated as a function of simulation time. For BRC, the conformation was determined according to ϵ defined in eq 11. The cutoff values for each conformation were determined from the location of the barriers in the pmf (see Supporting Information Figure S2). For Thrombin, the dihedral angle of the ligand atoms N6, C15, C16, and C17 was used to classify the In and Out states, where states above 0° were classified as In and states below 0° were classified as Out. The conformational probabilities were averaged over the four trials, and the standard deviation was used as an estimate of the error.

Results and Discussion

BRC Self-Transformation

The BRC self-transformation is simple, as the conformational landscape of the system is easily characterizable in one dimension by the ε coordinate, yet it is an illustrative example of the problems encountered in free energy calculations when the conformational space contains kinetically trapped states. The reference probabilities for the equatorial and axial conformations were obtained from the potential of mean force calculations, using the ε coordinate, as shown in Supporting Information Figure S2a. The population of the equatorial conformation was calculated as a function of simulation time at λ = 0 and λ = 1. The underlying potential is symmetric with respect to λ, so the populations at λ = 0 and λ = 1 should converge to the same value.

The population probabilities are summarized in Figure 3a. For the AcclAIM simulations, the populations at both end points converged, within error, to the reference population. The populations from the AIM simulations, on the other hand, in which no scaling was applied, did not overlap between the two end points nor with the reference value. Because a high ≈10 (kcal/mol) barrier is associated with chair flipping, the molecule was essentially trapped in the same chair conformation, regardless of the value of λ. Since the Br was fully coupled in an axial orientation at λ = 1 and fully coupled in an equatorial orientation at λ = 0, the resulting populations were inaccurate when scaling was not applied (AIM).

Figure 3

Results for the BRC self-transformation. Values represent averages over four independent trials. (a) Population of equatorial conformation sampled as a function of simulation time. The horizontal green line indicates the reference population as determined using US. Populations are shown at both λ = 0 (purple series) and λ = 1 (red series). Due to the symmetry in this self-transformation, the populations at both λ values should converge to the same value. Error bars represent the standard deviation over these trials. (b) Free energy as a function of simulation time. The free energy for a self-transformation should converge to 0 (kcal/mol), as indicated by the horizontal green line. Error bars represent the uncertainty, which was propagated from each trial using eq 13.

As expected, the disparities in the configurational sampling at the end points observed with conventional AIM translated into inaccurate, i.e. nonzero, free energies for the BRC self-transformation, whereas the free energies obtained with AcclAIM converged close to zero (Figure 3b). For the AIM simulations, the resulting free energy was −0.8 to −0.9 (kcal/mol), which is consistent with the difference in the free energy between the equatorial and axial conformations, according to the calculated pmf (Supporting Information Figure S2a). This agreement between US and AIM calculations not only highlights the inadequate sampling but also shows that the alternative softcore potential implemented for these simulations is valid for the determination of the free energy difference between different conformations. The BRC self-transformation demonstrates clearly how the inclusion of the scaled potential at intermediate λ values in the AIM framework can improve sampling of kinetically trapped conformations and correspondingly result in accurate free energies in cases where the AIM method alone is insufficient. A nearly identical BRC self-transformation, though without any partial atomic charges, has been used in previous work,[21] in which Hamiltonian replica exchange was used to enhance sampling. The charges included in the present work both improve the representation of BRC in water and are expected to make the test transformation more challenging. Nevertheless, we found that the AcclAIM method achieved comparable accuracy in the calculated free energy as the replica exchange approach. The latter approach requires multiple replicas at each λ value whereas the AIM framework permits the calculation of free energy from a single simulation that adaptively samples different Hamiltonians and λ values.

BRC Charge Removal

In addition to the self-transformation of BRC, we considered the asymmetric transformation, with an associated nonzero change in free energy, of removing the charges from BRC. Because equilibrium properties are independent of the initial state of the system, the conformational populations and the final free energy calculated from simulations starting from different initial conformations should converge to the same value, and we consider this overlap as a metric for convergence. The population of the equatorial conformation, calculated at λ = 0, is shown in Figure 4a (corresponding data from λ = 1 is provided in Supporting Information Figure S5). At both end points, as observed in the symmetric self-transformation, the equilibrium population was correctly sampled with AcclAIM, but the BRC molecule was confined to its initial conformation with the conventional AIM approach.

Figure 4

Results for the BRC charge removal simulations. Values are shown for simulations started from the equatorial (purple series) or axial (red series) conformations, which should converge to the same value. (a) Population of equatorial conformation at λ = 0 as a function of time, from simulations with either initial conformation. (b) Free energy as a function of time, from simulations with either initial conformation.

The free energy as a function of simulation time is shown in Figure 4b. For the AcclAIM simulations, the free energy values calculated using simulations starting from the two different conformations converged, within the estimated error, whereas the results from the AIM simulations depended on their initial conformations, demonstrating the impact of the kinetically trapped states. These results show that AcclAIM improves conformational sampling, which, in turn, can improve convergence of the free energy for alchemical transformations.

Thrombin Ligand Transformation

The final alchemical transformation, determining the relative binding free energy of the ligands CDA and CDB bound to Thrombin, involves a more complex system with substantially more degrees of freedom than BRC. On the basis of previous work,[24] it was assumed that sampling the two conformations of the P1 pyridine ring was critical in obtaining an accurate relative binding free energy. The population of the P1 pyridine ring in the In conformation was used to determine convergence of conformational sampling. The results for λ = 0 are shown in Figure 5a and for λ = 1 in Supporting Information Figure S7. These results show that the AIM simulations do not sample both conformations of the P1 pyridine ring, because they stay in their initial conformations. With AcclAIM, we were able to obtain convergence in the sampled populations between simulations, regardless of their different initial conformations.

Figure 5

Results for the transformation of CDA to CDB while bound to Thrombin. Values are shown for simulations started with the ligand P1 pyridine ring starting in the In conformation (purple series) or in the Out conformation (red series). (a) Population of the P1 pyridine ring in the In conformation at λ = 0 as a function of time. (b) Relative binding free energy as a function of simulation time, calculated using eq 12.

The corresponding free energy as a function of time is shown in Figure 5b. With AIM, there is a difference of 2 (kcal/mol) in the relative binding free energy, depending on the initial conformation. With AcclAIM, the free energies converge rapidly (within 25 ns). Interestingly, the free energies obtained here are in agreement with the experimental value reported for this transformation and are quite close to the free energy obtained using replica exchange with solute tempering (REST), despite the difference in force fields used.[24] If the relevant conformations are unknown a priori, we find that the computational efficiency of this type of free energy calculation is greater with AcclAIM than with REST (Supporting Information Figure S7, α = 0.168 16, compared to Figure 4b of ref (24)); however, the use of multiple initial configurations in a REST calculation would diminish such a difference.

We found that the convergence of the free energy was independent of the parameter α, within a tested range of 0.05–0.25 (Supporting Information Figure S8). Interestingly, we observed this marked improvement in the relative binding free energy in spite of a fairly large variance in the fraction of In population sampled. These results show that for the relative binding transformation of CDA and CDB to Thrombin, even a somewhat modest enhancement in the exploration of configurational space can translate to dramatic improvements in the accuracy of the free energy, and we anticipate this to be true for other biomolecular systems as well.

When applying AcclAIM to other systems, we would recommend starting with a lower scaling factor such as α = 0.05 and scaling the D, 4, and N terms. Then, the protocol described in Distribution of λ Values can be used to determine the number of λ values. If the number of windows is too large then the protocol should be repeated after increasing the scaling factor by 0.05. If additional information about the barriers in the system is known, selective scaling of the D, 4, or N terms could also be used to reduce the number of windows, as was done with the BRC transformations (Supporting Information Figures S3–S6). Based on these results, we would recommend a maximum number of λ values to be between 13 and 15. A larger number of λ values would indicate that the perturbation is too large, which will hinder convergence of the free energy.

Impact of the Biasing Factor

We introduced an alternative biasing factor MC2 (eq 5) for the MC step in AcclAIM. The motivation for modifying the biasing factor was that the difference in free energy between adjacent values of λ is expected to be more pronounced when the Hamiltonians are differentially scaled according to λ. Initial estimates of the free energy difference using TI, as with MC1 (eq 4), can be inaccurate, which causes a large difference between the potential energy and the biasing factor ΔU – ΔF. This makes moves between adjacent values of λ more difficult, decreasing the total number of λ round-trips. By estimating the free energy difference directly from the potential energy, as with MC2, the discrepancy between the potential energy and the biasing factor is reduced, which improves the probability of accepting MC moves for short sampling times. At longer sampling times, the efficiency of switching between λ values will depend on how well the system stays in equilibrium, following from the approximation inherent to the alternative biasing factor.

To evaluate the impact of the biasing factor on the number of round-trips in λ space, we performed AIM and AcclAIM simulations of the relative solvation free energy of CDA to CDB with the original and alternative biasing factors. As is shown in Figure 6a, the number of round-trips in λ space was not affected by the choice of biasing factor in conventional AIM simulations. However, as shown in Figure 6b, with the AcclAIM simulations we found a significant increase in the number of round-trips in λ space due as a result of the alternative biasing factor. For these simulations, the λ values were determined for MC2, and the same values were used for MC1 so that the only difference was the biasing factor used in the MC step. To verify that the gain in number of round-trips was not simply due to the distribution of λ values, additional simulations were run, in which MC1 was used in the protocol for determining the number of λ values. With these simulations, the total number of λ values increased, and the total number of λ round-trips decreased, suggesting that using the λ values determined with MC2 was more appropriate. It appears that for this transformation, the approximation underlying the alternative biasing factor was reasonable, because we did not see a decrease in the number of λ round-trips at longer times. Thus, for cases in which this approximation holds, the alternative biasing factor introduced here improves the effectiveness of scaling the potential at intermediate values of λ by improving the resulting number of λ round-trips. For cases in which this approximation does not hold, moves between adjacent values of λ could become more difficult, and thus for these systems we would recommend using the original biasing factor.

Figure 6

Number of λ round-trips as a function of time for the solvation free energy of the ligands that bind to Thrombin. Shown here are the results for (a) AIM and (b) AcclAIM simulations run using the original (MC1) or new biasing factors (MC2). The number of round-trips was averaged over four independent trials. The error was estimated by taking the standard deviation over these trials.

Conclusions

Using both a simple model system, Bromocyclohexane, as well as more complex biomolecule, Thrombin, we have demonstrated that the introduction of a scaled potential at intermediate λ values to the AIM method results in improved sampling of relevant conformational states at the end points of an alchemical transformation. Importantly, we have shown that for the simple case of BRC, the AcclAIM method is effective for both symmetric and asymmetric alchemical transformations. The challenge in using AcclAIM to study binding to complex biomolecules, however, remains being able to identify the slow (and most relevant) degrees of freedom that should be accelerated in the system.

We have also considered an alternative biasing method in AIM and AcclAIM that improves the likelihood of transitions in λ space when intermediate λ values are being scaled (e.g., in AcclAIM), and is still valid in the absence of scaling (e.g., in AIM). In practice, the mixing in λ space, e.g. the frequency of round-trips made in an AcclAIM simulation, can be conveniently used to gauge the level of acceleration applied to the chosen atoms in the system. Furthermore, the ability to specifically accelerate a ligand or residue with slowly exchanging conformational states, without having to identify an optimal physical reaction coordinate, is a clear strength of the AcclAIM approach. Because the potential energy at the end points is unbiased, the resulting free energies require no statistical reweighting, although the issue of poor exchanges between λ values still makes overacceleration a relevant issue.

Future work may consider alternative methods of determining λ values as well as the degrees of freedom that would benefit from enhanced sampling. Future consideration of alternative forms for the scaling factor may also optimize transitions in the presence of scaled potentials.