Alchemical Binding Free Energy Calculations in AMBER20: Advances and Best Practices for Drug Discovery

Abstract

Predicting protein-ligand binding affinities and the associated thermodynamics of biomolecular recognition is a primary objective of structure-based drug design. Alchemical free energy simulations offer a highly accurate and computationally efficient route to achieving this goal. While the AMBER molecular dynamics package has successfully been used for alchemical free energy simulations in academic research groups for decades, widespread impact in industrial drug discovery settings has been minimal because of the previous limitations within the AMBER alchemical code, coupled with challenges in system setup and postprocessing workflows. Through a close academia-industry collaboration we have addressed many of the previous limitations with an aim to improve accuracy, efficiency, and robustness of alchemical binding free energy simulations in industrial drug discovery applications. Here, we highlight some of the recent advances in AMBER20 with a focus on alchemical binding free energy (BFE) calculations, which are less computationally intensive than alternative binding free energy methods where full binding/unbinding paths are explored. In addition to scientific and technical advances in AMBER20, we also describe the essential practical aspects associated with running relative alchemical BFE calculations, along with recommendations for best practices, highlighting the importance not only of the alchemical simulation code but also the auxiliary functionalities and expertise required to obtain accurate and reliable results. This work is intended to provide a contemporary overview of the scientific, technical, and practical issues associated with running relative BFE simulations in AMBER20, with a focus on real-world drug discovery applications.

Figure 1:

Performance of AMBER20 for standard MD and thermodynamic integration (TI) on GeForce 1080Ti and 2080Ti graphics cards compiled on CUDA 9.1 using a Monte Carlo barostat, Langevin thermostat, and 4 fs time step.

Figure 2:

Hydration free energies of molecules containing alcohol functional groups. The bespoke force field uses the same parameters for bond stretch, bond angle, and van der Waals interactions as GAFF2. It derives the partial charges by fitting to the electrostatic potential computed using restricted Hartree-Fock with the 6–31G* basis set. The torsional parameters are then optimized to fit the potential energy surface computed by B3LYP/6–31G** for conformations generated at different torsional values of the rotatable bonds. Refitting the torsional parameters improves the agreement between the predicted and experimental hydration free energies.

Figure 3:

Hydration free energies of molecules containing aromatic halogens and aromatic nitrogens. Including virtual sites on the halogen and nitrogen atoms improves fitting to the electrostatic potential calculated by quantum chemistry and the agreement between the predicted and experimental hydration free energies.

Figure 4:

An example of ligand poses (purple carbons) docked A) not using core-restraints and B) using core-restraints. Employing core-constraints ensures that the binding mode is conserved between all of ligands in a congeneric series.

Figure 5:

Potential mappings using A) 2D or B) 3D information related to differing ortho substitutions on a terminal phenyl ring.

Figure 6:

Representative TI integrands for different λ schedules used in binding free energy calculations. The shape of the curve is highly dependant on the use of a single step versus multistep protocol. For simple charge changing transformations the curves may have a near linear character (fit dashed lines). The number of atoms being transformed (e.g. RBFE versus ABFE) also has a strong effect (left and right columns, respectively). The specific transformations are the Tyk2 ejm-47 (ABFE) and p38 2v→3fhm (RBFE) perturbations from the Wang, et al. data set. The specific coupling protocols are as outlined in Section 2.3. Note that decharge and recharge use opposite conventions for direction (λ=1 is fully coupled for decharge and λ = 0 is fully coupled for recharge).

Figure 7:

Comparison of 〈dU/dλ〉 curves from (leftmost column:) the original SSC(2) scheme with (m=2, n=6), from (middle column:) the modified SSC(2) scheme with (m=2, n=2), and from (rightmost column:) the modified SSC(2) scheme with (m=1, n=1) (defined in Eq.(20)). The β value is 12 Ų and α^Coul is 1. The molecular systems are upper row: the absolute hydration free energy of diphenyltoluene; middle row: the relative hydration free energy between the Factor Xa ligand L51h and L51c lower row: the absolute hydration free energy of a single Na⁺ ion.

Figure 8:

Comparison of 〈dU/dλ〉 curves from (leftmost column:) the original SSC(2) scheme with (m=2, n=6), from (middle column:) the modified SSC(2) scheme with (m=2, n=2), and from (rightmost column:) the modified SSC(2) scheme with (m=1, n=1) (defined in Eq.(20)). The β value is 50 Ų and α^Coul is 4. The molecular systems are upper row: the absolute hydration free energy of diphenyltoluene; middle row: the relative hydration free energy between the Factor Xa ligand L51h and L51c lower row: the absolute hydration free energy of a single Na⁺ ion.

Figure 9:

Seven types of virtual sites will be made available in a future release of AMBER. (a) aromatic halogens: fixed-distance VS from 2-atom frame; (b) aromatic halogens: flexible-distance VS from 2-atom frame; (c) aromatic nitrogens: flexible-distance VS from 3-atom frame; (d) aromatic nitrogens: fixed-distance VS from 3-atom frame; (e) aromatic nitrogens: fixed-distance-with-angle VS from 3-atom frame; (f) aromatic carbons: out-of-plane VS from 3-atom frame; (g) amines: in(out-of)-pyramid VS from 4-atom-frame. The virtual sites are shown as cyan beads; P,fn (n=1,2,3) are parent and frame atoms to define virtual sites. The relative positions of the virtual sites are specified by the illustrated geometric parameters.

Figure 10:

RBFE results for CDK2 (16 ligands). The left pane computes each edge of the RBFE network from independent MBAR optimizations. The center pane simultaneously optimizes all edges in the network, coupling the results through 22 cycle closure constraints. The right pane further includes a constraint that forces the RBFE of “ligand 28” to match experiment.