Interdependence of Inhibitor Recognition in HIV-1 Protease.

Molecular recognition is a highly interdependent process. Subsite couplings within the active site of proteases are most often revealed through conditional amino acid preferences in substrate recognition. However, the potential effect of these couplings on inhibition and thus inhibitor design is largely unexplored. The present study examines the interdependency of subsites in HIV-1 protease using a focused library of protease inhibitors, to aid in future inhibitor design. Previously a series of darunavir (DRV) analogs was designed to systematically probe the S1' and S2' subsites. Co-crystal structures of these analogs with HIV-1 protease provide the ideal opportunity to probe subsite interdependency. All-atom molecular dynamics simulations starting from these structures were performed and systematically analyzed in terms of atomic fluctuations, intermolecular interactions, and water structure. These analyses reveal that the S1' subsite highly influences other subsites: the extension of the hydrophobic P1' moiety results in 1) reduced van der Waals contacts in the P2' subsite, 2) more variability in the hydrogen bond frequencies with catalytic residues and the flap water, and 3) changes in the occupancy of conserved water sites both proximal and distal to the active site. In addition, one of the monomers in this homodimeric enzyme has atomic fluctuations more highly correlated with DRV than the other monomer. These relationships intricately link the HIV-1 protease subsites and are critical to understanding molecular recognition and inhibitor binding. More broadly, the interdependency of subsite recognition within an active site requires consideration in the selection of chemical moieties in drug design; this strategy is in contrast to what is traditionally done with independent optimization of chemical moieties of an inhibitor.

Introduction

Human immunodeficiency virus type 1 (HIV-1) protease is a retroviral aspartyl protease that is an essential enzyme required for processing viral polyproteins and maturation of the virus and therefore a key therapeutic target. Highly active antiretroviral therapy (HAART), the current treatment standard, has significantly improved mortality and morbidity rates of patients infected with HIV-1.[1−5] HAART is a combination therapy consisting of three or more drugs from two or more classes. Protease inhibitors (PIs) have become a vital component of HAART and key to treatment of HIV-1 infections. The emergence of resistant viruses threatens the efficacy of current PIs and can lead to treatment failure.

Currently, there are eight FDA approved PIs. Darunavir (DRV), the latest PI approved by the FDA, is the most potent antiretroviral drug thanks to a high antiviral activity and high genetic barrier to the development of resistance (https://www.fda.gov/). Multiple mutations throughout the protease are needed to confer significant levels of resistance to DRV. Understanding the driving forces underlying the superior resistance profile of DRV compared to other PIs not only aids the future design of PIs but also due to the wealth of structural information HIV-1 protease is an excellent system to test general design principles that can be applied to other systems.

HIV-1 protease is a 99 amino acid homodimer (Figure A). The active site of HIV-1 protease can be characterized as a channel that has eight subsites (S4–S1 and S1–S4′). Each subsite position corresponds to an amino acid of the substrate (P4–P1 and P1′–P4′ from N to C terminus) with the scissile bond between the P1–P1′ positions.[6] DRV occupies four subsites (S2 to S2′), with P2, P1, P1′, and P2′, making contacts with hydrophobic residues and several aspartic acid residues including catalytic D25 and D25′ (Figure B). Because protease contains two identical monomers, by convention the monomer binding the C terminal side of substrates and containing subsites S1′ to S4′ is referred to as the prime monomer. The aniline moiety of DRV by analogy of peptidomimetics corresponds to P2′, while the bis-THC moiety is assigned to P2 (Figure B).

Figure 1

A) HIV-1 protease bound to DRV (PDB ID: 1T3R). Protease is a symmetric homodimer (monomers in purple and green) where each monomer is comprised of 99 amino acids but binds asymmetric substrates and inhibitors. By convention the residues and subsites within the monomer that binds the C-terminal side of peptide substrates or the aniline moiety of DRV are referred to as the prime residues to distinguish between the identical monomers. B) Chemical structure of DRV with P2 to P1 and P1′ to P2′ moieties labeled. Subsites S2 to S2′ correspond to the residues that make up the binding cleft for each ligand moiety P2 to P2′ (indicated by purple and light green as in panel A). HIV-1 protease active site residues that interact with DRV are depicted in green for hydrophobic, red for acidic, and gray for glycine.

