Metadynamics-Based Approaches for Modeling the Hypoxia-Inducible Factor 2α Ligand Binding Process

Molecular Docking of the KG-721 Ligand

Conformational analysis of the ligand structure was performed using Macromodel⁷⁴ with the OPLS_2005⁷⁵ force field. The obtained global minimum was used as the starting point for molecular docking calculations using Glide⁷⁶ XP⁷⁷ (Extra Precision). In particular, Glide uses a flexible ligand-rigid protein approach, in which a series of hierarchical filters are applied to find the possible positions and conformations of the ligand in the binding cavity (poses). The properties of the protein are represented on a grid that provides gradually more accurate scores. The initial screenings are deterministically performed over the complete phase space of the ligand to identify the most promising poses. From the selected poses, the ligand is then refined in the torsional space in the receptor field. To take into account the flexibility of the protein, the ensemble-docking approach was used, which involves ligand docking to multiple receptor conformations. These can be derived either experimentally or computationally (e.g., by MD simulations).⁵ The conformational ensemble here selected consisted of the crystallographic structures of the HIF-2α PAS-B in complex with artificial ligands available in PDB (3F1O,⁵⁰3H82,⁵³3H7W,⁵³4GS9,⁵¹ and 4GHI(52)). The results showed that the best XP score is the one related to the KG-721 ligand in the 4GHI structure.

Steered MD Simulations (SMD)

In SMD simulations,^26,27,78 a time-dependent external force is applied to the ligand to aid its unbinding from the protein. In particular, the transition between the bound and the unbound states is achieved by adding a harmonic time-dependent potential, acting on a descriptor (or collective variable), to the standard Hamiltonian. During the transition, the external work performed on the system can be calculated using the Jarzynski equation.⁷⁹ All the SMD simulations were performed using the PLUMED 2.4.6^80,81 plugin integrated in GROMACS 2018.6.⁶⁶ We chose the ligand-protein distance as the pulling variable. This was defined as the distance between the center of mass of selected atoms at the bottom of the binding cavity (different for the two pathways, see Figure S2) and the center of mass of the ligand heavy atoms. The spring constant was set to the value of 10.0 kcal/mol·Å², and the ligand was pulled from the initial value of CV to 35 Å in 25 ns with a resulting pulling velocity of 0.984 Å/ns. We ran 50 independent replicas, and the time length for each simulation was 25 ns, which ensured the achievement of a complete solvation of the ligand in the unbound state. The starting point of each replica was derived from an ensemble of states extrapolated at regular time intervals of 0.2 ns from the last 10 ns of the unbiased simulation.

Metadynamics (MetaD) and Path Collective Variables (PCVs)

The central idea of the metadynamics method^19,30 is to bias the system along a set of CVs using a history-dependent potential. To achieve this, a Gaussian-shaped potential is added to bias the system at the current position of the CVs, at regular time intervals. This allows the system to escape from any local minimum and to visit new regions in the CVs space. In metadynamics, to push the system to visit even high free-energy regions, the Gaussian-shaped potential has constant height. On the contrary, in the well-tempered metadynamics⁸² approach, used in this work, the height of the Gaussian is decreased with the amount of bias already deposited according to

where w₀ is an initial Gaussian height, ΔT an input parameter with the dimension of a temperature, k_B is the Boltzmann constant, and τ_G is the time interval at which Gaussians are deposited.⁸²

The path CV formalism^32,83 has been widely used to investigate biological processes, to compute their free-energy surfaces, and to characterize their kinetic behavior.^39,34 In this work, PCVs were used to study the transition between the bound and the unbound states in the unbinding process of some HIF2-α ligands. We described the transition pathway with a set of frames derived from the SMD simulations: 12 frames were used for the THS-020 ligand and 11 frames for the KG-721 ligand. For the first part of the path, the frames (from 1 to 7 for THS-020 and from 1 to 6 for KG-721) were obtained from the SMD simulations with the lowest value of external work performed on the system. Frames were selected to be equally spaced (2 Å). For the second part of the path, the frames were obtained with linear interpolation (see the Supporting Information, Supplementary Text for the details). Following the procedure proposed by Branduardi et al.,³² we introduced two collective variables: s(R), the progress along the reference path, and z(R), the distance orthogonal to the reference path. The λ value was set to 33.0 nm². The distance between the instantaneous conformational state during the simulation and the reference coordinates in the path was evaluated by the RMSD metric.⁸⁴ In particular, the RMSD along the entry/exit pathway was calculated between a selection of protein atoms and all the ligand heavy atoms (see Figure S3). In all simulations, the Gaussian-shaped potentials were deposited every 500 simulation steps, the initial height was set to 1 kJ/mol, and the decay corresponding to a bias factor of 10 was chosen. The Gaussian widths (σ) for the s(R) and z(R) variables were set to 0.05 and 0.007, respectively. Widths were set so that they are about 1/3 of the CV standard deviations observed in the unbiased MD simulation. The two variables, s(R) and z(R), were constrained to be less than 12 and 0.2 nm², respectively.

