Funnel metadynamics as accurate binding free-energy method

Abstract

A detailed description of the events ruling ligand/protein interaction and an accurate estimation of the drug affinity to its target is of great help in speeding drug discovery strategies. We have developed a metadynamics-based approach, named funnel metadynamics, that allows the ligand to enhance the sampling of the target binding sites and its solvated states. This method leads to an efficient characterization of the binding free-energy surface and an accurate calculation of the absolute protein-ligand binding free energy. We illustrate our protocol in two systems, benzamidine/trypsin and SC-558/cyclooxygenase 2. In both cases, the X-ray conformation has been found as the lowest free-energy pose, and the computed protein-ligand binding free energy in good agreement with experiments. Furthermore, funnel metadynamics unveils important information about the binding process, such as the presence of alternative binding modes and the role of waters. The results achieved at an affordable computational cost make funnel metadynamics a valuable method for drug discovery and for dealing with a variety of problems in chemistry, physics, and material science.

Fig. 1.

(A) Schematic representation of the ligand/protein binding process and the funnel restraint potential used in FM calculations. The shape of the funnel can be customized on the target by setting a few parameters. In particular, given z, the axis defining the exit-binding path of the ligand, z_cc is the distance where the restraint potential switches from a cone shape into a cylinder. The α-angle defines the amplitude of the cone, and R_cyl is the radius of the cylindrical section. (B) The funnel restraint potential applied to trypsin (Upper Left) and COX-2 (Upper Right) enzymes with the ligands considered in the study, benzamidine (Lower Left) and SC-558 (Lower Right). In the trypsin case, α is 0.55 rad and z_cc is 18 Å (Table S1). In the COX case, α is 0.6 rad and z_cc is 44 Å (Table S2). In both cases, R_cyl is set to 1 Å.

Fig. 2.

The FES of the benzamidine/trypsin binding is computed using a reweighting algorithm (50) as a function of the projection on the z axis of the ligand center of mass and a torsion CV, where z is the axis of the funnel (Table S1). Isosurfaces are shown every 1 kcal/mol. (Insets A–C) Conformations representing the bound poses, basin A (Inset A) and basin B (Inset B). C shows one of the isoenergetic conformations representing the unbound state. In the unbound state, z > 20 Å, the ligand can assume a large number of orientations, which are represented by states with similar energy values, as shown in the FES.

Fig. 3.

The FES of the SC-558/COX-2 binding is computed using a reweighting algorithm (50) as a function of the projection on the z axis and distance from z axis of the center of mass of the ligand, where z is the axis of the funnel (Table S2). Isosurfaces are shown every 3 kcal/mol. (Insets A and B) Binding modes corresponding to the two deepest energy basins, basin A (X-ray) and basin B (alternative pose). (Upper) Ligand position relative to the enzyme during FM simulations. The spheres represent the ligand center of mass and are colored according to their corresponding free-energy values. For clarity, only selected frames are shown.