The substrate envelope hypothesis has been shown to be an effective strategy to explain the molecular mechanism of drug resistance mutations and guide the design of new resistance-proof inhibitors.[7−14] Where inhibitors exceed the substrate envelope are probable sites for resistance mutations to occur, and inhibitors designed to fit within the volume that the endogenous substrates occupy are less prone to resistance.[15−18] DRV has been evaluated against the substrate envelope and fits well within this constrained volume.[12−14] Previous studies found that the isopropyl functional group at the P1′ position of DRV did not fully fill the substrate envelope, and thus larger P1′ moieties could be explored.[12] Further, free energy decomposition studies of HIV-1 protease mutants bound to DRV showed that replacement of the isopropyl (P1′) and aniline (P2′) groups could improve potency and the resistance profile.[19,20] Resistance mutations V82A and I84V increase the size of the active site because the substituted amino acids are smaller in size, and in response the P1′ and P2′ groups adopt an alternate conformation that leads to a loss of van der Waals contacts.

A series of DRV analogs was designed to optimize the DRV scaffold by altering the functional groups at the P1′ and P2′ positions while maintaining the constraints of the substrate envelope (Table ). Two functional groups, methyl-butyl and ethyl-butyl, were evaluated at the P1′ position with five different moieties at the P2′ position for a total of ten DRV analogs. Structure, enzymatic activity, antiviral activity, and pharmacokinetic properties of these inhibitors were evaluated with wild-type protease and resistant variants.[13] All inhibitors, like DRV, have low picomolar enzymatic potency and antiviral activity superior to DRV, and several inhibitors also exhibited good microsomal stability. Additionally, high-resolution X-ray crystal structures of the DRV series bound to wild-type HIV-1 protease have been determined, and the inhibitors were found to fit well within the substrate envelope as designed. Previous design and characterization of these DRV analogs (Table ) provide a unique opportunity to determine the coupled subsite relationships within the active site of HIV-1 protease in molecular recognition and inhibition.

Table 1

DRV and 10 Analogs Designed To Optimally Fill the Substrate Envelope and Evaluated in This Studya

The two series, UMASS1–5 and UMASS6–10, have the same R1 group at the P1′ position, and 5 congeneric pairs within the series have the same R2 group at the P2′ position.

There are numerous examples where protease subsites exhibit an interdependence in substrate specificity, i.e. subsites cannot be decoupled but must be optimized in a concerted fashion.[21−27] The interdependence of substrate specificity of subtilisin, a serine protease,[22,24] has been previously analyzed. Contributions to kcat/Km from each subsite of subtilisin were not additive; amino acids in one site shield adverse effects of less favorable amino acids at other distal positions. In Factor VIIa, another serine protease, the subsites P3 and P4 were found to be coupled in terms of substrate specificity as nonhomologous amino acid classes were preferred in the P3 subsite depending on the amino acid present in the P4 subsite.[23] In the protease elastase, the effects of contact complementary optimization in the P1 subsite were negated if the contacts in the P4 subsite were optimized independently.[25−27] Each of these examples, as well as many others, focuses on modifying or optimizing substrate specificity or processing. While the interdependence of substrate specificity is well-known, these relationships do not necessarily translate to design and optimization of resistance-proof inhibitors. Little has been done to evaluate interdependence of subsites in inhibitor binding. One example in HIV-1 protease showed a compensatory effect between the P1–P2 and P1′–P2′ subsites in the evolution of resistance mutations against known PIs.[28]