Identification of Unbinding Pathways for the THS-020 Ligand by the SMD Method

Starting from the X-ray structure for the HIF-2α PAS-B domain in complex with THS-020, the possible ligand unbinding pathways were investigated using the SMD approach. Steered molecular dynamics is a popular method for studying ligand-protein unbinding events⁸⁹⁻⁹¹ and can provide both a qualitative description of the pathways and a quantitative estimate of the free-energy difference between the bound and unbound states.

To calculate the free-energy difference using the Jarzynski equality, it is necessary to have a large number of SMD replicas. To this aim, a 20 ns unbiased MD simulation was performed starting from the X-ray structure, generating an ensemble of 50 slightly different states of the complex (Figure S4a), extracted from the last 10 ns. These were then used as starting points for the 50 SMD simulations. The structural convergence of the unbiased simulation was assessed by calculating the RMSD matrix on the protein Cα atoms and on the ligand heavy atoms (Figure S4b). Each SMD replica was 25 ns long. To verify if the values chosen for the SMD simulation parameters (see the Methods section) are appropriate, the RMSD plot on Cα atoms (Figure S5) and the secondary structure conservation graphs (Figure S6) were calculated for the 50 replicas along path 1. As shown in Figure S6, no significant distortions of the protein structure (except for a slight deformation of the Fα helix upon ligand unbinding) were observed during simulations, thus confirming the validity of the proposed protocol.

Other authors identified two entry/exit pathways for solvent water by MD simulations of the apo HIF-2α PAS-B⁵³: path 1 gets through the Fα helix and the Gβ strand, while path 2 gets through the Fα helix, the short Eα helix, and the AB loop (Figure ​Figure22). On this basis, for each of these two pathways, we calculated a CV allowing us to pull the ligand out of the binding cavity, by selecting an appropriate set of residues at the bottom of the binding cavity (see the Methods section).

Open in a separate window

Figure 2

THS-020 unbinding pathways. In pathway 1 (left), the ligand passes through Fα and Gβ while in pathway 2 (right) through Fα, Eα, and the AB loop. The starting protein structure is represented as gray cartoons, the ligand conformations in the first and last frames of the trajectory as blue sticks, and the conformations of the ligand in the intermediate frames as transparent sticks.

The work profiles resulting from the 50 replicas, in Figure S7, show an increase in the work value during the initial part of the unbinding process followed by relatively settled work values, indicating the absence of interaction between the ligand and the protein.

All the curves show a similar profile, but a qualitative comparison of the total work reveals that less work is required for unbinding following pathway 1. Moreover, a broader range of values is observed for the replicas following pathway 2, indicating that higher barriers can occur in some of the replicas associated to this pathway. For a quantitative analysis, the minimum and maximum work values (W_min and W_max), the minimum value of the maximal force (F_max) among the replicas, and the free-energy difference between the unbound and the bound states (ΔF_unbind) were extracted for each pathway (Table 1). The results show that the W_min necessary to pull the ligand outside the cavity along path 1 is about 10 kcal/mol less than that required for path 2; a similar trend is observed for the values of W_max and ΔF_unbind. A difference of 75.59 pN between the F_max values in the two paths is observed, which confirms a clear preference of path 1 over path 2.

Therefore, steered MD allowed us to compare the two unbinding pathways of THS-020 and to select the preferred one by identifying a higher-energy barrier along pathway 2. However, SMD provided a value for the ΔF_unbind, estimated by means of the Jarzynski equality, that was about 4 times higher than the experimental value. It is indeed known that these nonequilibrium simulations generally undersample the relevant protein-ligand states across the unbinding pathway, leading to errors in the computed binding free energy.⁹² Moreover, replicas with less work done on the system have an enormous weight compared to all the other trajectories, which makes the method extremely sensitive to insufficient sampling.⁹³ For this reason, we then applied the PCV MetaD approach that was recently proposed as a valuable method for computing absolute binding free energies in ligand binding.^34,39,35

