Funnel metadynamics as accurate binding free-energy method

Abstract

Studying the molecular interactions between a drug and its target helps in understanding the target functional mechanism and offers the possibility for exogenous control of its physiological activity. In recent years, a vast experimental and computational effort has revealed in ever-more-precise detail the ligand/target recognition mechanism (1, 2). In this context, an accurate estimation of the ligand-binding affinity is in great demand because it would facilitate many steps of the drug discovery pipeline, such as structure-based drug design and lead optimization; this is not, however, a simple task. In fact, an accurate estimation of the binding affinity or, equivalently, the absolute protein–ligand binding free energy, requires an accurate description of the ligand/protein interactions, their flexibility, and the solvation process. Many methods have been proposed to tackle this problem. For instance, docking protocols are widely used to generate and rank candidate poses based on empirical scoring functions, either physically or statistically based (3–5). These techniques have been proven to be highly efficient in screening a large number of compounds in a short time (6); this, however, at the price of limited accuracy in estimating affinities (7).

Alternatively, a variety of methods to describe ligand/protein interactions in a more accurate way at higher computational cost have been proposed. These techniques can be grouped in two categories: (i) endpoint and (ii) pathway methods. The former group is composed of those techniques that sample ligand and protein in unbound and bound states and compute the protein–ligand binding free energy by taking the difference between the absolute free energy of these two states. Examples include microscopic linear response approximation (8), linear interaction energy (9, 10), protein dipoles Langevin dipoles (11), as well as molecular mechanics Poisson–Boltzmann surface area, and generalized Born surface area (12).

At variance with endpoint methods, in pathway methods, the ligand is gradually separated from the protein. The binding free energy is then obtained by summing different contributions coming from a discretized path that connects the initial and final state. This class includes methods in which the ligand/protein interactions are gradually switched off, such as thermodynamic integration (13), free-energy perturbation (14, 15), double-decoupling method (16), and double-annihilation method (17). Techniques such as steered molecular dynamics (SMD) (18) and umbrella sampling (19), where the ligand and the protein are physically separated from each other, also belong to this group. While in SMD, the ligand is dragged out from the protein using a moving restraining potential, in umbrella sampling, the path from the bound to the unbound state is divided in a finite number of windows, which are independently sampled.

Though these methods have been successfully used to compute the ligand binding free energy in many cases (20–22), the requirement of knowing in advance the binding mode hampers a more general applicability. The intensity of the efforts in developing these methods reflects both the great potential of these calculations and their difficulties. In particular, the difficulties arise mainly from the fact that the ligand/protein binding process is a rare event, difficult to sample with standard techniques such as molecular dynamics (MD). Even the most ambitious efforts in this direction, though revealing precious details of the binding process (23, 24), have not been able to determine accurately the binding energy. To achieve this result, the use of enhanced sampling methods is mandatory.

Among the emerging techniques, metadynamics (25) has proven to be very useful in studying long-timescale processes (26, 27), particularly in complex ligand/protein binding cases (28–30). Metadynamics works by adding an external history-dependent potential that acts on few degrees of freedom, named collective variables (CVs). In such a way, the sampling is accelerated, and the free-energy surface (FES) of the process can be calculated from the added potential. Unfortunately, only a qualitative estimation of the protein–ligand binding free energy could be obtained for the binding processes studied so far (28, 31). In fact, once the ligand leaves the binding pocket, it has difficulty finding its way back, and starts exploring all of the possible solvated states. These conformations represent a vast part of the configuration space that cannot be sampled thoroughly in a limited computation time. Therefore, once out, the ligand does not again find the binding site, and multiple binding/unbinding events, which are the key to an accurate determination of the binding free energy in metadynamics, cannot be observed.