In the current study, we aimed to determine if there is an interdependence between subsites of HIV-1 protease in inhibitor binding, using the focused library of DRV analogs where the functional groups at P1′ and P2′ positions are systematically altered (Table ). The DRV analogs, along with the available high-resolution structural data, provided an ideal set for evaluation of subsite interdependence. Molecular Dynamics (MD) simulations starting from cocrystallized structures were performed to systemically analyze the atomic fluctuations, energetics, hydrophobic and electrostatic interactions, and water structure for HIV-1 protease bound to DRV and 10 DRV analogs and determine how different functional groups influence these properties. The inhibitor interactions at the P1′ position were found to influence the van der Waals (vdW) contacts in the P2′ position, as well as hydrogen bonding between the inhibitor and key structural features, namely the catalytic residue D25 and the highly conserved flap water. Thus, subsites in differential recognition of inhibitors within the active site of HIV-1 protease are highly coupled and interdependent.

Results

MD trajectories starting from high-resolution structures of HIV-1 protease cocrystallized with DRV and DRV analogs were generated and analyzed to reveal subsite interdependence and interactions between the inhibitor and protein. Each complex was simulated in triplicate for 100 ns each. Correlated fluctuations between the protein and inhibitor, vdW contact energies, electrostatic interactions, and water structure were evaluated and compared between DRV and the analogs to reveal subsite relationships.

Asymmetric Inhibitors Propagate Asymmetric Fluctuations Throughout the Enzyme

The asymmetry in B-factors between the monomers in the crystal structures of the inhibitor complexes indicated a link between the chemical moieties present in the substrate binding sites and the dynamics of the monomer (Figure S1 and analysis in the Supporting Information: Crystallographic B-factor comparison). To identify the regions of the homodimeric HIV-1 protease that fluctuate the most highly in concert with atomic fluctuations of the bound inhibitor cross-correlational analysis was used. Specific regions of the protease are highly correlated/anticorrelated with the atomic fluctuations of DRV (Figure ) and DRV analogs (Figures S2 and S3). Residues D25, T26, G27, G28, D29, G49, and I50 of the protease were the most correlated with the fluctuations of the inhibitor, but the correlations of each monomer with the inhibitors were asymmetric. The nonprime monomer containing the P2 pocket was always more highly correlated with the inhibitor motions (Figure B, red and orange colors), while the elbow region in the prime monomer was the most anticorrelated with the inhibitor (Figure B, blue). Distance difference plots also showed differential results for the two monomers as a result of the asymmetry of the inhibitor, confirming the results obtained from analysis of the crystal structures (analysis in the Supporting Information: Distance difference matrices show asymmetric inhibitors impact homodimeric protease monomers differentially and Figures S4–S7). Thus, the effects of the asymmetric inhibitor are propagated in an asymmetric manner to distal protein residues.

Figure 2

A) Pearson cross-correlations between DRV inhibitor atoms and C-alpha positions of HIV-1 protease residues. B) Average cross-correlation intensities by residue determined in panel A mapped onto the protease structure.

Alterations of P1′ Impact P2′ van der Waals Contacts but Not Vice Versa

The interdependency of subsites was investigated by evaluating how different functional groups at P1′ and P2′ positions of the inhibitor alter vdW contacts across subsites. By comparing DRV with UMASS1 and UMASS6, where the P1′ increases in size by one and then two methyl groups relative to DRV, respectively (Figure ), the interdependency between S1′ and the other subsites was evaluated. As the P1′ moiety increased in size, vdW contacts at the S1′ subsite became more favorable as expected, but while no change was observed at the S1 or S2 subsites, the corresponding contacts at S2′ became less favorable due to loss of vdW contacts (Figure ).

Figure 3

van der Waals contact energies for DRV, UMASS1, and UMASS6, with the same P2′ moiety but differing hydrophobic substitutions at P1′. Error bars indicate errors determined by block averaging over the concatenated three trajectories.

Specifically, with the addition of each carbon atom, extending the P1′ moiety, as more favorable contacts with L23, I84, V82, and I50′ were made (see Figure A) the vdW contacts decreased at P2′. These coupled changes occurred as the dihedral angle between the P1′ and P2′ moieties was altered, leading to less favorable contacts between the aniline group of the inhibitor and protease residues A28′, D29′, and D30′ (Figure B). This trend persisted when other congeneric pairs of DRV analogs with the same P2′ moiety were compared (Figure S8).