Metadynamics Simulations and Free-Energy Profiles with PCVs for THS-020

For a detailed mechanistic interpretation of the ligand binding/unbinding process, we used well-tempered metadynamics simulations with the PCV approach^83,35,94 (see the Methods section). This allowed us to characterize the relevant states along the preferred path obtained with SMD simulations (the one with lower values of total work obtained from the SMD simulation, path 1), as well as to estimate the binding free-energy value. The key points of this method are the choice of appropriate CVs and the construction of a set of equally spaced frames along the CVs in terms of RMSD between adjacent snapshots. This frameset represents a reference path for investigating the process. As CVs, we used the progress along the path, s(R), and the distance orthogonal to the reference path, z(R). We want to underline the importance of the path construction phase, especially in a case with a buried binding site like the one presented in this work. Here, we decided to include both ligand and protein atoms in the frameset that represents the path to better describe the protein conformational changes during the process (mouth opening through side-chain conformational changes and small backbone adjustment). Only protein atoms involved in the conformational changes, highlighted by SMD simulations, were included. Moreover, a hybrid approach that combines frames from SMD simulations and linear interpolation was used for the inner and outer parts of the path, respectively. Details about the construction of the reference path can be found in the Supporting Information, Supplementary Text. The reference path obtained with this approach is represented in Figure S8. The RMSD matrix of the frameset (Figure S8) is a symmetric matrix with a typical gull wing shape, indicating that the frames are correctly equally spaced.

We collected a total of 1.8 μs of metadynamics simulation in which we observed several binding/unbinding events, as shown in Figure ​Figure33A. The binding free energy (ΔF_bind), calculated as the free-energy difference between the deepest minima in the bound state and the flat plateau in the unbound state, turns out to be equal to −11.8 kcal/mol. The free-energy profile during the simulation, shown in Figure ​Figure33B, indicates that the simulation reaches a constant value of ΔF_bind after about 1200 ns. The convergence was also monitored by plotting the hill heights as a function of the simulation time (Figure S9).

Open in a separate window

Figure 3

(A) Plot of the CV1 (s(R)) against simulation time during THS-020 MetaD simulation: the lowest values of s(R) correspond to the bound state while the highest to the unbound ones; (B) one-dimensional projection of the binding free-energy values associated to path 1 during the metadynamics simulation.

The free-energy surface (FES) for the binding/unbinding process, as a function of CV1 (s(R)) and CV2 (z(R)) in 1.8 μs of metadynamics simulation, is shown in Figure ​Figure44 together with the relevant minima found along the pathway (labeled as A to G). The coordinates and the binding free-energy values of each minimum are reported in Table S1. A cluster analysis of the conformations belonging to each minimum hole was then performed (see the Methods section for details). In each minimum, the first cluster is the most populated (Figure S10), and its centroid was used as the representative structure for that minimum. Looking at the FES (Figure ​Figure44), three different regions can be identified following CV1: the bound state, with s(R) values between 1 and 3 (minima A–C), the intermediate states, with s(R) values between 3 and 7 (minima D–G), and the unbound state, with s(R) values from 7 onward. In the deepest minimum (A), the ligand is oriented in the same way as in the X-ray structure (ligand RMSD = 1.36 Å). This geometry is stabilized by a hydrogen bond between the NH group of the ligand and the H248 residue as well as by a transient hydrogen bond between the oxygen atom of furan and the S246 residue.

Open in a separate window

Figure 4

Free-energy surface obtained from the PCV approach for the binding/unbinding of the THS-020 ligand. The isolines are drawn using a 1.5 kcal/mol spacing. The 3D structures of the centroids of the main minima are reported with different colors: the protein is represented as cartoons and the ligand as sticks. The black lines indicate the corresponding minima in the FES.

Moving up to higher CV2 values, alternative binding geometries can be detected. In minimum B, the ligand is shifted toward the exit of the cavity, and breaking of the hydrogen bonds that stabilize minimum A causes a lower stability. Instead, in minimum C, the ligand is located on the mouth of the binding cavity and is even turned 180°, with the CF₃ group oriented toward the bottom of the cavity. Also, in this minimum, the amino group of the ligand forms a hydrogen bond with the S292 residue. This last conformation represents the second minimum in energy.

We acknowledge the CINECA as part of the agreement with the University of Milano-Bicocca and the CINECA award under the ISCRA initiative (Project IsB18 ADAMS) for the availability of high-performance computing resources.