Here, we present a metadynamics-based approach, named funnel metadynamics (FM), which overcomes all these limitations and allows an accurate estimation of the absolute protein–ligand binding free energy. In particular, in FM, a funnel-shaped restraint potential is applied to the system, reducing the space to explore in the unbound state. In such a way, the sampling of ligand-bound and -unbound states is highly enhanced, thus leading to an accurate estimation of the binding FES within a reasonable simulation time. We have used FM to study the binding process in two different systems: benzamidine/trypsin and SC-558/cyclooxygenase 2 (COX-2). The latter system represents a particularly challenging case study with the coexistence of protein motions, solvent effect, and complex binding pathways. In both cases, the X-ray pose turned out to have the lowest free-energy value, and the computed absolute protein–ligand binding free energy was in good agreement with experiments. Furthermore, from FM simulations, important information on the binding process has been disclosed, such as the role of waters in the benzamidine/trypsin case and the presence of an alternative ligand binding mode in COX-2.

Using FM, no a-priori information about the ligand binding mode is required, and the exploration capability of the original method remains unaltered. Thus, the exploration of buried binding sites is possible, while taking into account slow protein motion and solvent effects. The present protocol represents a most valuable method to study ligand/protein interaction, and its relatively low computational cost renders its use appealing even in industrial applications where speed is valued.

Theoretical Background

Funnel Idea.

Sampling ligand/protein binding with all-atoms MD simulations in explicit solvent is extremely attractive because it can provide molecular information at high resolution. Unfortunately, these processes usually have a long timescale and therefore cannot be sampled at reasonable computational costs. Nevertheless, the advance in computer power and the advent of graphics processor units have allowed the study of ligand/protein binding processes using plain MD simulations, however still at great computational cost and using dedicated hardware (23, 24). In fact, tens of microseconds of MD simulations were necessary to collect enough statistics to describe the ligand binding process.

An alternative method is to use enhanced sampling techniques to access long-timescale events within a reasonable computational time. One of these techniques is umbrella sampling, which has largely been used to study ligand/protein interactions and to compute the absolute protein–ligand binding free energy (20, 21). However, this technique fails in exploring thoroughly the fully solvated state of the ligand. This limitation has been remedied by using a cylindrical restraint potential to reduce the sampling space (32, 33). The effect of the external potential can be rigorously taken into account (34), and the binding constant K_b in the presence of the restraint is given by

Here, ΔG_site is the change in the free energy of the bound state caused by the presence of the restraint, β = (k_B T)⁻¹, k_B is the Boltzmann constant, T the temperature of the system, and S_u is the cross-section of the cylinder. The potential W(z) and its value in the unbound state, W_ref, can be derived from the potential of mean force (PMF). Eq. 1 provides an unbiased estimator of K_b, independent of the choice of the restraint potential. If the radius of the cylinder is chosen to be much larger than the lateral fluctuations of the ligand in the binding site, the restraint potential is not felt by the ligand in the binding site, thus ΔG_site = 0. Unfortunately, the choice of the radius is not simple. In fact, though a small radius value limits the exploration in the unbound state, it reduces also the exploration of the binding site. However, choosing a large radius is not advantageous because it increases the space to be sampled in the unbound state and hence the computational time. As a result, the PMF calculation might be affected by this choice and, consequently, the estimation of the protein–ligand binding free energy (21).

To overcome this limitation, we have developed a funnel-shaped restraint potential that can be applied to the target protein. This potential is a combination of a cone restraint, which includes the binding site, and a cylindric part, which is directed toward the solvent (Fig. 1). Using the funnel potential during the simulation, the system does not feel any repulsive bias when the ligand explores regions inside the funnel area. As the ligand reaches the edge of the funnel, a repulsive bias is applied to the system, disfavoring it from visiting regions outside the funnel. As can be seen in Fig. 1, if the shape of the funnel is properly chosen, the sampling at the binding site is not affected by the external bias, whereas in the bulk water the repulsive potential reduces the space to be explored to a cylindric region; this favors the observation of multiple binding/unbinding events leading to a faster convergence of the results.

Open in a separate window

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 Å.