Figure 4

Overlay of the HIV-1 protease structure in complex with DRV, UMASS1, and UMASS6. A) van der Waals contacts (L23, P81, V82, I84, and I50′) with the P1′ moiety increase, and B) the dihedral angle (cyan circles) between the P1′ and P2′ groups shifts with a hydrophobic extension of the P1′ moiety, changing the contacts between inhibitor and residues A28′, D29′, and D30′.

Upon binding, the conformational entropy of the inhibitor decreases as the inhibitor’s degrees of freedom are reduced. An additional penalty will be paid if the inhibitor’s bound state deviates from the low energy unbound state. MD simulations of the free (unbound) inhibitors were run to compare the P1′-P2′ torsion angle of DRV and DRV analogs in the bound and unbound states. The torsion angle for DRV was the same in both the bound and unbound states, while there was a 2 to 18 degree perturbation between the bound and unbound states for all other DRV analogs (Table ). The hydrophobic extension of the P1′ moiety not only leads to less favorable vdW contact energy at the S2′ subsite but the torsion angle between the P1′ and P2′ moieties in the bound state deviates from the lowest energy state.

Table 2

Comparison of the Dihedral Angle of the P2′ Moiety in the Bound and Unbound States for DRV and DRV Analogs

	P2′ torsion (degrees)
compound	bound	unbound	Δ torsion
DRV	–85	–85	0
UMASS1	–81	–90	9
UMASS2	–85	–74	11
UMASS3	–79	–81	2
UMASS4	–84	–92	8
UMASS5	–87	–75	12
UMASS6	–75	–93	18
UMASS7	–75	–81	6
UMASS8	–68	–81	13
UMASS9	–84	–80	4
UMASS10	–84	–69	15

Next, the analogs within the two series, which have the same P1′ but a differing P2′ moiety, were compared. In the P1′ methyl-butyl (UMASS 1–5) and ethyl-butyl (UMASS 6–10) series, as the chemical diversity was introduced at the P2′ position there was no significant change in the vdW contacts between the protein and inhibitor P1, P2, or P1′ positions (Figure ). Therefore, the observation that the P1′ moiety impacting the vdW contacts at the P2′ position is unidirectional, as the reverse is not true.

Figure 5

van der Waals contact energies by subsite between DRV analogs and HIV-1 protease. Error bars are the errors determined by block averaging over the concatenated three trajectories.

Alterations in P1′ Also Impact Key Direct and Water Mediated Hydrogen Bonds

All intermolecular hydrogen bonds between DRV and HIV-1 protease were determined along with the percentage of the time the bond existed during the MD simulations (Figure S9). Direct protein-inhibitor hydrogen bonds were observed between DRV and residues D25, D29, D30, D25′, and D30′. Additionally, there were water-mediated hydrogen bonds with I50, I50′, and D29′. The hydrogen bonding patterns between protease and DRV were compared to other DRV analogs to analyze the effect of different chemical moieties in the S1′ and S2′ subsites.

The two series with differing P1′ (UMASS1–5, UMASS 6–10) were compared to evaluate the impact of the hydrophobic extension of the P1′ moiety on the hydrogen bonding patterns and frequencies in each subsite (Table S1). No significant differences were observed in the hydrogen bonds of S2 or S1′ subsite with the P1′ extension. For all P2′ congeneric pairs, with the exception of DRV, UMASS1 and UMASS6, there was also no dependency between the S1′ and S2′ subsites (Expanded analysis in the Supporting Information: Interdependency of electrostatic interactions between the S1′ and S2′ subsites mediated by P2′ moiety). Specifically, the frequency of the transition state hydrogen bond between the protonated catalytic residue D25 and the oxygen occupying the P1 position of the inhibitor (Figure A) is more stable when a methyl-butyl moiety is in the P1′ position regardless of the P2′ moiety. The P1 hydrogen bond frequencies differed by 14% across the methyl-butyl series and by 39% across the ethyl-butyl series (Figure B). Thus, the P1′ extension introduced variability in the frequency of the P1 hydrogen bond, suggesting a dependency between the S1 and S1′ subsites. This P1–P1′ dependency indicates that a hydrophobic P1′ extension reduces the electrostatic interactions in the adjacent P1 subsite.

