HIV-1 Protease Dimerization Dynamics Reveals a Transient Druggable Binding Pocket at the Interface

The binding mechanism of HIV-1 protease monomers leading to the catalytically competent dimeric enzyme has been investigated by means of state-of-the-art atomistic simulations. The emerging picture allows a deeper understanding of experimental observations and reveals that water molecules trapped at the interface have an important role in slowing down the kinetics of the association process. Unexpectedly, a cryptic binding pocket is identified at the interface of the complex, corresponding to a partially bound dimer that lacks enzymatic function. The pocket has a transient nature with a lifetime longer than 1 μs, and it displays very favorable druggability features. Docking as well as MM-GBSA free-energy calculations further support the possibility to target the new binding site by means of inhibitors able to prevent the complete dimerization by capturing the inactive conformation. This discovery could open the way to the rational design of a new class of anti-HIV drugs.

1. Structural relaxation in the presence of soft restraints (1 kcal·mol -1 ·Å -2 ) on all the nonhydrogenous atoms of the protein and the ligand. In the second and third steps, the restraints were kept only on backbone and C atoms, respectively, and on the non-hydrogenous atoms of ligand. Finally, restraints were removed from the ligand and from a selection of residues having at least one atom within 4 Å from the ligand. In each step the structure of the solute from previous run was used as target for restraints, and up to 25.000 optimization steps were performed using the conjugategradients algorithm.
2. Next, annealing up to 340 K was performed in 2 ns, using the same setup as in the last step of the relaxation described at the previous point, and constant volume and temperature conditions (NVT ensemble). This was followed by quenching to 310K in 1 ns, and then by a 1-ns long equilibration with same setup as above, but in the NTP ensemble.
3. Finally, a partially restrained MD run of 5 ns, using the same setup as above was performed. Trajectories were saved every 50 ps. A time step of 2 fs was used, and periodic boundary conditions were employed, and electrostatic interactions were treated using the particle-mesh-Ewald (PME) method, with a real-space cutoff of 9 Å and a grid spacing of 1 Å per grid point in each dimension. The van der Waals interactions were modeled with a Lennard-Jones potential, using also a cutoff of 9 Å. The simulations were performed in the NPT ensemble and the temperature was kept at 310 K by applying the Langevin thermostat with a collision frequency set to 5 ps -1 . The pressure was kept at 1.013 bar using a Berendsen barostat (5).
The force field parameters of both ligands with formal charge equal to zero were taken from the GAFF force field (6), and the missing ones were generated using the antechamber/parmchk2 modules of AMBER14. Namely, atomic restrained electrostatic potential (RESP) charges were derived using the antechamber tool of AMBER, after structural optimization at the b3lyp/6-31G(d,p) level performed with Gaussian09 (7) in the presence of implicit solvent (PCM) in order to avoid overstabilization of intramolecular H-bonds. The parameters of the compounds investigated here are available upon request to the authors. The parm14SB force field was used to parametrize the protein, in conjunction with the TIP3P (8) model for water and the Joung&Cheatam modified parameters for ions (9).

Post-processing of trajectories.
The free energy of binding of the inhibitors to the protein was evaluated by means of the Molecular Mechanics -Generalized Born Surface Area (MM-GBSA) postprocessing method (10, 11) using the MMPBSA.py tool of the AmberTools package (12). According to the MM-GBSA theory (13,14), the free energy of binding Gbind is evaluated through the formula: Gcom, Grec, and Glig are the absolute free energies of complex, receptor, and ligand, respectively, averaged over the equilibrium trajectory of the complex (single trajectory approach). According to these schemes, the free energy difference can be decomposed as: where ΔEMM is the difference in the molecular mechanics energy (null in the singletrajectory approach, within which the structures of the apo-protein and the ligand are extracted from the dynamics of the complex), ΔGsolv is the solvation free energy, and TΔSconf is the solute conformational entropy. The first two terms were calculated with the following equations: EMM includes the molecular mechanics energy contributed by the bonded (Ebond, Eangle, and Etorsion) and non-bonded (Evdw and Eele, calculated with no cutoff) terms of the force field. ΔGsolv is the solvation free energy, which can be modeled as the sum of an electrostatic contribution (ΔGsolv,p, evaluated using the MM-GBSA approach) and a non-polar one (ΔGsolv,np = γΔSA + b, proportional to the difference in solvent-exposed surface area, ΔSA).
In the MM-GBSA approach, the electrostatic solvation free energy was calculated using the implicit solvent model in Ref. (15,16) (igb = 8 option in AMBER14) in combination with mbondi3 (17, 18) (for H, C, N, O, S elements) and intrinsic (19) radii. Partial charges were taken from the AMBER/GAFF force fields, and relative dielectric constants of 1 for solute and 78.4 for the solvent (0.1 M KCl water solution) were used. The non-polar contribution is approximated by the LCPO (20) method (20) implemented within the sander module of AMBER. In addition to being faster, the MM-GBSA approach furnishes an intrinsically easy way of decomposing the free energy of binding into contributions from single atoms and residues (21), which is alternative to the "alanine scanning" approach (4,(21)(22)(23)(24). Solvation free energies were calculated on the optimized poses and on 80 frames extracted from the last 4 ns of the partly restrained MD simulation. The solute conformational entropy contribution (TΔSconf) is composed by a rototranslational term, calculated through classical statistical mechanics formulas, and by a vibrational term, which has been estimated here through normal-mode analysis using the mmpbsa_py_nabnmode module of AMBER14. Solute entropies were calculated, Concerning the other residues, only those contributing more than 1 kcal/mol to stabilization of the complex are reported.