Furthermore, if the funnel restraint potential is used in free-energy calculations with techniques such as metadynamics, the free-energy difference between bound and unbound states depends exclusively on the free-energy value at the two states, independently from the path that connects one state to the other; this is advantageous with respect to other methods, such as umbrella sampling. In fact, in umbrella sampling, where the PMF is calculated using the different contributions coming from the simulation windows that connect the bound state to the unbound one, all of the simulation windows in the bound states must not feel the restraint potential to provide an easy estimate of the PMF (ΔG_site = 0). If this condition is not fulfilled, the PMF calculation and, consequently, the binding free-energy estimate, is more complex (ΔG_site ≠ 0) and depends on the chosen path.

Absolute Protein–Ligand Binding Free Energy.

In free-energy calculations, the absolute protein–ligand binding free energy ΔG_b⁰ is typically computed using the following formula:

where K_b is the equilibrium binding constant, and C⁰ = 1/1,660 Å⁻³ is the standard concentration. Using metadynamics calculations with the funnel restraint potential and using Eq. 1 with some rearrangements, as reported in Methods, Eq. 2 becomes

where, ΔG is the free-energy difference between bound and unbound states, and πR²_cyl is the surface of the cylinder used as restraint potential. With β and C⁰ being constant, the absolute protein–ligand binding free energy is equal to ΔG minus the analytical correction in Eq. 3 (see Methods for details). It is important to stress that the funnel restraint potential ensures a number of recrossing events between the different states visited by the system during the simulation, leading to a quantitatively well-characterized free-energy profile and a converged estimation of ΔG.

Results

We have used the funnel restraint potential in combination with well-tempered metadynamics (35), hereafter named FM, to study two ligand/protein binding cases and compute the absolute protein–ligand binding free energies. The first case is the trypsin/benzamidine complex, mainly used as a reference model, and the second is COX-2 in complex with the potent inhibitor SC-558. The latter, which has been previously studied by us (28), represents a challenging ligand/protein binding case.

Bound State.

The whole sampling took ∼0.5 µs of FM simulations. Looking at the FES computed as a function of the projection on z axis, where z is the axis of the funnel restraint potential and a torsion CV (Methods; Table S1), two main basins can be detected (Fig. 2). The deepest one, basin A in Fig. 2, corresponds to the benzamidine in its crystal pose, where a number of strong interactions with the enzyme are formed (Fig. 2; SI Text). It is interesting to note that in this pose, a water molecule is present in the binding pocket and forms an H-bond network with the side chains of Tyr228, Asp189, and Ser190. A water molecule is present in a similar position in many X-ray structures (e.g., PDB ID code 3atl) (38), thus suggesting a structural role for this molecule in the trypsin binding site (Fig. S1). The second energy minimum, basin B, is ∼1 kcal/mol (1 kcal = 4.18 kJ) higher than basin A. Here, the ligand is slightly rotated in the binding site. Overall, however, the interactions established in pose A are conserved (Fig. 2). In particular, the diamino group of benzamidine engages a direct H-bond with Ser190, and a water bridge interaction is established with the Asp189 side chain. Furthermore, the aromatic ring of the ligand is involved in an interaction with the sulfur atoms forming the disulfide bridge between Cys191 and Cys220. It is interesting to note that in basin B the diamino group of benzamidine engages H-bond interactions with the carbonyl oxygens of Val227, Val213, and Ser214 via two water molecules. Water molecules located at similar positions in the binding site can be found in the X-ray structures of trypsin in the apo and ligated form (PDB ID codes 1s0q and 3atl, respectively; Fig. S1).

Open in a separate window

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.

Our results show an important role played by waters during ligand binding, and a similar functional role has been reported also in other studies (29, 39). Therefore, the use of atomistic simulations with explicit solvent is mandatory to take into account the solvation effect and have an accurate estimation of ligand/protein interactions.

The stability of the basin B pose has been further assessed through an over 100-ns-long unbiased MD simulation. During the whole simulation, the ligand binding mode is stable conserving all of the interactions described above (SI Text; Fig. S2). The depth of basin B and the good stability of the ligand/protein interactions lead us to consider this pose the first binding event of benzamidine in the active site before reaching its final position in basin A. Alternatively, the basin B pose can be considered the first unbinding event of the ligand from the catalytic site of the enzyme. It is important to stress that both basin A and B are within the cone region, and their exploration is not affected by the funnel restraint potential (Fig. S3).