Figure 6

A) The hydrogen bond between the inhibitor and D25. B) The change in the frequency of the hydrogen bond shown in panel A for each congeneric pair. C) Water mediated hydrogen bond between the inhibitor P2′ moiety and I50′. D) The change in frequency of the water mediated hydrogen bond shown in panel C for each congeneric pair. Frequencies are the sum of two flap water mediated hydrogen bonds with I50′ and the inhibitor.

Within each series (UMASS1–5, UMASS6–10) the impact of the P2′ position on other subsites was evaluated. The same hydrogen bonding patterns were observed for UMASS1 through UMASS5 and UMASS6 through UMASS10 within the S2, S1, and S1′ subsites but not the S2′. The P2′ capping group influenced the hydrogen bonding patterns and frequencies between the P2′ moiety and the S2′ subsite residues. Namely, the P2′ hydroxyl capping group of UMASS3 and UMASS8 interacted differently with D29′ and D30′ compared to other DRV analogs, forming three hydrogen bonds. The heteroatoms of the fused ring system of UMASS4, UMASS5, UMASS9, and UMASS10 formed only one hydrogen bond in the S2′ subsite, while UMASS2 and UMASS7, with a terminal methoxy, formed that same hydrogen bond but at a much lower frequency (19% and 14%, respectively). While the chemical functionality of the P2′ moiety greatly influenced the hydrogen bonding pattern within the S2′ subsite itself, there was little effect on other subsites.

Water-mediated hydrogen bond patterns and frequencies between the inhibitors and the protease were also examined. Inhibitors with a polar heteroatom in the para position at the P2′ position (DRV, UMASS1, UMASS3, UMASS6, and UMASS8) formed a hydrogen bond with D29′. This water mediated D29′ hydrogen bond is within the S2′ subsite and not unexpected since the P2′ moiety of these inhibitors is capable of hydrogen bonding. As expected, the tetra-coordinated flap water was observed for DRV (Figure S4) and all DRV analogs with high frequency (Table S2). This flap water is known to play a critical role in molecular recognition by HIV-1 protease. While the flap water mediated hydrogen bonding patterns were similar in the S1 and S1′ subsites across all UMASS analogs, the frequency of the hydrogen bond with I50′ (Figure C) was more variable within the ethyl-butyl series than the methyl-butyl series: 61% versus 17% change, respectively (Figure D). The variability in the hydrogen bonding frequency of the flap water suggests a dependency between the P2′ moiety and the water mediated hydrogen bond of the P1′ moiety with I50′.

To summarize, the evaluation of intermolecular hydrogen bonding interactions shows an interdependency between the P1′ position and hydrogen bonds in the S1 subsite. Specifically, extending the P1′ moiety led to a reduction in the frequency of hydrogen bonding with the catalytic D25. Further, flap water mediated hydrogen bonding was impacted by introducing chemical diversity in the S2′ subsite, depending on the P1′ moiety. Thus, key protease−inhibitor interactions are intimately linked to the chemical moiety in the S1′ subsite.

The P1′ Moiety Impacts the Surrounding Water Structure

Water structure is often ignored when evaluating molecular recognition and binding events. The congeners DRV, UMASS1, and UMASS6 provide an opportunity to investigate the impact of increasing hydrophobicity at the P1′ position on the structured water molecules around the protein-inhibitor complex. To elucidate the changes in the solvent surrounding the active site in response to variations in the P1′ moiety, the explicit water in the MD simulations for DRV, UMASS1, and UMASS6 was compared. The crystal structures of HIV-1 protease in complex with these three inhibitors were previously determined to 1.2, 1.95, and 1.5 Å resolution, respectively,[13] and showed significant variability in water content, rendering the analysis of the crystallographic waters unfeasible. The location and mobility of water molecules were evaluated in terms of occupancy and conservation of water sites during the simulations, and alterations in water sites were examined across the three inhibitors. While the number of water molecules in the proximity of the HIV-1 protease active site varied significantly across inhibitors, water sites with high occupancy were conserved in the crystal structures (Figure S10).

