Improved prediction of HIV-1 protease-inhibitor binding energies by molecular dynamics simulations

Correlation between experimentally determined and calculated binding energy

The primary objective of this study was to determine whether protein relaxation would improve prediction of protein-ligand binding energies. Table 1 shows that the correlation coefficient of the experimentally determined and calculated binding energies from AutoDock for the twenty-five protease-inhibitor complexes after protein-rigid docking was 0.38. The correlations after 0.01, 0.1, 1, and 10 picoseconds (ps) of MD simulations were 0.53, 0.87, 0.79, and 0.74 respectively. These correlations are plotted in Figure 1. One outlier (from 1hvi) was noticed in Figure 1B; after it is eliminated, the correlation coefficient changes from 0.53 to 0.72 after 0.01 ps of simulation.

When a constant value of 6.5 is subtracted from the predicted energies after 0.1 ps MD, the binding energies of almost all of the predictions were within 2.0 kcal/mol of the experimental values. Three complexes, 1hvi, 1hvr, and 1hte had poor predicted energies, with an average error of 4.34, 4.31, and 5.09 kcal/mol, respectively. Among these, one (1hvi) had better predicted energies when compared to the results of protein-rigid docking.

Influence of the protease flap movement on calculated binding energy

The beta-strand flap is the most flexible region in the HIV-1 protease. It is normally 7 Å RMSD from the active site and is in an open conformation in the native state [22, 23]. The protease undergoes significant structural changes on binding to an inhibitor. The two flaps fold over the inhibitor to form a tunnel-shaped active site and are held in this close position by hydrogen bonding from Ile50 and Ile50' NH groups of the enzyme to a water molecule, which in turn is hydrogen bonded to the P2 and P1' CO groups of the inhibitor [24]. The bonding stabilizes the flaps in a closed position and inhibits the activities of the enzyme.

MD simulation has been used to study the movement of the flap region of HIV-1 protease with a ligand [25–30]. The flaps initially opened to an all-atom RMSD of 25 Å within 200 ps and became completely open at the end of a 10 ns simulation. In this study (Figure 2), the flaps opened up and moved away from the x-ray structure from 0.54 Å at 0.1 ps to 3.30 Å RMSD at 10 ps (the flap RMSD was calculated from residue 40 to 60 of each protein chain). These movements, after 0.1 ps of simulation, are inversely correlated with the quality of the binding energy prediction. As shown in Table 1, the correlation coefficient significantly decreased from 0.87 at 0.1 ps to 0.74 at 10 ps as the all-atom flap RMSD increased from 0.54 to 3.30 at 0.1 and 10 ps, respectively.

Complementarity between the ligand and the binding site is the basic concept behind ligand binding. This is manifest as steric complementarity, i.e. the shape of the ligand is mirrored in the shape of the binding site, allowing molecular interactions between two molecules [31]. MD simulations allow rearrangement of the protease side chain, especially on the active site surface, which improves the interacting surface complementarities of the complex. As shown in Table 2, after 0.1 ps, the time scale that produced the best correlation coefficient, the average all-atom RMSD of the complex was only 0.35 Å, while the average all-atom RMSD of the flap region was 0.54 Å.

Influence of MD simulation duration on ligand binding

A small N value indicates strong specificity of binding, with all of the solutions resembling one of only a small number of binding conformations and orientations. On the other hand, if N is large, the experiment indicates a low specificity of binding, since the solutions are composed of many different binding conformations or orientations.

In this study, the clustering result of ligands had an average N of 10.72 for docking without any MD, 8.36, 4.12, 5.88, and 7.16 for docking with 0.01, 0.1, 1, and 10 ps MD simulation, respectively. On average, the lower number of docking solutions in the cluster for all MD simulation time scales indicates that the ligands bind to their binding pocket with high specificity. The docked ligands after 10 ps MD simulation exhibited a wide range of RMSDs (0.87 Å to 3.86 Å) with an average RMSD of 2.52 Å, which indicates the failure of the ligands to recognize and specifically bind to the binding site with the protease flaps open.

Molecular Dynamics simulations