Finally, we note that in a recent study by Söderhjelm et al. (37), a binding mode highly similar to pose B was found, and that in our simulations the ligand often occupies, along its way to the binding site, positions close to the states described in this study.

A movie showing the benzamidine binding/unbinding to trypsin under the action of FM is provided (Movie S1).

SC-558/COX-2 System.

Encouraged by the results obtained on the benzamidine/trypsin system, we decided to use FM in a more complex case study: the binding of SC-558 to COX-2. This study is particularly challenging because the binding site is buried, and substantial motion of the protein scaffold is needed to allow the ligand in. Furthermore, the ligand needs to dispose of its solvation water to access to the binding site. This system has already been studied by us, and the existence of a second binding mode different from the crystallographic one has been revealed (28). In this work, we could not give an accurate estimation of the binding energy due to the difficulty of exploring the solvated state. However, the binding cavity and the free-energy difference between the two poses were properly and fully described.

We used a CV setting similar to that of Limongelli et al. (28). In particular, the projection on the z axis of the funnel of the ligand center of mass and a torsion have been chosen as CVs to describe the ligand motion with respect to the protein. Furthermore, a path CV has been used to take into account the motion of the helices that form the gate to access the active site (Methods; Tables S2 and S3). As in the trypsin/benzamidine case, the funnel parameters (z_cc and α) were chosen so as not to alter the potential seen by the ligand inside the binding site (Fig. 1; Methods).

The FES was converged after ∼350 ns, and for the part inside, the binding cavity was very similar to that described in ref. 28. In fact, as in ref. 28, two binding poses were found, basin A, which corresponds to the X-ray pose, and basin B, that was first described in ref. 28 and corresponds to an alternative binding mode of SC-558 in COX-2 (Fig. 3; SI Text). As can be seen in Fig. 3, these two basins are within the cone region, and their exploration is not affected by the funnel restraint potential. Contrary to ref. 28, we can now thoroughly sample the unbound state, which is reached at z > 44 Å from the binding site. Here, the ligand is outside the interaction range of the protein and can assume a large number of orientations. As seen in the trypsin case, in the unbound state the FES is flat, corresponding to different states of the system with similar free-energy values.

Open in a separate window

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.

The free-energy difference between bound and unbound states, ΔG, of SC-558 is equal to −14.9 kcal/mol, which, considering the analytical correction of 3.8 kcal/mol, leads to an estimate of the absolute binding free energy of −11.1 ± 1.5 kcal/mol (Fig. S5). This value is in good agreement with the experimental one of ∼−12 kcal/mol, derived from the ligand IC₅₀ value (42) via the Cheng–Prusoff relation ΔG_b⁰ = −K_b T ln IC₅₀/(1 + (S/K_m)) (43) using the substrate concentration (S) and the concentration of substrate at which enzyme activity is at half maximal (K_m), reported in ref. 44.

These results show that FM preserves the exploration capability of metadynamics in the binding cavity. At the same time, sampling of bound and unbound states is considerably enhanced (Fig. S6), and an accurate estimation of the absolute protein–ligand binding free energy can finally be computed. One might note that COXs are monotopic membrane proteins, and that in our model we do not consider the membrane. This represents an approximation that can be however tolerated, because COXs binding site is hydrated with the waters playing a key role in the cyclooxygenase site (45). Furthermore, most COX inhibitors (e.g., aspirin, ibuprofen, flurbiprofen) are not lipophilic enough to cross the membrane, suggesting that they can follow an extracellular binding pathway to reach the catalytic site. Therefore, this approximation does not affect the estimate of the ligand/protein binding free energy, and the agreement with the experimental value supports the reliability of our model. In this difficult case, in which the binding site is buried and the use of a larger number of CVs is necessary, FM has proven to be very powerful in studying the binding process and accurately computing the protein–ligand binding free energy.

Discussion

