Skip to main content

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



The accurate prediction of enzyme-substrate interaction energies is one of the major challenges in computational biology. This study describes the improvement of protein-ligand binding energy prediction by incorporating protein flexibility through the use of molecular dynamics (MD) simulations.


Docking experiments were undertaken using the program AutoDock for twenty-five HIV-1 protease-inhibitor complexes determined by x-ray crystallography. Protein-rigid docking without any dynamics produced a low correlation of 0.38 between the experimental and calculated binding energies. Correlations improved significantly for all time scales of MD simulations of the receptor-ligand complex. The highest correlation coefficient of 0.87 between the experimental and calculated energies was obtained after 0.1 picoseconds of dynamics simulation.


Our results indicate that relaxation of protein complexes by MD simulation is useful and necessary to obtain binding energies that are representative of the experimentally determined values.


The human immunodeficiency virus type 1 aspartic protease (HIV-1 PR) is an important enzyme due to its key role in viral maturation. Inactivation of the enzyme causes the production of immature, noninfectious viral particles. The enzyme therefore is an attractive target in anti-AIDS drug design, and the effect of binding various inhibitors on the protease structure is currently the focus of intensive research [13].

To obtain information about the position and energy of binding between an inhibitor and the corresponding protein, several automated docking programs have been developed [46]. Given recent improvements in search algorithms and energy functions, computational docking methods have become a valuable tool to probe the interaction between an enzyme and its inhibitors. These methods can contribute significantly to the understanding of structural and energetic basis of enzyme-substrate interactions [79].

Protein-ligand docking methods aim to predict the binding energy of the protein-ligand complex given the atomic coordinates. In such calculations, both the protein and ligand can be treated as rigid bodies [10, 11]; alternately, the ligand, the protein, or both molecules, can be completely or partially flexible [12, 13].

One advantage of incorporating flexibility is that it enables a search without bias introduced by the initial model. This bias normally influences both the orientation and conformation of the ligand in the active site, which usually represents a local minimum conformation [14]. More importantly, the lock and key concepts used to evaluate enzyme-substrate binding, in reality, refer to flexible locks and keys that are both in constant dynamic (thermal) motion [15]. For example, an analysis of conformational changes upon complex formation for a representative set of 39 pairs of ligand-free and ligand-bound structures showed that 50% of these proteins undergo substantial main-chain and side chain conformational changes when the ligand is bound [15].

Several techniques have been developed to predict the binding energies of HIV-1 protease-inhibitor complexes [1621]. Head RD et al., 1996 [16] used an approach that calculates physicochemical properties of the ligands and the receptor-ligand complexes to estimate the free energy of binding. The enthalpy of binding is calculated by molecular mechanics, while properties such as complementary hydrophobic surface area are used to estimate the entropy of binding through heuristics. Gohlke H et al., 2000 [17] developed DrugScore, a knowledge-based scoring function, to discriminate between well-docked ligand binding modes and those largely deviating from the native structure.

Schapira M et al., 1999 [20] used the finite difference Poisson-Boltzmann implementation of the electrostatic term, in conjunction with appropriate surface and entropy terms to predict the binding energy of 13 HIV-1 protease-inhibitor complexes. The predicted binding energies had a correlation coefficient of 0.66 with the experimental data. Recently, Österberg F et al., 2002 used AutoDock 3.0, a ligand flexible docking program, together with combining 21 protease structures into a single representative grid of interaction energies for incorporating the side chain motion into a docking simulation protocol [21]. The correlation coefficient between experimental and calculated binding energies produced by this technique was 0.79.

Since most current docking programs treat the receptor and/or ligand as rigid objects, this represents an important cause of failure to predict the correct binding enzyme-substrate energies [18]. We address this problem by performing Molecular Dynamic (MD) simulations on HIV-1 protease-inhibitor complexed and using the resulting structures to calculate the binding energies by AutoDock, a ligand flexible docking program.

Results and discussion

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.

Figure 1