MD simulations of all complexes were carried out with the NAMD software version 2.5b1 [53] using the X-PLOR force field [54]. Missing atoms in 1heg, 1hte and 5hvp were added by using the guesscoord command in NAMD. The van der Waals, bond, angle, dihedral, and improper dihedral parameters for all the ligands were adopted from the Hetero-compound Information Centre-Uppsala HIC-Up http://xray.bmc.uu.se/hicup[55]. The water molecule under the flaps present in all complexes was included in all steps throughout this study except for 1hvr, where the structural water was removed in the preparation of the protease to be docked with the ligand (the ligand of 1hvr, which is a urea-based inhibitor, does not bind a water molecule in this position). The water molecules were added to the binding site of 1hvs because all the water molecules were absent in the experimental structure. Protein protonation states were modeled as in the related HIV-1 protease MD simulation study [56]. All protein residues were modeled in their charged state except for one of the two aspartic acid groups (Asp 25 and Asp 25') in the active site since previous studies [57–59] have shown that at least one of these two aspartic acids is protonated. We used a protonated Asp 25' and deprotonated Asp 25 for all HIV-1 protease-ligand complexes. The terminal residues of both monomers were also protonated (Pro 1, Pro 1', Phe 99 and Phe 99').

The complexes were immersed in a 20 Å radius sphere of TIP3-water using the program SOLVATE [60] to allow the protein-ligand complexes to relax in an aqueous environment. One hundred steps of energy minimization of the protein-ligand-water complex were initially performed, followed by 10 ps MD simulation at 300 K, with an atom-based shifted distance-dependent dielectric constant, ε = 4r; a switch function on van der Waals interaction, and a time step of 1 femtosecond (fs). The nonbonded interaction list was updated every 20 time steps. The van der Waals interactions were truncated at a distance of 12 Å. The simulations were repeated with five different starting seeds. The structures at 0.01 ps, 0.1 ps, 1 ps and 10 ps were recorded and processed in the docking step.

Preparation of protease-inhibitor complexes

To calculate the binding energy with AutoDock, we first prepared the protein complexes by separating each snapshot recorded from MD simulations into one file containing the protease and the water molecules, and one file containing only the ligand. Polar hydrogens were then added to the protein coordinates with the PROTONATE utility from AMBER [61]. AMBER united atom force field charges were assigned, and solvation parameters were added using the ADDSOL utility. The 3D affinity grid fields were created using the auxiliary program AutoGrid. The center of protein mass was chosen as the grid center. In this stage, the protein was embedded in the 3D grid and a probe atom was placed at each grid point. The affinity and electrostatic potential grid was calculated for each type of atom in the protease molecule. The number of grid points in x, y, z-axis was 60 × 110 × 60 with grid points separated by 0.375 Å.

Ligands that had a peptide-like N- or C-terminal end were assigned a charge. Hydrogen atoms were added to fill all empty valences, and Kollman united-atom charges [62] were also assigned to the ligands. Rotatable dihedrals in the ligands were assigned using the program AutoTors and were allowed to rotate freely. The nonpolar hydrogens were removed and the partial charges from these were added to the carbon that held the hydrogen. The atom type for the aromatic carbons was reassigned to be handled by the aromatic carbon grid map. These preparations were done for each ligand using the AutoTors module.

Automated docking

Docking calculations were carried out using AutoDock, version 3.0.5 [4]. Three binding energy terms were taken into account in the docking step: the van der Waals interaction represented as a Lennard-Jones 12–6 dispersion/repulsion term, the hydrogen bonding represented as a directional 12–10 term, and the Coulombic electrostatic potential.

Docking runs were performed using the Larmarckian genetic algorithm (LGA) [4] as previously described [63, 64] with some modifications of the docking parameters. The LGA describes the relationship between the protein and the ligand by the translation, orientation, and conformation of the ligand. These "state variables" are the ligand's genotype, and the resulting atomic coordinates together with the interaction and intermolecular energies are the ligand's phenotype. The environmental adaptation of the ligand's phenotype was reverse transcribed into its genotype and became heritable traits.

Docking began with a population of random ligand conformations in random orientations and at random translations. Each docking experiment was derived from 100 different runs that was set to terminate after a maximum of 1,500,000 energy evaluations or 27,000 generations, yielding 100 docked conformations. The population size was set to 50. The elitism number, the rate of gene mutation and the rate of gene crossover were 1, 0.02 and 0.8 respectively. A pseudo-Solis and Wets local search was then used to minimize energy of the population. The probability that docking solution in the population would undergo a local search was set to 0.06 and the constraint was set to a maximum of 300 iterations per search. The maximum number of successes or failures before changing the size of local search space (rho) were both set to 4. The starting conformations of the ligand were set to random positions. Translations were set to have a maximum limit of 2 Å/step and the orientation step size for the angular component and the torsions had a maximum limit at 50 degrees/step.

At the end of a docking job with multiple runs, AutoDock performed cluster analysis. Docking solutions with ligand all-atom RMSDs within 1.0 Å of each other were clustered together and ranked by the lowest energy representative. The lowest-energy solution of the lowest ligand all-atom RMSD cluster was accepted as the calculated binding energy.