All water sites within 6 Å of the inhibitor were determined. A water site was defined as having at least twice the density of bulk water for the duration of the MD trajectories (see Materials and Methods for further details). For DRV, UMASS1, and UMASS6, 31, 33, and 32 water sites were identified, respectively. Twenty-seven sites were conserved between at least two of the inhibitor complexes. On average, all conserved water sites were occupied at least 50% of the time (Figure A). Water sites 4, 19, 23, and 25 had the highest variability in occupancy (Figure B); two of these water sites are within the active site, and two are adjacent to the active site (Figure C). For the congeners DRV, UMASS1, and UMASS6, 10 unique water sites were found (i.e., no equivalent was found in any other simulation). As observed with the conserved water sites, these unique sites were located both proximal and distal to the active site. This result suggests that not only do the changes in the P1′ position impact water sites directly contacting the active site but also the effect propagates to distal sites as well. Thus, the P1′ extension impacted the occupancy of the water sites both within and distal to the active site.

Figure 7

A) The occupancy of the 27 conserved water sites as determined over a 300 ns MD simulation of HIV-1 protease bound to DRV, UMASS1, and UMASS6. B) The standard deviation in occupancy of each water site shown in panel A. Sites 4, 19, 23, and 25 (ordered by highest occupancy) have the largest deviation in occupancy for the series DRV, UMASS1, and UMASS6, suggesting that these water sites are influenced by the P1′ moiety. C) Water sites with the largest standard deviations are either in direct contact with the ligand (sites 4 and 25) or at distal sites (sites 19 and 23). The DRV-bound HIV-1 protease structure (PDB ID: 1T3R) is displayed to orient the reader to the relative location of conserved water sites.

Discussion

While subsite couplings in substrate specificity have been well documented for various proteases, an understanding of how these interdependencies may translate to inhibitor binding has been lacking. Many active site couplings are not revealed from static structures. Because protein flexibility and conformational change are critical to understanding the binding mechanism of HIV-1 protease, MD simulations of a focused library of DRV and analogs bound to protease were performed and systematically analyzed to reveal the effect of P1′ and P2′ groups on the interactions of moieties in the other subsites. Two major findings from these analyses are that (1) the P1′ moiety at the S1′ subsite has extensive effects on other subsites and (2) asymmetry of the inhibitor impacts each monomer differently.

The S1′ subsite plays a central role, propagating changes to other subsites (Figure ). An increase in hydrophobic contacts of the P1′ moiety decreased contacts in the S2′ subsite, but the dependency was unidirectional. The P1′ moiety also impacted the frequency of the hydrogen bond between the inhibitor and catalytic residue D25, the hydrogen bonds mediated by the conserved flap water, and the occupancy of water sites both proximal and distal to the active site. The DRV analogs in the focused library used were originally designed to increase the vdW contacts at the S1′ site by extending the P1′ group and optimally filling the substrate envelope. However, coupling between S1′ and S2′ subsites and a reciprocal decrease of interactions at the S2′ substite indicate these moieties should be optimized simultaneously, in context with each other.

Figure 8

Interdependencies between subsites mapped onto the DRV structure.

The asymmetry of the inhibitors was propagated to alter dynamics throughout the protease, differentially impacting the two protease monomers. The most highly concerted atomic fluctuations between the inhibitors and protease were always in the 30s and 50s loops of the nonprime monomer, which interacts more strongly with the inhibitor compared to the prime monomer. Several of our own previous studies have also shown that asymmetric inhibitors differentially affect the two monomers and disrupt the symmetry of the HIV-1 protease and the resistance mutations that are present in both monomers.[29] MD simulations of wild-type and a drug resistant variant of HIV-1 protease bound to DRV revealed one side of the active site widening while the other side constricting. A common mechanism of drug resistance for HIV-1 protease may be through altered protein dynamics, shifting the ensemble and altering the flexibility of the enzyme.[30,31] In addition, mutations distal from the active site were found to asymmetrically disrupt the hydrogen-bonding network, which propagated changes through the two monomers differently to the active site.