Plots of experimentally-determined and calculated binding energies for the twenty-five HIV-1 protease-inhibitor complexes. The correlation coefficient improves from 0.38 without any MD simulation (A), to 0.53 at 0.01 ps of MD simulation (B), peaking at 0.87 after 0.1 ps (C), and dropping off to 0.79 at longer time scales of to 1 ps or more (D).

Table 1 Comparison of experimentally-determined and calculated binding energies for twenty-five HIV-1 protease-inhibitor complexes.

All these correlations represent a significant improvement over the protein-rigid docking results. The best results were obtained from the structures at the 0.1 ps MD simulation time point. These structures also have low ( 0.3 Å) average all-atom root mean square deviations (RMSD) relative to the experimental results. Table 2 shows the all-atom RMSD between each simulated complex and the corresponding protease-inhibitor x-ray structure. The average all-atom RMSD for the complexes increases from 0.18 Å at 0.01 ps to 2.92 Å at the end of 10 ps of simulation time. Similar results are consistently observed regardless of the five starting seeds used in the MD simulations.

Table 2 All-atom root mean square deviation (RMSD) of the protein-ligand complexes relative to their corresponding x-ray structures.

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 [2530]. 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.

Figure 2

Superposition of the Cα traces of part of the HIV-1 protease x-ray structure 4phv before (dark line) and after (light line) 10 ps of MD simulation bound to the inhibitor (space-fill). The flap region (above the inhibitor) moved away from the x-ray structure during the simulation, with all-atom RMSDs of 0.54 Å at 0.1 ps, 0.95 Å at 1 ps and 3.30 Å at 10 ps, respectively. Generally, after 0.1 ps, as the all-atom RMSD increased, the correlation coefficient of the experimentally-determined and calculated binding energies decreased.

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

Ligand docking revealed a consistent set of recurring binding modes. For all MD time scales, well-clustered docking results could be obtained. Generally, the lowest binding energy clusters are associated with the lowest all-atom RMSDs of the ligands. The best results in terms of lowest binding energy are located in a similar position of the x-ray structure at the active site. Table 3 shows the number of docking solutions in a cluster (N) along with the all-atom ligand RMSD for each MD time scale.

Table 3 Ligand all-atom RMSD (Å) and the number of docking solutions (N) in the cluster from 100 Larmarckian genetic algorithm (LGA) docking runs of twenty-five protease-inhibitor complexes.

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.


In this study, we illustrate the importance of taking dynamics into account to predict the structure and energetics of protein-ligand interactions. It is clear that relaxation of HIV-1 protease for 0.1 ps MD simulations is enough for rearrangement of the surface-binding pocket to produce good correlations between calculated and experimental binding energies. The binding energies of all protease molecules bound to different inhibitors were influenced by the movement of the flap regions with the correlation coefficient decreasing as the flaps moved away from the experimental structure. The differences in these correlations may reflect biological features of the dynamics of HIV-1 protease-inhibitor complexes.

Future work with larger data sets, different energy functions, different docking and binding energy evaluation methods, and more starting seeds, is necessary to determine the optimal parameters to robustly predict protein-substrate binding energies in silico.


A set of twenty-five HIV-1 protease-inhibitor complexes with experimentally determined structures and binding energies chosen from the Protein Data Bank (PDB), were used in this study [3250]. All complexes were determined by x-ray crystallography with resolution and R-factor less than 3.0 and 0.19, respectively. The PDB codes of the complexes along with the experimental binding energies are listed in Table 1.

The experimental binding energies used in this study were converted from the experimental inhibition constants (K i ), using the equation:

ΔG0 = -RT ln K i

