Enhanced ligand sampling for relative protein-ligand binding free energy calculations

Abstract

Free energy calculations are used to study how strongly potential drug molecules interact with their target receptors. The accuracy of these calculations depends on the accuracy of the molecular dynamics (MD) force field as well as proper sampling of the major conformations of each molecule. However, proper sampling of ligand conformations can be difficult when there are large barriers separating the major ligand conformations. An example of this is for ligands with an asymmetrically substituted phenyl ring, where the presence of protein loops hinders the proper sampling of the different ring conformations. These ring conformations become more difficult to sample when the size of the functional groups attached to the ring increases. The Adaptive Integration Method (AIM) has been developed, which adaptively changes the alchemical coupling parameter λ during the MD simulation so that conformations sampled at one λ can aid sampling at the other λ values. The Accelerated Adaptive Integration Method (AcclAIM) builds on AIM by lowering potential barriers for specific degrees of freedom at intermediate λ values. However, these methods may not work when there are very large barriers separating the major ligand conformations. In this work, we describe a modification to AIM that improves sampling of the different ring conformations, even when there is a very large barrier between them. This method combines AIM with conformational Monte Carlo sampling, giving improved convergence of ring populations and the resulting free energy. This method, called AIM/MC, is applied to study the relative binding free energy for a pair of ligands that bind to thrombin and a different pair of ligands that bind to aspartyl protease β-APP cleaving enzyme 1 (BACE1). These protein-ligand binding free energy calculations illustrate the improvements in conformational sampling and the convergence of the free energy compared to both AIM and AcclAIM.

Figure 1

Two ligands bound to BACE1 (gray ribbon) in (a) the initial state at λ = 0 and (b) the final state at λ = 1. Ligand 17a has a 2,5-dichlorophenyl group, while ligand 24 has a propynyl group attached to a pyridine ring. In the initial state, ligand 17a fully interacts with the rest of the system (opaque spheres), while the propynyl group of ligand 24 is decoupled (transparent spheres) and does not interact with the rest of the system. In the final state, the chlorine atoms of ligand 17a are decoupled (transparent spheres), while the propynyl group fully interacts (opaque spheres). The acceptance probability after a rotation when the molecule is in the initial state will be high, because of the symmetry of 2,5-dichlorophenyl. In this state, the propynyl atoms are decoupled and do not have charge or vdW interactions with the protein or water molecules, preventing steric clashes. The acceptance probability after rotation in the final state will be much lower, due to a higher chance of steric clashes between the propynyl group and other molecules in the system. Tachyon in visual molecular dynamics was used for rendering.

Figure 2

Structures of the systems studied in this work. (a) Thrombin (gray ribbon) bound to the ligands CDA and CDB. The ligands differ in that CDB has a methyl group, represented here in orange, attached to the P1 pyridine ring. This ring has two possible conformations, one where the methyl group points In toward the protein as shown here and the other where the ring has flipped out. (b) BACE1 (gray ribbon) bound to the ligands 17a and 24. Ligand 17a has a 2,5-dichlorophenyl group, while ligand 24 has a 5-(Prop-1-yn-1-yl)pyridin-3-yl group. The atoms unique to ligand 24 are shown as transparent orange. This ring has two possible conformations, either In as shown or Out where it is flipped 180°. Tachyon in visual molecular dynamics was used for rendering.

Figure 3

Fraction of the thrombin ligand P1 pyridine ring in the In conformation as a function of simulation time at (a) λ = 0, corresponding to CDA and (b) λ = 1, corresponding to CDB. The purple series shows the results for the simulations starting with the In conformation, and the red series shows the results for the simulations starting with the Out conformation. The AIM and AcclAIM results have been reported previously and are reproduced here for comparison. Error bars represent the standard deviation over four independent trials.

Figure 4

Relative binding free energy for the transformation of thrombin ligand CDA to CDB as a function of simulation time. The purple series shows the results for the simulations starting with the In conformation, and the red series shows the results starting with the Out conformation. The AIM and AcclAIM results have been reported previously and are reproduced here for comparison. Error bars represent the propagation of the error according to eq 5.

Figure 5

Fraction of the BACE1 ligand ring in the In conformation as a function of simulation time at (a) λ = 0, corresponding to ligand 17a and (b) λ = 1, corresponding to ligand 24. The purple series shows the results for the simulations starting with the In conformation, and the red series shows the results starting with the Out conformation. Error bars represent the standard deviation over four independent trials.

Figure 6

Relative binding free energy for the transformation of BACE1 ligand 17a to 24 as a function of simulation time. The purple series shows the results for the simulations starting with the In conformation, and the red series shows the results starting with the Out conformation. Error bars represent the propagation of the error according to eq 5.