The findings in this work are consistent with previous studies that found structural couplings between the P1–P2 and P1′–P2′ subsites[28] when optimizing amino acids at these positions for peptide-based inhibitor design. The current work extended the analysis to the whole inhibitor-protease-solvent system, as well as determining the underlying molecular mechanisms, and intermolecular interactions and energetics behind subsite couplings. Since previous studies revealed that altering the P1′ and P2′ moieties could improve inhibitor properties,[19] the focus was on the corresponding subsites. Future studies could include systematic evaluation of other subsites, different inhibitor pharmacophores, and extend to drug resistant variants to understand how/if subsite interdependency is altered.

The ability to design robust and potent therapeutics is limited by our understanding of mechanisms of inhibition. While subsite interdependency has been investigated by means of substrate optimization, little had been done on inhibitor binding. Elucidation of the couplings between the subsites of HIV-1 protease in binding of inhibitors not only will aid in future drug design but also will provide a proof of concept for application to other therapeutic targets. We have shown that the HIV-1 protease active site is comprised of many highly interdependent relationships that cannot be further deconstructed. These relationships can be used to guide the design of new inhibitors and provide a template to elucidate subsite interdependencies for other important therapeutic targets.

Materials and Methods

Protein Structure Preparation

Crystallographic structures of DRV and DRV analogs bound to HIV-1 protease were previously determined (PDB ID: 1T3R, 3O99, 3O9A, 3O9B, 3O9C, 3O9D, 3O9E, 3O9F, 3O9G, 3O9H, 3O9I).[13] All structures were prepared for simulation using the Protein Preparation Wizard from Schrödinger;[32,33] hydrogen atoms were added, multiple residue conformations were resolved, protonation states were determined, and any structural violations were resolved. All crystallographic waters were maintained, but crystallization buffers and salts were removed. The hydrogen bond network was optimized by exhaustive sampling of water orientations and flipping the side chain orientation of Asn, Gln, His, and terminal hydroxyl and thiol hydrogen as appropriate. The structures were minimized with heavy atoms restraints until an RMSD of 0.3 Å was achieved using the OPLS2005 force field (which was also used to parametrize the inhibitors).[34−36]

Molecular Dynamics Simulations

All MD simulations of the HIV-1 protease bound complex and free inhibitor were performed using Desmond[36−39] as implemented by Schrödinger using the OPLS2005 force field. Systems were prepared by solvating the structure in a cubic box that extended at least 10 Å beyond the nearest solute atom in all directions using the TIP3 solvation model.[40] Sodium chloride was added to the equivalent of 150 mM to simulate physiological conditions. The system was neutralized by adding counterions as needed (Na+ or Cl–).

A rigorous pre-equilibration model was employed as described elsewhere.[41] Briefly, a series of restrained minimization were performed to gradually relax the system. Initially all heavy solute atoms were restrained with a force constant of 1000 kcal mol–1 Å–2 for 10 steps with steepest descent followed by up to 2000 steps using the LBFG method to a convergence of 50 kcal mol–1 Å–2. Restraints were removed from side chains using LBFG for 5000 steps or until a convergence of 50 kcal mol–1 Å–2. The restraints on the backbone were gradually removed in procession using the following force constants 1000, 500, 250, 100, 50, 10, 1, and 0 kcal mol–1 Å–2 using the LBFG method to a convergence of 50 kcal mol–1 Å–2.