where ΔG0 is the Gibbs free energy of binding (kcal/mol), R is the gas constant (1.987 cal K-1 mol-1) and T is the absolute temperature (room temperature, 300 K). Assuming that the K i values of each complex were obtained under the quasi-equilibrium conditions of Michaelis-Menton kinetics [51, 52], we expect that a plot of ΔG0 against ln K i is a straight line, with a slope given by the thermal energy, RT = 0.6 kcal/mol (at 300 K). Since the K i values were determined by different experimental techniques, some level of noise might exist in the underlying data. We expect our method to behave robustly in this regard.

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[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 [5759] 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.


  1. 1.

    Zhu J, Fan H, Liu H, Shi Y: Structure-based ligand design for flexible proteins: application of new F-DycoBlock. J Comput Aided Mol Des 2001, 15: 979–996. 10.1023/A:1014817911249

    CAS  Article  PubMed  Google Scholar 

  2. 2.

    Jayatilleke PR, Nair AC, Zauhar R, Welsh WJ: Computational studies on HIV-1 protease inhibitors: influence of calculated inhibitor-enzyme binding affinities on the statistical quality of 3D-QSAR CoMFA models. J Med Chem 2000, 43: 4446–4451. 10.1021/jm9905357

    CAS  Article  PubMed  Google Scholar 

  3. 3.

    Sham YY, Chu ZT, Tao H, Warshel A: Examining methods for calculations of binding free energies: LRA, LIE, PDLD-LRA, and PDLD/S-LRA calculations of ligands binding to an HIV protease. Proteins 2000, 39: 393–407. 10.1002/(SICI)1097-0134(20000601)39:4<393::AID-PROT120>3.3.CO;2-8

    CAS  Article  PubMed  Google Scholar 

  4. 4.

    Morris GM, Goodsell DS, Halliday RS, Huey R, Hart WE, Belew RK, Olson AJ: Automated Docking Using a Lamarckian Genetic Algorithm and Empirical Binding Free Energy Function. J Comput Chem 1998, 19: 1639–1662. Publisher Full Text 10.1002/(SICI)1096-987X(19981115)19:14<1639::AID-JCC10>3.0.CO;2-B

    CAS  Article  Google Scholar 

  5. 5.

    Henry AG, Richard MJ, Michael JES: Modelling Protein Docking using Shape Complimentarity, Electrostatics and Biochemical Information. J Mol Biol 1997, 272: 106–120. 10.1006/jmbi.1997.1203

    Article  Google Scholar 

  6. 6.

    Kramer B, Rarey M, Lengauer T: Evaluation of the FlexX incremental construction algorithm for protein-ligand docking proteins. Structure Functions and Genetics 1999, 37: 228–241. Publisher Full Text 10.1002/(SICI)1097-0134(19991101)37:2<228::AID-PROT8>3.0.CO;2-8

    CAS  Google Scholar 

  7. 7.

    Verkhivker GM, Rejto PA, Gehlhaar DK, Freer ST: Exploring the energy landscapes of molecular recognition by a genetic algorithm: 10 analysis of the requirements for robust docking of HIV-1 protease and FKBP-12 complexes. Proteins 1996, 25: 342–353. Publisher Full Text 10.1002/(SICI)1097-0134(199607)25:3<342::AID-PROT6>3.3.CO;2-3

    CAS  Article  PubMed  Google Scholar 

  8. 8.

    Sham YY, Chu ZT, Tao H, Warshel A: Examining methods for calculations of binding free energies: LRA, LIE, PDLD-LRA, and PDLD/S-LRA calculations of ligands binding to an HIV protease. Proteins 2000, 39: 393–407. 10.1002/(SICI)1097-0134(20000601)39:4<393::AID-PROT120>3.3.CO;2-8

    CAS  Article  PubMed  Google Scholar 

  9. 9.

    Pegg SC, Haresco JJ, Kuntz ID: A genetic algorithm for structure-based de novo design. J Comput Aided Mol Des 2001, 15: 911–933. 10.1023/A:1014389729000

    CAS  Article  PubMed  Google Scholar 

  10. 10.

    Vakser IA: Evaluation of GRAMM low-resolution docking methodology on the hemagglutinin-antibody complex. Proteins 1997, Suppl 1: 226–230. Publisher Full Text 10.1002/(SICI)1097-0134(1997)1+<226::AID-PROT31>3.3.CO;2-0

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Vakser IA, Matar OG, Lam CF: A systematic study of low-resolution recognition in protein – protein complexes. Proc Natl Acad Sci USA 1999, 96: 8477–8482. 10.1073/pnas.96.15.8477

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  12. 12.

    Schaffer L, Verkhivker GM: Predicting structural effects in HIV-1 protease mutant complexes with flexible ligand docking and protein side chain optimization. Proteins 1998, 33: 295–310. 10.1002/(SICI)1097-0134(19981101)33:2<295::AID-PROT12>3.0.CO;2-F

    CAS  Article  PubMed  Google Scholar 

  13. 13.

    Verkhivker GM, Rejto PA: A mean field model of ligand-protein interactions: implications for the structural assessment of human immunodeficiency virus type 1 protease complexes and receptor-specific binding. Proc Natl Acad Sci USA 1996, 93: 60–64. 10.1073/pnas.93.1.60

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  14. 14.

    Ota N, Agard DA: Binding mode prediction for a flexible ligand in a flexible pocket using multi-conformation simulated annealing pseudo crystallographic refinement. J Mol Biol 2001, 30: 607–617. 10.1006/jmbi.2001.5147

    Article  Google Scholar 

  15. 15.

    Betts MJ, Sternberg MJ: An analysis of conformational changes on protein-protein association: implications for predictive docking. Protein Eng 1999, 12: 271–283. 10.1093/protein/12.4.271

    CAS  Article  PubMed  Google Scholar 

  16. 16.

    Head RD, Smythe ML, Oprea TI, Waller CL, Green SM, Marshall GR: VALIDATE: A new method for the receptor-based prediction of binding affinities of novel ligands. J Am Chem Soc 1996, 118: 3959–3969. 10.1021/ja9539002

    CAS  Article  Google Scholar 

  17. 17.

    Gohlke H, Hendlich M, Klebe G: Knowledge-based Scoring Function to Predict protein-ligand interactions. J Mol Biol 2000, 295: 337–356. 10.1006/jmbi.1999.3371

    CAS  Article  PubMed  Google Scholar 

  18. 18.

    Akaho E, Morris GM, Goodsell DS, Wong D, Olson AJ: A study on docking mode of HIV protease and their inhibitors. J Chem Software 2001, 7: 103–114.

    CAS  Article  Google Scholar 

  19. 19.

    Klebe G: Recent developments in structure-based drug design. J Mol Med 2000, 78: 269–281. 10.1007/s001090000084

    CAS  Article  PubMed  Google Scholar 

  20. 20.

    Schapira M, Totrov M, Abagyan R: Prediction of the binding energy for small molecules, peptides and proteins. J Mol Recognit 1999, 12: 177–190. 10.1002/(SICI)1099-1352(199905/06)12:3<177::AID-JMR451>3.0.CO;2-Z

    CAS  Article  PubMed  Google Scholar 

  21. 21.

    Österberg F, Morris GM, Sanner MF, Olson AJ, Goodsell DS: Automated docking to multiple targets structures: Incorporation of protein mobility and structure water heterogeneity in AutoDock. Proteins 2002, 46: 34–40. 10.1002/prot.10028

    Article  PubMed  Google Scholar 

  22. 22.

    Harte WE Jr, Swaminathan S, Mansuri MM, Martin JC, Rosenberg IE, Beveridge DL: Domain communication in the dynamical structure of human immunodeficiency virus 1 protease. Proc Natl Acad Sci USA 1990, 87: 8864–8868.

    PubMed Central  CAS  Article  PubMed  Google Scholar 

  23. 23.

    Harte WE Jr, Swaminathan S, Beveridge DL: Molecular Dynamics of HIV-1 protease. Protiens 1992, 13: 175–194.

    CAS  Article  Google Scholar 

  24. 24.

    York DM, Darden TA, Pedersen LG, Anderson MW: Molecular dynamics simulation of HIV-1 protease in a crystalline environment and in solution. Biochemistry 1993, 32: 1443–1453.

    CAS  Article  PubMed  Google Scholar 

  25. 25.

    Navia MA, Fitzgerald PM, McKeever BM, Leu CT, Heimbach JC, Herber WK, Sigal IS, Darke PL, Springer JP: Three dimensional structure of aspartyl protease from human immunodeficiency virus HIV-1. Nature 1989, 337: 615–620. 10.1038/337615a0

    CAS  Article  PubMed  Google Scholar 

  26. 26.

    Wlodawer A, Miller M, Jaskolski M, Sathyanarayana BK, Baldwin E, Weber IT, Selk LM, Clawson L, Schneider J, Kent SB: Conserved folding in retroviral proteases: crystal structure of a synthetic HIV-1 protease. Science 1989, 245: 616–621.

    CAS  Article  PubMed  Google Scholar 

  27. 27.

    Miller M, Schneider J, Sathyanarayana BK, Toth MV, Marshall GR, Clawson L, Selk L, Kent SB, Wlodawer A: Structure of complex of synthetic HIV-1 protease with a substrate-based inhibitor at 2.3 Å resolution. Science 1989, 246: 1149–1152.

    CAS  Article  PubMed  Google Scholar 

  28. 28.

    Scott WR, Schiffer CA: Curling of flap tips in HIV-1 protease as a mechanism for substrate entry and tolerance of drug resistance. Structure Fold Des 2000, 8: 1259–1265. 10.1016/S0969-2126(00)00537-2

    CAS  Article  PubMed  Google Scholar 

  29. 29.

    Collins JR, Burt SK, Erickson JW: Flap opening in HIV-1 protease simulated by 'activated' molecular dynamics. Nat Struct Biol 1995, 2: 334–338.

    CAS  Article  PubMed  Google Scholar 

  30. 30.

    Collins JR, Burt SK, Erickson JW: Activated dynamics of flap opening in HIV-1 protease. Adv Exp Med Biol 1995, 362: 455–460.

    CAS  Article  PubMed  Google Scholar 

  31. 31.

    Rastall RA: Protein, structure and function: Ligand binding.[]

  32. 32.

    PDB code: Crystal structures of complexes of a peptidic inhibitor with wild-type and two mutant HIV-1 proteases. Biochemistry 1996, 35: 10627–10633. 10.1021/bi960481s

    Article  Google Scholar 

  33. 33.

    PDB code: A check on rational drug design: crystal structure of a complex of human immunodeficiency virus type 1 protease with a novel gamma-turn mimetic inhibitor. J Med Chem 1995, 38: 3246–3252.

    Article  Google Scholar 

  34. 34.

    PDB code: Comparative analysis of the X-ray structures of HIV-1 and HIV-2 proteases in complex with CGP 5 a novel pseudosymmetric inhibitor. Structure 3820, 3: 381–389.

    Google Scholar 

  35. 35.

    PDB code: Crystal structure of a complex of HIV-1 protease with a dihydroxyethylene-containing inhibitor: comparisons with molecular modeling. Protein Sci 1992, 1: 1061–1072.

    Article  Google Scholar 

  36. 36.

    PDB code: Crystal-structure of HIV-1 protease in complex with VX-478, A potent and orally bioavailable inhibitor of the enzyme. J Am Chem Soc 1995, 117: 1181–1189.

    Article  Google Scholar 

  37. 37.

    PDB code: X-ray crystallographic studies of a series of penicillin-derived asymmetric inhibitors of HIV-1 protease. Biochemistry 1994, 33: 8417–8427.

    Article  Google Scholar 

  38. 38.

    PDB code: 1hvi, 1hvj, 1hvk, 1hvl, Hosur MV, Bhat TN, Kempf DJ, Baldwin ET, Liu BS, Gulnik S, Wideburg NE, Norbeck DW, Appelt K, Erickson JW: Influence of stereochemistry on activity and binding modes for C2 symmetry-based diol inhibitors of HIV-1 protease. J Am Chem Soc 1994, 116: 847.

    Article  Google Scholar 

  39. 39.

    PDB code: Structural basis of drug resistance for the V82A mutant of HIV-1 proteinase. Nat Struct Biol 1995, 2: 244–249.

    Article  Google Scholar 

  40. 40.

    PDB code: An orally bioavailable HIV-1 protease inhibitor containing an imidazole-derived peptide bond replacement: crystallographic and pharmacokinetic analysis. Biochemistry 1994, 33: 11671–11677.

    Article  Google Scholar 

  41. 41.

    PDB code: A novel, picomolar inhibitor of human immunodeficiency virus type 1 protease. J Med Chem 1996, 39: 392–397. 10.1021/jm9507183

    Article  Google Scholar 

  42. 42.

    PDB code: Structure-based design of novel HIV protease inhibitors: carboxamide-containing 4-hydroxycoumarins and 4-hydroxy-2-pyrones as potent nonpeptidic inhibitors. J Med Chem 1995, 38: 3624–3637.

    Article  Google Scholar 

  43. 43.

    PDB code: Rational design of potent, bioavailable, nonpeptide cyclic ureas as HIV protease inhibitors. Science 1994, 263: 380–384.

    Article  Google Scholar 

  44. 44.

    PDB code: 4phv, Bone R, Vacca JP, Anderson PS, Holloway MK: X-Ray Crystal Structure of the HIV Protease Complex with L-700,417, an Inhibitor with Pseudo C2 Symmetry. J Am Chem Soc 1991, 113: 9382.

    Article  Google Scholar 

  45. 45.

    PDB code: The crystal structures at 2.2-A resolution of hydroxyethylene-based inhibitors bound to human immunodeficiency virus type 1 protease show that the inhibitors are present in two distinct orientations. J Biol Chem 1992, 267: 22770–22778.

    Google Scholar 

  46. 46.

    PDB code: Rational design, synthesis, and crystallographic analysis of a hydroxyethylene-based HIV-1 protease inhibitor containing a heterocyclic P1'-P2' amide bond isostere. J Med Chem 1994, 37: 3100–3107.

    Article  Google Scholar 

  47. 47.

    PDB code: Structure of complex of synthetic HIV-1 protease with a substrate-based inhibitor at 2.3 Å resolution. Science 1989, 246: 1149–1152.

    Article  Google Scholar 

  48. 48.

    PDB code: Crystallographic Analysis of a Complex between Human Immunodeficiency Virus Type 1 Protease and Acetyl-Pepstatin at 2.0 Angstroms Resolution. J Biol Chem 1990, 265: 14209–14219.

    Google Scholar 

  49. 49.

    PDB code: Structure at 2.5-A resolution of chemically synthesized human immunodeficiency virus type 1 protease complexed with a hydroxyethylene-based inhibitor. Biochemistry 1991, 30: 1600–1609.

    Article  Google Scholar 

  50. 50.

    PDB code: Design, Activity and 2.8 Angstroms Crystal Structure of a C2 Symmetric Inhibitor Complexed to HIV-1 Protease. Science 1990, 249: 527–533.

    Article  Google Scholar 

  51. 51.

    Stryer L: Biochemistry. WH Freeman and Co., New York 3 Edition 1999, 187–195.

    Google Scholar 

  52. 52.

    Creighton TE: Proteins: Structures and Molecular Properties. WH Freeman and Co., New York 2 Edition 1993, 386–389.

    Google Scholar 

  53. 53.

    Kale L, Skeel R, Bhandarkar M, Brunner R, Gursoy A, Krawetz N, Phillips J, Shinozaki A, Varadarajan K, Schulten K: NAMD2: Greater scalability for parallel molecular dynamics. J Comput Phys 1999, 151: 283–312. 10.1006/jcph.1999.6201

    CAS  Article  Google Scholar 

  54. 54.

    Brunger AT: X-PLOR version 3.1, A system for X-ray crystallography and NMR. Yale University Press, New Haven, CT 1992.

    Google Scholar 

  55. 55.

    Kleywegt GJ, Jones TA: Databases in protein crystallography. Acta Cryst 1998, D54: 1119–1131.

    CAS  Google Scholar 

  56. 56.

    Geller M, Miller M, Swanson SM, Maizel J: Analysis of the structure of HIV-1 protease complexed with a hexapeptide inhibitor. Part II: Molecular dynamic studies of the active site region. Proteins 1997, 27: 195–203. 10.1002/(SICI)1097-0134(199702)27:2<195::AID-PROT5>3.3.CO;2-1

    CAS  Article  PubMed  Google Scholar 

  57. 57.

    Piana S, Carloni P: Conformational flexibility of the catalytic Asp dyad in HIV-1 protease: an ab initio study on the free enzyme. Proteins: Struct Funct Genet 2000, 39: 26–36. 10.1002/(SICI)1097-0134(20000401)39:1<26::AID-PROT3>3.0.CO;2-N

    CAS  Article  Google Scholar 

  58. 58.

    Hyland LJ, Tomaszek TA Jr, Meek TD: Human immunodeficiency virus-1 protease 2. Use of pH rate studies and solvent kinetic isotope effects to elucidate details of chemical mechanism. Biochemistry 1991, 30: 8454–8463.

    CAS  Article  PubMed  Google Scholar 

  59. 59.

    Yamazaki T, Nicholson LK, Torchia DA, Wingfield P, Stahl SJ, Kaufman JD, Eyermann CJ, Hodge CN, Lam PYS, Ru Y, Jadhav PK, Chang CH, Weber PC: NMR and X-rays evidence that the HIV protease catalytic aspartyl groups are protonated in the complex formed by the protease and a non-peptide cyclic urea-based inhibitor. J Am Chem Soc 1994, 116: 10791–10792.

    CAS  Article  Google Scholar 

  60. 60.

    Grubmuller H, Heymann B, Tavan P: Ligand binding: molecular mechanics calculation of the streptavidin-biotin rupture force. Science 1996, 271: 997–999.

    CAS  Article  PubMed  Google Scholar 

  61. 61.

    Pearlman DA, Case DA, Caldwell JW, Ross WS, Cheatham TE, DeBolt S III, Ferguson D, Seibel G, Kollman P: AMBER, a package of computer programs for applying molecular mechanics, normal mode analysis, molecular dynamics and free energy calculations to simulate the structural and energetic properties of molecules. Comp Phys Commun 1995, 91: 1–41. 10.1016/0010-4655(95)00041-D

    CAS  Article  Google Scholar 

  62. 62.

    Weiner SJ, Kollman PA, Case DA, Singh UC, Ghio C, Alagona G, Profeta S, Weiner P: A new force field for molecular mechanical simulation of nucleic acids and proteins. J Am Chem Soc 1984, 106: 765–784.

    CAS  Article  Google Scholar 

  63. 63.

    Huang X, Xu L, Luo X, Fan K., Ji R, Pei G, Chen K, Jiang H: Elucidating the inhibiting mode of AHPBA derivatives against HIV-1 protease and building predictive 3D-QSAR models. J Med Chem 2002, 45: 333–343. 10.1021/jm0102710

    CAS  Article  PubMed  Google Scholar 

  64. 64.

    Rockey WM, Laederach A, Reilly PJ: Automated docking of α-(1>4)- and α-(1>6)-linked glucosyl trisaccharides and maltopentose 14 into the soybean β-amylase active site. Protein 2000, 40: 299–309. Publisher Full Text 10.1002/(SICI)1097-0134(20000801)40:2<299::AID-PROT100>3.0.CO;2-G

    CAS  Article  Google Scholar 

Download references


This work was supported in part by a Searle Scholar Award to Ram Samudrala. We thank members of the Samudrala group for helpful comments.

Author information



Corresponding author

Correspondence to Ram Samudrala.

Additional information

Author's contributions

EJ performed Molecular Dynamics simulations, docking, evaluated the results, and drafted the manuscript. RS helped with evaluation of the results produced, refining the manuscript, and provided intellectual guidance and mentorship. All authors read and approved the final manuscript.

Authors’ original submitted files for images

Rights and permissions

Reprints and Permissions

About this article

Cite this article

Jenwitheesuk, E., Samudrala, R. Improved prediction of HIV-1 protease-inhibitor binding energies by molecular dynamics simulations. BMC Struct Biol 3, 2 (2003).

Download citation


  • Molecular Dynamic Simulation
  • Root Mean Square Deviation
  • Calculated Binding Energy
  • Docking Solution
  • Flap Region