Here, we have developed a metadynamics-based approach, named FM, which allows us to compute the absolute protein–ligand binding free energy, leaving unaltered the exploration capability of the method. Thanks to a funnel-shaped restraint potential, the sampling of bound and unbound states is highly enhanced, leading to a quantitatively well-characterized binding free-energy surface within a reasonable simulation time. We have used FM to sample the ligand binding process in two different systems: benzamidine/trypsin and SC-558/COX-2. In both cases, the X-ray pose was the lowest-energy pose, and the estimated binding free energy was in good agreement with experiments. In both systems, the whole simulation of the binding process took less than 0.5 µs, which is definitely less demanding than other MD-based approaches that take tens of microseconds of simulations (23, 24). Furthermore, FM allows us to compute the absolute protein–ligand binding free energy by overcoming some of the difficulties related to other methods. For instance, with respect to umbrella sampling, FM allows the use of a larger number of collective variables, as in the COX case, still achieving convergence of the results at an affordable computational cost; this allows us to describe the slow degrees of freedom of the system, which is often necessary in complex ligand/protein interactions. Furthermore, using FM, information about the ligand binding mode is not required in advance, although the approximate location of the binding pocket in the target structure should be known. However, if this information is not available, FM can be combined with other methods, such as virtual screening protocols, or more advanced techniques (37), that are more efficient in finding candidate binding sites. Furthermore, FM can also be combined with methods, such as free-energy perturbation, which are able to estimate the ligand binding affinity for a series of analogs (46); and whereas the former helps in finding the correct ligand binding mode with an accurate energy estimation of the binding event, the latter can be used to assess changes in the ligand/protein binding affinity with respect to different ligand substitutions.

We note that the application of FM can go beyond ligand/protein binding studies. In fact, it can be exploited in many other research fields where two-body interactions are important. For instance, FM simulations can be used to study the interactions between atoms in crystal growth (47). Furthermore, the funnel-shaped restraint potential can be used in combination with other enhanced sampling methods, such as bias exchange, parallel tempering, or replica exchange (26), to achieve convergence in even shorter computational time. Finally, one might customize on the studied system the shape of the restraint potential, changing the cone part and leaving unaltered the cylindrical region of the potential for the unbound state.

The results achieved here at a reasonable computational cost make FM a valuable method to tackle a wide variety of problems in chemistry, physics, and material science.

Methods

FM Theory.

The absolute protein–ligand binding free energy ΔG_b⁰ can be defined in terms of the equilibrium binding constant K_b, as seen in Eq. 2. Following the approach of Allen et al. (32), we can write K_b in terms of the PMF W(r) as

where r' is a reference ligand position in the bulk water, and H_site(r) is equal to 1 inside the binding site and 0 elsewhere. We want to express K_b in terms of the 1D PMF W(z), which is obtained by integrating over the allowed lateral displacement in the xy space. We assume that the binding site is contained in the range z_min ≤ z ≤ z_max. The allowed displacement is determined by a hybrid truncated conical and cylindrical restraint H_R(x,y,z) with axis aligned to the z axis, as follows:

with the condition z_cc > z_max. The truncated conical restraint is defined as

with R(z) = R_cyl + tan (α) (z_cc − z), and α modulating the aperture of the cone. The cylindrical restraint is defined as

The 1D PMF W(z) is given by

We define the 1D PMF to be zero in the bulk water at z = z_bulk, with z_bulk > z_cc. In this region, H_R(x,y,z) = H_cyl(x,y), and the constant C in Eq. 8 can be written as

Because in the bulk water, W(x,y,z) is independent from x and y, we can write

and the constant C can be written as

By using the expression above, we can then rewrite Eq. 8 as

Let us assume that in the binding site, the (truncated conical) restraint does not act on the ligand, i.e., the free-energy cost introduced by the restraint is zero when the ligand is in the binding site. We can write

where

If we substitute Eqs. 13 and 12 in Eq. 4, we obtain an expression for the binding constant in terms of the PMF W(z):

FM simulations have been carried out with NAMD code (48) using the PLUMED plugin (49). Details on FM simulations are reported in SI Text.

Supplementary Material

Acknowledgments

Footnotes

References