A series of short preproduction MD simulations were performed to equilibrate the system. A 10 ps NVT ensemble with restraints on solute heavy atoms of 50 kcal mol–1 Å–2 using a Berendsen thermostat at 10 K was performed. A 1 fs time step for bonded and short-range interaction up to 9 Å and a 3 fs time step for long-range electrostatic interactions were employed. A 10 ps MD simulation using an NPT ensemble with a Berendsen thermostat followed, run at 10 K with a 2 ps window for bonded and short-range interactions and 6 fs for long-range electrostatics. Over 50 ps, the temperature of the system was increased to 300 K with restraints on heavy solute atoms followed by a 10 ps simulation where all harmonic restraints were removed. Production simulations were performed using an NPT ensemble with a Nose-Hover thermostat and a Martyna-Tuckerman-Klein (MTK) barostat. Simulations were carried out with no harmonic restraints for 100 ns at 300 K and 1 bar. The cutoff for nonbonded interactions was 9 Å using Particle Mesh Ewald (PME) methodology. All simulations were performed in triplicate to check for reproducibility and ensure sufficient sampling, each with a different initial velocity for a total production time of 300 ns.

van der Waals Contact Analysis

The intermolecular van der Waals contact potential for each protease-inhibitor complex was calculated using a modified Lennard-Jones potential, V(r), equivalent to 4ε[(σ/r)12 – (σ/r)6] where r is the distance between atom pairs i of the protease and j of the inhibitor. The values of ε and σ are the depth of the energy well and the diameter of atoms when treated as a hard sphere. All protease-inhibitor nonbonded pairs within 6 Å were included. When the distance was less than ε the potential was set to ε.[42] The error associated with van der Waals energies was determined using block averaging over a total of 300 ns of concatenated trajectories.

Hydrogen Bonding Analysis

Identification and frequency of hydrogen bonds between protease and inhibitor (direct and water mediated) were determined using an in-house script built off of the Schrödinger API.[43] A hydrogen bond was determined to be present if the distance between the hydrogen and acceptor atoms was less than 2.5 Å. Additionally the following geometric restraints were required: 1) the angle between the donor-hydrogen-acceptor atoms had to be at least 120 degrees and 2) the angle between the hydrogen, acceptor, and atom bonded to the acceptor atom had to be at least 90 degrees. These criteria were applied to each frame of the trajectories to determine the frequency of the hydrogen bond. Where multiple hydrogen bonds were observed with a protein residue, such as with D25′, D29′, I50′, D30, and I50, the frequencies were added up and treated as a single electrostatic interaction.

Water Site Analysis

The spatially resolved water density from MD trajectories was calculated using an in-house Python script using a grid based approach, similar to one previously described.[44] Modules from the Schrödinger python API were utilized for trajectory handling and structural alignment. Each frame of the MD trajectory was fitted on a common reference frame. Thereafter the coordinates of the water oxygen atoms within 6 Å of the protein were mapped on a three-dimensional rectangular grid. The distance between the adjacent grid points was set to 0.5 Å. The grid spacing represents a compromise between maintaining good spatial resolution and a relatively low statistical error.

Following grid generation, local maxima in the discretized density distribution were identified. As an initial filter, only grid points with an average density of at least twice that of bulk water (0.066 Å–3) were considered. A grid point qualified as a water site if its density exceeded both this initial cutoff value and the density of all adjacent grid points. Finally, if two or more local maxima were located within 1.5 Å of each other, these grid points were considered to contribute to the same water site, and the coordinates of the water site were calculated based on the density weighted values of all contributing maxima.

The water-site occupancy was determined by comparing its coordinates with the coordinates of the water oxygen atoms in the prealigned trajectory. If the oxygen atom of a TIP3P water-molecule was located within 1.5 Å of the water-site, the molecule was considered to occupy the water site. If at any point during the course of the simulation more than one water-molecule was located within 1.5 Å, the water site was still considered to be occupied by only one water-molecule. For all subsequent calculations the water-molecule closer to the site was considered to be the ”occupying” water.

Cross-Correlations

Atomic positional fluctuations were determined for each simulation, and cross correlations between atom pairs were calculated usingwhere ΔR is the change in the position of site i, and ΔR is the change in the position in site j. The range of the correlation varies from minus one to one, minus one meaning the atom pair is completely anticorrelated, one meaning a complete correlation, and zero meaning that there is no correlation. Results were generated and plotted using in-house FORTRAN scripts and MatLab.