Systematic Molecular Dynamics, MMâPBSA, and Ab Initio Approaches

Article pubs.acs.org/JPCB

Systematic Molecular Dynamics, MM−PBSA, and Ab Initio Approaches to the Saquinavir Resistance Mechanism in HIV‑1 PR Due to 11 Double and Multiple Mutations Haralambos Tzoupis,*,†,‡ Georgios Leonis,*,† Aggelos Avramopoulos,† Thomas Mavromoustakos,‡ and Manthos G. Papadopoulos*,† †

Institute of Biology, Pharmaceutical Chemistry and Biotechnology, National Hellenic Research Foundation, 48 Vas. Constantinou Avenue, Athens 11635, Greece ‡ Department of Chemistry, National and Kapodistrian University of Athens, Panepistimioupolis Zographou, Athens 15771, Greece S Supporting Information *

ABSTRACT: Mutations in the human immunodeficiency virus (HIV) enable virus replication even when appropriate antiretroviral therapy is followed, thus leading to the emergence of drug resistance. In a previous work, we systematically examined seven single mutations that are associated with saquinavir (SQV) resistance in HIV-1 protease (Tzoupis, H.; Leonis, G.; Mavromoustakos, T.; Papadopoulos, M. G. J. Chem. Theory Comput. 2013, 9, 1754−1764). Herein, we extend our analysis, which includes seven double (G48VV82A, L10I-G48V, G48V-L90M, I84V-L90M, L10I-V82A, L10I-L63P, A71V-G73S) and four multiple (L10I-L63PA71V, L10I-G48V-V82A, G73S-I84V-L90M, L10I-L63PA71V-G73S-I84V-L90M) SQV−HIV-1 PR mutant complexes, in an attempt to generalize our findings and formulate the main elements of the SQV resistance mechanism in the protease. On the basis of molecular dynamics (MD), molecular mechanics Poisson−Boltzmann surface area (MM−PBSA), and ab initio computational approaches, we identified specific features that constitute the HIV-1 PR mechanism of resistance at the molecular level: the low flexibility of SQV in the binding cavity and the preservation of hydrogen bonding (HB) and van der Waals interactions between SQV and several active-site (Gly27/27′, Asp29/29′/30/30′, especially Asp25/25′) and flap (Ile50/50′, Gly48/48′) residues of the protease contribute significantly to efficient binding. The total enthalpy loss in all mutants is mostly due to the loss in enthalpy of the active-site region. Furthermore, it was observed that mutation accumulation may induce stabilization to SQV and to the flaps through enhanced HB interactions that lead to improved inhibition (e.g., accumulation of mutations in complexes containing L10I, G48V, L63P, I84V, or L90M single mutations). It was also concluded that permanent flap closure is obtained independently of mutations and SQV binding is mostly driven by van der Waals, nonpolar, and exchangeenergy contributions. Importantly, it was indicated that the optimal positioning of SQV and the structure of the binding cavity are tightly coupled, since small changes in geometry may affect the binding energy greatly. The results of our theoretical approaches are in agreement with experimental evidence and provide a reliable description of SQV resistance in HIV-1 PR.

I. INTRODUCTION One important feature characterizing infectious diseases is the emergence of drug-resistant pathogens. The continuous administration of drugs leads to mutations in the genomes of viruses and bacteria that confer resistance. Human immunodeficiency virus 1 (HIV-1) is a typical example of the emergence of drug resistance through mutations in various genes.1 Extensive clinical research has identified two major drug targets for HIV-1 treatment, namely, enzymes reverse transcriptase (RT) and viral protease (PR).2,3 Currently, there are 16 commercially available inhibitors for RT and 9 for HIV-1 PR.4 Modern treatment regimens for HIV-1 introduce a complement of RT and PR inhibitors.5,6 Despite the effectiveness of such combinatorial cocktails in fighting HIV-1 infection, their © 2014 American Chemical Society

action diminishes over time, mostly due to the appearance of drug-resistant viral strains.7,8 HIV-1 PR, being one of the most important drug targets, has been extensively studied. The enzyme belongs to the family of aspartic proteases and consists of two identical chains of 99 amino acids each.9 As with all members of this family, the active site is found in a binding cleft, confined by two entrancemodulating flaps. In HIV-1 PR, the two catalytic aspartic acids (Asp) are parts of two triads comprising Asp25/25′−Thr26/ 26′−Gly27/27′, while the flaps are defined by residues at Received: March 18, 2014 Revised: July 18, 2014 Published: July 18, 2014 9538

dx.doi.org/10.1021/jp502687q | J. Phys. Chem. B 2014, 118, 9538−9552

The Journal of Physical Chemistry B

Article

Figure 1. Map of mutations in HIV-1 PR related to the present study. SQV (green) is found in the binding cavity of the protease. Two water molecules have been implicated in mediating hydrogen bonds between SQV and Asp29/Arg8′ in several HIV-1 PR mutants.50

positions 44-55/44′-55′.10,11 Experimental and theoretical studies have revealed the role of the flap region in direct interactions with drugs that lead to the effective inhibition of HIV-1 PR.10−12 This region has also been implicated in modulating the entrance of the natural substrate or inhibitors into the binding cleft, due to its highly flexible nature in the apo (unbound) form of the enzyme.13,14 One reason for the emergence of drug-resistant strains can be attributed to the pliable nature of HIV-1 PR.15 Studies have shown that approximately one-third of PR’s amino acids can be mutated, and the enzyme still retains most of its functionality.16−18 The pattern for the mutations in HIV-1 PR involves the appearance of a single mutation that leads to low-level drug resistance. This is closely followed by the accumulation of other mutations that subsequently increase drug resistance.15,18,19 The high level of resistance conferred by multiple mutations is not necessarily associated with a specific protease inhibitor (PI) but may be observed upon binding to a variety of inhibitors.20−23 Another problem is the low oral availability of certain drugs, such as saquinavir (SQV).24,25 Thus, to effectively overcome this obstacle, high drug doses are administered to patients.25 However, this strategy leads not only to side effects such as lipodystrophy and hypertension but also to increased drug resistance.26 The importance of HIV-1 PR inhibition is extensively documented in the literature, in both computational27−29 and experimental studies.30−32 Clinical studies have mapped in great detail the different mutation sites.16,17,33−36 Consequently, they have matched different mutations that cause resistance against specific drugs. One of the first-generation protease inhibitors (PIs) is saquinavir, which has been shown to mimic the transition state of the cleavage site (Phe-Pro amide bond) in the natural substrate of the protease.24,37−41 Among the different mutations associated with SQV resistance, the most important is considered to be the one at positions 48/48′ (Gly to Val).39,42 These are located at the flaps of the enzyme, and

several other SQV-related mutations have been reported in the literature.43 It has been suggested by Lefebvre et al. that cross resistance in HIV-1 PR is developed due to the fact that different inhibitors occupy similar spaces and have several groups in common, which interact with particular residues.44 This hypothesis has been supported by computational studies by Alcaro et al., which state that amino acids such as Gly27/27′ and Asp29/29′ develop common interactions.45 It has been suggested from the superimposition of different inhibitors with the natural substrate of the enzyme that drug resistance mutations occur in positions where the inhibitor protrudes from the active-site cavity.25,30,46−48 These observations suggest that effective inhibition of the protease should focus on enhancing drug−enzyme interactions.49 Identifying common features in different mutant strains would offer insight into such interactions and thus assist in designing more efficient drugs.44 Following our previous efforts to investigate the SQV resistance pattern in seven single HIV-1 PR mutants,50 here we complement our studies by focusing on seven additional double mutants (G48V-V82A, L10I-G48V, G48V-L90M, I84VL90M, L10I-V82A, L10I-L63P, A71V-G73S) and four multiple mutants (L10I-L63P-A71V, L10I-G48V-V82A, G73S-I84VL90M, and L10I-L63P-A71V-G73S-I84V-L90M) in an attempt to obtain a thorough description of the SQV resistance mechanism in HIV-1 PR (Figure 1). The combination of various computational techniques, such as molecular dynamics (MD), molecular mechanics Poisson−Boltzmann surface area (MM−PBSA), and ab initio energy decomposition analysis of interactions offers insight into the dynamic behavior, principal interactions, and binding mode of SQV inside the protease. Importantly, experimental studies48,51−54 have consistently validated relevant theoretical results obtained from our group regarding binding and resistance in HIV-1 PR (i.e., previous studies including SQV in single mutants,50 four drugs in the I50V mutant,55 and darunavir,56 aliskiren,57 and fullerene 9539



Article

Table 1. Comparison of Average RMSD Values (in Å) for WT and Mutant HIV-1 PR Complexesa WT G48V-V82A L10I-G48V G48V-L90M I84V-L90M L10I-V82A L10I-L63P A71V-G73S L10I-L63PA71V L10I-G48VV82A G73S-I84VL90M L10I-L63PA71VG73S-I84VL90M

flap RMSDb

protease RMSDb

complex 0.58 1.22 1.30 1.07 1.04 1.25 1.09 1.52 0.98

[0.54/1.23]a [0.54/1.11] [0.54/0.80] [0.92/0.80] [1.11/1.23] [1.11/1.47] [1.55/1.12] [1.11/1.47/1.55]

0.52 0.99 1.37 1.19 0.94 1.23 0.92 1.66 0.94

active site RMSDb 0.60 0.61 0.62 0.62 0.91 0.50 0.50 0.63 0.49

[0.55/0.97] [0.55/0.99] [0.55/0.65] [0.86/0.65] [0.99/0.97] [0.99/1.39] [1.45/0.68] [0.99/1.39/1.45]

[0.50/0.90] [0.50/0.49] [0.50/0.51] [0.66/0.51] [0.49/0.90] [0.49/0.53] [0.62/0.70] [0.49/0.53/0.62]

SQV RMSDb 1.15 2.00 1.44 1.24 1.01 0.97 0.93 1.80 0.77

[1.43/0.96] [1.43/0.99] [1.43/1.09] [0.92/1.09] [0.99/0.96] [0.99/1.28] [1.55/0.99] [0.99/1.28/1.55]

1.17 [0.54/1.11/1.23]

0.97 [0.55/0.99/0.97]

0.57 [0.50/0.49/0.90]

1.34 [1.43/0.99/0.96]

1.03 [1.12/0.92/0.80]

1.09 [0.68/0.86/0.65]

0.85 [0.70/0.66/0.51]

0.82 [0.99/0.92/1.09]

1.35 [1.11/1.47/1.55/1.12/0.92/0.80]

1.19 [0.99/1.39/1.45/0.68/0.86/0.65]

0.78 [0.49/0.53/0.62/0.70/0.66/0.51]

1.50 [0.99/1.28/1.55/0.99/0.92/1.09]

a

Flap region: residues 44−55/44′−55′. Active site: residues 21−31/21′−31′. bThe average RMSD values for the respective single mutants are provided in brackets.

analogs58 in WT HIV-1 PR). Therefore, it is shown that such computations reliably complement the experimental evidence by providing valuable information. This work is a contribution toward the complete description of SQV resistance in HIV-1 PR at the molecular level.

the last 5 ns of each drug−protease trajectory. The total binding free energy (ΔGbind) is computed using the following equation

II. METHODS II.1. Systems. As the initial template, the high-resolution (1.16 Å) crystal structure of the WT HIV-1 PR−saquinavir complex (PDB code: 3OXC) was used.32 According to HIV-1 PR treatment in our previous studies, the active site has been considered to be monoprotonated at position 25.54,56 II.2. Molecular Dynamics Simulations in Water. MD simulations in explicit water have been carried out with the SANDER module of AMBER 11.59 For the construction of the protease, the ff99SB force field has been used.60 Missing hydrogen atoms were added to SQV with Chimera 1.7.61 For the construction of the mutant complexes, we replaced the sidechain atoms of the amino acids at the respective positions with the LEaP program in AMBER, and next we visually inspected the areas where changes were made to exclude the possibility of unwanted steric clashes. The parameters for the drug were described with the general Amber force field (GAFF), and the ANTECHAMBER routine was used to assign AM1-BCC partial atomic charges.62,63 The solvation of the complexes in water was performed using the TIP3P water model.64 The MD simulation in all systems was performed at 300 K and for 20 ns in the NPT ensemble, with Langevin dynamics temperature scaling.65 Additionally, five 150 ns MD runs have been carried out for selected mutants. All hydrogen atoms were constrained to their equilibrium distance,66 and the particle mesh Ewald (PME) method was used to calculate long-range electrostatic interactions.67 For electrostatic and van der Waals interactions, a cutoff of 10 Å was applied. For a detailed description of the MD protocol, see the Supporting Information. Clustering calculations were performed with the ptraj program of AMBER. The hierarchical method68 was used for grouping, with the RMS as the distance metric (cutoff 2.5 Å). II.3. Molecular Mechanics Poisson−Boltzmann Surface Area Calculations. The calculation of the binding energy was performed with the MM−PBSA method69 during

where Gcomplex, GHIV‑1 PR, and GSQV are the energies for the complex, the protease, and the drug, respectively. The binding energy is expressed as a combination of enthalpy and entropy terms:

ΔG bind = Gcomplex − (G HIV ‐ 1PR + GSQV )

ΔG bind = ΔH − T ΔS

(1)

(2)

The enthalpy of binding can be further decomposed into protein−ligand and solvation free energy contributions (Supporting Information). The entropy term −TΔS (eq 2) was calculated with normalmode analysis using the NMODE module70,71 in AMBER. A complete description of the MM−PBSA scheme is presented in the Supporting Information. The binding energy is also reported in terms of the inhibition constant, Ki. Ki is an equilibration constant that describes the tendency of a compound to bind to a target molecule (e.g., protein). Specifically, the inhibition constant is the inhibitor concentration that is required to produce half maximal inhibition. Therefore, the smaller the inhibition constant, the smaller the amount of drug required to inhibit the activity of the protein. The conversion from ΔG to Ki was based on the equation ΔG = RT ln Ki.72,73 II.4. Ab Initio Computation of Intermolecular Interaction Energies. We have computed the interaction energy, ΔE, of certain pilot systems employing the method of Su and Li.74 ΔE is given by ΔE = E el + E ex + Erep + Epl + E disp

(3)

el

where E is the electrostatic term which involves the classical Coulomb interactions between monomers, Eex is the exchange contribution which contains the exchange terms between monomers, Erep is the repulsion term, and Epl is the polarization component defined as the orbital relaxation energy from the monomer Hartree−Fock (HF) spin orbitals to those of the supermolecule.74 The dispersion energy, Edisp, corresponds to 9540



Article

Figure 2. RMSD for SQV with respect to its initial conformation in selected mutant complexes. A description of the remaining complexes is provided in Figure S2. For simplicity, only the average value is shown for WT.

was observed that structural changes for SQV in HIV-1 PR are comparable for all combinations, regardless of single, double, or multiple mutations. However, six of the (double and multiple) mutations considered herein induce greater conformational changes in the structure of SQV (in comparison to the change observed in WT), while in another five mutants the drug undergoes smaller changes (Table 1). All RMSD values have been calculated with respect to the initial SQV conformations and fall within a reasonable range: 0.77 Å < RMSD < 2.00 Å. Interestingly, triple mutant L10I-L63P-A71V (avg. RMSD = 0.77 Å, Figure 2) has a smaller average RMSD than each of its single mutations alone. Moreover, beginning with L10I-L63P (avg. RMSD = 0.93 Å), the third added mutation, A71V, alters the conformation of SQV to less of an extent than the double mutant does. This example shows that the accumulation of mutations may decrease their effect on the structure of SQV. Similarly, SQV in L90M is more mobile than in G48V-L90M and I84V-L90M (Figure S1), thus indicating the formation of interactions that stabilize the drug in the double mutants. Of course, since each of the three mutations confers resistance to SQV, this is a mere indication regarding the drug’s conformational diversity among mutants. There may be cases, however, such as for G48V-V82A, where the RMSD of SQV is greater than its respective value in each single mutant. Also, the increased mobility of SQV in L10I-L63P-A71V-G73S-I84VL90M (as expressed by the broad range of RMSD values in Figure 2) may have induced flexibility in the protease that leads to the destabilization of interactions between SQV and HIV-1 PR.

the difference between the MP2 and HF interaction energies. The decomposition of ΔE identifies significant contributions and reveals their physical origins. The method of Su and Li is an extension and modification of techniques proposed by Kitaura and Morokuma,75,76 Ziegler and Rauk,77 and Hayes and Stone.78 At the HF level, the Eel term is the same as that computed in the Kitaura−Morokuma (KM) energy decomposition scheme.75,76,79 Also, the sum of Eex and Erep is identical to the exchange−repulsion (Exr) term in the KM approach, while the Epl term is the sum of the polarization, the charge transfer (Ect), and the mixing terms of the KM analysis. The decomposition analysis was carried out with the GAMESS80 quantum chemistry code, employing the 6-31G* basis set for all atoms. The basis set superposition error (BSSE) has been corrected using the counterpoise method of Bernardi and Boys.81

III. RESULTS AND DISCUSSION In this section, a comprehensive structural analysis of each protease complex, including important regions such as the active site and the flaps along with SQV in the binding cavity, has been attempted to elucidate the dynamics associated with specific mutations. Next, the investigation of hydrogen bond interactions will be complemented with ab initio and MM− PBSA energy analyses to rationalize the structural information and to provide a thorough description of the SQV resistance mechanism in HIV-1 PR. III.1. Conformational Analysis of Saquinavir in the Binding Site of WT and Mutated Forms of HIV-1 PR. It 9541



Article

Previous considerations suggested that conformational changes of the drug may not reduce binding, provided that the structure of SQV eventually stabilizes (as in V82A and A71V single mutants);50 similar observations involved six inhibitors in I50V complexes.55 In the present work, the drug’s conformations in double and multiple mutant strains vary from stable to flexible, thus indicating weakened binding in several cases (Figure S2). The RMSD profiles of L10I-V82A and A71V-G73S indicate that 20 ns simulations may not be sufficient to sample the conformational space adequately. For instance, in A71V-G73S a conformational change of the drug at ∼11 ns (Figure 2) indicates a rearrangement of SQV in the double mutant. To question this, we prolonged each MD run for 20 ns and also initiated one more simulation (with different random seeds) for each mutant. The results (Figure S3) demonstrated that the RMSD profiles after the initial ∼10 ns converge to similar values. Additionally, two different MD runs50,55 for the same SQV-WT HIV-1 PR complex show that these simulations also present similar profiles (Figure S3). In order to further ascertain the stability of the systems and eliminate their dependence on the initial conditions, we performed five additional MD runs (150 ns each) for the following complexes of SQV: WT, L10I-V82A, G73S-I84V-L90M, L10I-G48VV82A, and L10I-L63P-A71V-G73S-I84V-L90M. It was revealed that the mutated systems are practically stabilized after the first 15 ns (Figures S4 and S5). This observation suggests that the protease retains a relatively stable conformation early in each of the different MD simulations, despite the accumulation of mutations. III.1.1. Clustering Analysis. To enrich our observations, we performed clustering for SQV inside the binding cavity. In each single mutant, SQV showed only one representative structure throughout the simulation (Figure S6). This feature is also observed in the double and multiple mutants. The representative SQV structures were visualized in Chimera 1.7,61 and the different conformations were superimposed with the respective representative structure of SQV in WT (Figure 3,

Table 2. RSMD Values between the Representative SQV Conformations after Clustering and the Conformation of the Drug in WT−HIV-1 PR complexes L10I-G48V G48V-L90M G48V-V82A L10I-L63P L10I-V82A I84V-L90M A71V-G73S G73S-I84V-L90M L10I-G48V-V82A L10I-L63P-A71V L10I-L63P-A71V-G73S-I84V-L90M

SQV RMSD compared to WT (Å)a 2.10 2.02 4.30 2.01 1.26 1.37 3.94 1.54 4.17 1.75 4.32

[3.16/1.67] [3.16/1.26] [3.16/4.21] [1.67/4.03] [1.67/4.21] [1.29/1.26] [3.99/1.73] [1.73/1.29/1.26] [1.67/3.16/4.21] [1.67/4.03/3.99] [1.67/4.03/3.99/1.73/1.29/1.26]

a

Average RMSD values for the respective single mutants are given in brackets.

WT (avg. RMSD = 4.3 Å) and group b resembling the structure of the drug in WT (avg. RMSD = 1.6 Å). The most profound differences in SQV conformations are found in two mutants containing the G48V-V82A combination (namely, G48V-V82A and L10I-G48V-V82A) as well as in L10I-L63P-A71V-G73S-I84V-L90M. These conformational changes are mirrored in the RMSD values between WT and the mutants for the drug, as shown in Table 2. Similar RMSD values have been calculated for V82A (4.21 Å) and for G48V-V82A (4.30 Å), indicating that the V82A mutation is responsible for the greatest effect on the SQV conformation. Also, for mutations at positions 71/71′ the same pattern was observed between A71V and A71V-G73S. Representative SQV structures before and after the sudden conformational change observed above (complex A71S-G73S, Figure 2) are shown in Figure S3f. The RMSD between them is 1.9 Å, less than the cutoff for the clustering calculations (2.5 Å). The structural jump at 11 ns may be attributed mostly to a repositioning of the phenyl and quinoline rings. Additionally, the RSMD of SQV in the complex with L10I-G48V-V82A (4.17 Å) is close to the RMSD in V82A, as opposed to G48V and L10I single mutations (Table 2). This, in conjunction with the above observation, further supports the finding that V82A has a greater effect on the drug’s conformation than the changes in positions 48/48′ or 10/10′. III.1.2. Strain Energies. We have optimized the representative structure from the clustering of SQV by employing the B3LYP/6-31G* method and the polarizable continuum model (PCM) and conductor-PCM (C-PCM) solvent methodologies.83,84 The energy of this structure is denoted by E1. We have also computed (with B3LYP/6-31G*) the energy of the SQV conformation (E2), without optimization, from each double and multiple mutant. The difference E2 − E1 is defined as the strain energy (Es).85,86 Es presents the effect of the environment (protein and water) on the conformation of the ligand (Table 3).86 Consequently, the energy of the ligand in the cavity of the protein could deviate from its local (or global) minimum.85,87 The results of Table 3 show that the solvent and the interactions between SQV and the protease in WT have the most profound effect on the restriction of the drug’s conformation inside HIV-1 PR. This leads to the assumption that SQV inside WT may form several HBs and presents favorable interactions with different amino acids of the

Figure 3. Representations of the two major conformational groups (a and b) for SQV in the active site of HIV-1 PR. In red, the conformation of SQV in WT is shown; superimposed structures (blue) represent the SQV conformations in the different multiple mutant strains.

red). The matchmaker function82 of Chimera was used to compare the different conformations. The RMSD values from these comparisons are presented in Table 2. Thus, the representative SQV structures were divided into two main groups (Figure 3a,b), with group a deviating mostly from the 9542



Article

III.2. Conformational Analysis of the Protease. RMSD calculations for the Cα atoms of HIV-1 PR in all mutant complexes showed that the structure of the protease is eventually stabilized during the simulations (Figure S7). A summary of the results is presented in Table 1. We observed that all mutants involving two, three, or six mutations have RMSD values which are larger than that of WT. Similarly, as shown in our previous studies, the WT protease structures are more stable than the structures of any single mutant bound to either SQV50 or various inhibitors, as in the case of I50V with six drugs.55 Also, greater structural changes in HIV-1 PR were observed for the double, triple, and sextuple mutants than for individual single mutations. However, some exceptions denote that the accumulation of mutations does not necessarily produce incremental effects: in G48V-V82A, L10I-L63P, L10I-G48V-V82A, and L10I-L63P-A71V mutants, the RMSDs for the proteases are generally reduced compared to those for individual mutations V82A, L63P, and A71V (Figure S7). Furthermore, point mutations V82A, L63P, and A71V present significantly greater variations in RMSD values, thus inducing greater mobility in the protease than their aforementioned double and triple mutants. This observation could be partially rationalized, since positions 82/82′ are located close to the interface of the two monomers and, as has been suggested, mutations at this site interrupt the interactions between the two monomers and therefore cause a loss in the functionality of the enzyme.19

Table 3. Strain Energies (Es) of SQV (kcal/mol) complexes

strain energya

L10I-G48V G48V-L90M G48V-V82A L10I-L63P L10I-V82A I84V-L90M A71V-G73S G73S-I84V-L90M L10I-G48V-V82A L10I-L63P-A71V sextuple mutant WT−SQV

97.263 (96.700)b 92.243 79.066 89.106 86.596 75.928 97.263 89.106 102.284 (101.790)b 97.263 93.498 113.579

a

Values were computed with the PCM model; solvent: water. bCPCM model; solvent: water.

protease. The lower strain energy in the mutated complexes suggests that the drug is not as tightly bound; therefore, it may explain the increased resistance of the mutated complexes compared to WT. Additionally, we performed C-PCM calculations for two mutants (L10I-G48V and L10I-G48VV82A) to account for any differences between methodologies. We note that the difference between the PCM and C-PCM approaches has a minor effect on the computed strain energies (Table 3).

Figure 4. (a) RMSD based on Cα atoms for the flap region of HIV-1 PR in mutant complexes L10I-L63P and L10I-L63P-A71V and a comparison with the respective single mutants. (b) Atomic fluctuations for protease residues in HIV-1 PR complexes. 9543



Article

Figure 5. Evolution of the main hydrogen bonds in SQV−HIV-1 PR complexes.

of a fullerene derivative, which presented poor binding to HIV1 PR.55 Similarly, flap closing occurs in the double, triple, and sextuple complexes, therefore further supporting the aforementioned statement. III.4. Atomic Fluctuations. The increased mobility of active-site and flap residues in the single-mutation forms of the complexes is generally diminished upon accumulation of mutations. The active site is significantly stabilized in all double and multiple strains, except in the case of mutations containing L90M, where the active site remains practically unchangeable. The flap region is very stable in most double and multiple mutants yet is more flexible than the active site (Figure 4b and Table S1). A noticeable exception to the stabilization of the flaps is L10I-G48V, where the residues of the double mutants show greater fluctuations than the corresponding residues of L10I. It was mentioned earlier that mutants such as L10I-L63P and L10I-L63P-A71V present lower total and flap RMS deviations than single mutations L10I, L63P, and A71V. This trend is also reflected in the mobility of the corresponding residues (Figure 4). III.5. Hydrogen-Bonding Interactions. Saquinavir forms an extensive HB network with the binding cavity of WT HIV-1 PR, involving five residues of the protease in several HB interactions (Table S2 and Figure 5). HBs in the mutant strains involve fewer amino acids than in WT, thus supporting the observation that a significant number of residues must participate in drug interactions to maintain efficient binding. The importance of specific residues in the interactions between WT and SQV has also been reported in previous theoretical studies.29,53 As with the single mutations,50 the amino acids of double and multiple mutant strains that are associated with the main HB interactions are Arg8, Gly27/27′, Asp29/29′, Asp30/30′, Ile50/ 50′, and Gly(Val)48/48′. This observation is in accordance with other computational studies reported in the literature.27,45

Finally, despite the fact that mutation accumulation may produce smaller conformational changes in certain double and triple mutants, it is suggested that a significant number of mutations (e.g., L10I-L63P-A71V-G73S-I84V-L90M) would certainly increase the mobility of the complex and induce great conformational changes. In Figure S7, the greater changes and mobility of the sextuple mutant compared to those of the respective triple mutants are depicted by the RMSD plots. III.3. Conformational Analysis of the Active Site and Flap Regions. The active site (residues 21/21′−31/31′) of HIV-1 PR is relatively stable in all double and multiple mutant complexes (Table 1). However, in I84V-L90M, G73S-I84VL90M, and L10I-L63P-A71V-G73S-I84V-L90M a noticeable structural rearrangement of the binding cavity (RMSD ≥ 0.8 Å) was observed. Interestingly, we note that the I84V-L90M motif is present in these three cases. Also, the aforementioned observation regarding the greater structural changes on the V82A protease compared to those on double or multiple (V82A-containing) mutants is verified by the RMSD calculations for the active site. Mutation accumulation on the V82A strain induced a significant stabilization in the active-site region (Table 1). The results of Table 1 for the flap region (residues 44/44′− 55/55′) suggested the following: (i) All mutants have average RMSD values larger than the corresponding value of WT. This shows that the flaps undergo more significant structural changes in the mutant forms than in WT. (ii) L10I-L63P and L10IL63P-A71V have average RMSD values which are smaller than those of single mutants (i.e., L10I, L63P), whereas the flaps of complexes with single mutations are generally less affected than the flaps of double and multiple mutants (Figure 4). In our previous work, we observed that all single mutants induced a closed conformation in the flaps,50 thus suggesting that a closed structure does not necessarily imply efficient binding. This hypothesis originated from an earlier study that indicated the approach of the protease flaps upon entrapment 9544



Article

All of these residues are found at or near the active site and the flap regions of the protease. This highlights the importance of the particular regions in SQV binding. In our previous work, we showed that the majority of the single mutants present HB interactions with amino acids in the active site region and at the flaps.50 The only exceptions were mutants V82A, I84V, and L63P, which created HBs with SQV and amino acids only in the active-site region. The double and triple mutants showed similar patterns, interacting with both active-site and flap residues, except G48V-V82A and L10I-L63P, where SQV−flap interactions are absent (Table S2). This implicates that point mutations 63 and 82 are crucial to the reduction of interactions compared to the WT and support their role as destabilizing factors. As suggested in the conformational analysis section, however, the accumulation of mutations reduces the structural flexibility of single mutants L63P, V82A, and A71V. The above observation regarding mutations at positions 63 and 71 is further supported, as complexes L10I-L63P and L10I-L63PA71V showed more HBs than complex L63P and complex L10I-L63P-A71V has many more HBs than does A71V. An enhancement in interactions is also observed in double and triple mutants, which contain G48V or I84V, compared to single mutants G48V and I84V. For comparison, SQV-HIV-1 PR HB interactions are given in Table S2c. Further analysis and comparison of HBs with the single mutants revealed that mutant G48V-V82A presents only two HBs with residues Gly27 and Asp29 in the active-site region with small occupancy. These mutations have been reported as conferring major resistance to SQV.53,88,89 This observation denotes the importance of mutations in positions 48/48′ as well as the detrimental effect of V82A mutation on SQV binding by increasing the level of drug resistance in patients; compared to the other multiple mutants, only strain L10IL63P-A71V-G73S-I84V-L90M shows a similar pattern. Although SQV shows great conformational diversity in double and triple mutants (Figure 3), we do not observe a clear correlation between the number of interactions and conformational changes in the drug. This suggests that mutations of certain residues (e.g., amino acids 48/48′ and 82/82′) are more important than the accumulation of a large number of mutated residues. Following the HB calculations (Figure 5 and Table S2), we assume that drugs designed to retain HB interactions with several residues at or near the active site as well as with residues at the flaps of the protease would lead to effective inhibition. Most significantly, it appears that interactions between the drug and active-site residues Asp25/25′ must be present. Mutants that do not involve SQV-Asp25/25′ HBs (e.g., L10I-L63P, L10I-G48V-V82A, Figure 5) may lead to significantly diminished binding. The importance of the interaction between Asp25/25′ and the drug has been documented in different studies.27,40,45,47,90 Another feature of SQV binding is the presence of watermediated HBs between the enzyme and the drug.50,53 In a previous work, we showed that the WT and several single mutants (L10I, L63P, A71V, and G73S) also presented these interactions.50 Although this was not the case for G48V, V82A, and I84V, all of their respective double mutants formed water bridges with residues at positions 8′ and 29 (Table 4). Of particular importance was the finding that in A71V and G73S, water molecules stabilized SQV in the cavity, but this type of interaction vanishes in all double and triple mutants containing A71V or G73S, except L10I-L63P-A71V.

Table 4. WT and Mutant Complexes of HIV-1 PR with Water-Mediated Interactions between SQV and the Protease complexes WT G48V-V82A L10I-G48V G48V-L90M I84V-L90M L10I-V82A L10I-L63PA71V

water-mediated interactions Asp29 Arg8′ Asp29 Asp29/Arg8′ Asp29 Asp29 Asp29/Gly48′

water

saquinavir

occupancy (%) 39 30 29 20/16 15 13 10/14

III.6. Energy Analysis. In our previous work regarding the impact of single mutations in SQV binding,50 we found that G48V confers medium-level resistance to SQV (ΔΔGMM−PBSA = 1.15 kcal/mol, with respect to the WT). Similarly, this happened in single mutants L63P (ΔΔGMM−PBSA= 1.36 kcal/ mol) and I84V (ΔΔGMM−PBSA = 1.35 kcal/mol). The mutation at positions 48/48′ is considered to be of major importance for SQV resistance since it is located at the flaps of the protease which modulate the entrance of ligands into the binding cavity of the enzyme. All other single mutations presented lower resistance to the drug with ΔΔG values ranging from 0.21 kcal/ mol (for A71V) to 0.72 kcal/mol (for L10I). The resistance due to the G48V mutation is further enhanced when combined with the L10I mutation (ΔΔGMM−PBSA = 1.39 kcal/mol with respect to the WT) or with the L90M mutation (ΔΔGMM−PBSA = 1.36 kcal/mol).27,89 In the double mutant G48V-V82A (ΔΔGMM−PBSA = 2.03 kcal/ mol) we observed an even higher impact in comparison to WT. Therefore, the observation that mutations at positions 82/82′ induce great effects in the SQV conformation (Figures 2 and 3) indicates that the drug rearranges inside the binding cavity. Along with L10I-G48V-V82A and the sextuple mutant, complexes G48V-V82A, L10I-L63P, and A71V-G73S showed the most significant reduction in binding affinity compared to the WT (Tables 5 and 6). This supports the above observation (HB analysis section) regarding the necessity of HB interactions between the drug and Asp25/25′ to maintain effective binding. The lack of these interactions in the sextuple, L10I-G48V-V82A, L10I-L63P, and G48V-V82A complexes may rationalize their diminished binding affinities for SQV. Importantly, the available experimental data (for L10I-G48VV82A and L10I-L63P-A71V-G73S-I84V-L90M) provide additional support for the adequacy of our computational approach (Table 6). The reduced binding energy in complex L10I-V82A may be justified by the nonstabilization of SQV in the cavity, as observed in Figure S2. On the other hand, SQV in A71V-G73S eventually stabilized despite the significant conformational change at ∼11 ns (Figure 2). However, this conformational change of the drug may still affect its binding energy and may explain the reduced binding (ΔΔGMM−PBSA = 2.03 kcal/mol). Also, SQV structures in mutants that acquire a more bent conformation compared to SQV structures in WT (Figure 3) appear to have reduced binding. Table 7 displays the relative values of the inhibition constants (Ki) with respect to the WT for each complex. The calculated values for the inhibition constant are in agreement with values reported in the literature.91,92 Lines 1−9 clearly demonstrate the negative effect of each single mutation on SQV binding (Ki/Kj > 1). Double (lines 10−16) and multiple (lines 17−20) mutations 9545



Article

Table 5. MM−PBSA Results in WT HIV1-PR−SQV and Double Mutant Complexes HIV-1 PR−SQV energy analysis (kcal/mol)a

WT

L10I-G48V

G48V-V82A

G48V-L90M

L10I-L63P

A71V-G73S

I84V-L90M

L10I-V82A

ΔEele ΔEvdw ΔEMM ΔGNP ΔGPB ΔGsolv ΔH(MM+solv) −TΔStot ΔGMM−PBSA Ki (nM)c ΔΔGMM−PBSAb

−41.28 −60.17 −101.45 −8.98 73.23 64.25 −37.20 24.87 −12.33 1.04(0.42)d

−40.27 −56.83 −97.10 −9.06 72.77 63.71 −33.39 22.45 −10.94 10.7 1.39

−43.62 −56.27 −99.89 −9.21 74.54 65.33 −34.56 24.28 −10.30 31.4 2.03

−43.27 −56.93 −100.20 −8.40 73.41 65.01 −35.19 24.22 −10.97 10.2 1.36

−44.55 −55.02 −99.57 −7.63 73.85 66.22 −33.35 23.19 −10.16 39.7 2.17

−42.56 −60.09 −102.65 −8.50 77.26 68.76 −33.89 23.59 −10.30 31.4 2.03

−40.39 −59.68 −100.07 −8.67 74.40 65.73 −34.34 23.00 −11.34 5.5 0.99

−42.48 −58.77 −101.25 −8.27 73.38 65.11 −36.14 25.45 −10.69 16.3 1.64

Calculations were performed on the last 5 ns of each trajectory. bΔΔGMM−PBSA is the difference between each mutant and WT. Error bars for the energy calculations are shown in Table S8. The description of MM−PBSA terms is found in the Supporting Information. cThe conversion of calculated ΔG values to Ki was obtained at 300 K. dExperimental value in parentheses.53 a

Table 6. MM−PBSA Results in HIV1-PR−SQV Multiple Mutant Complexes

Table 7. Inhibition Constants (Ki) and Relative Values (Ki/ Kj) in Different HIV-1 PR Complexesa,b

HIV-1 PR−SQV complex

energy analysis (kcal/mol)a

L10IL63PA71V

G73SI84VL90M

L10IG48VV82A

L10I-L63P-A71VG73S-I84V-L90M

ΔEele ΔEvdw ΔEMM ΔGNP ΔGPB ΔGsolv ΔH(MM+solv) −TΔStot ΔGMM−PBSA Ki (nM)b ΔΔGMM−PBSA

−40.62 −59.27 −99.89 −9.01 70.76 61.75 −38.14 26.38 −11.76 2.71 0.57

−43.74 −57.65 −101.39 −7.68 73.12 65.44 −35.95 25.30 −10.65 17.4 1.68

−42.68 −57.62 −100.30 −6.79 72.45 65.66 −34.64 25.14 −9.50 120(90)c 2.83

−38.22 −59.84 −98.06 −7.12 59.67 66.79 −31.27 21.34 −9.93 58.4(78)c 2.40

1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17

a

Calculations were performed on the last 5 ns of each trajectory. Conversion of calculated ΔG values to Ki was obtained at 300 K. c Experimental values are in parentheses.91 Error bars for the energy calculations are shown in Table S9. b

18

may have relative inhibition constants that vary from improved (Ki/Kj < 1) to weakened (Ki/Kj > 1) binding, depending on which mutant the subsequent mutation accumulation has acted. It may be concluded that the accumulation of mutations on single or double mutant motifs is most likely to produce higherorder mutants (triple and sextuple mutations) with less effective SQV binding. Nevertheless, we observed that in some cases mutation accumulation does not necessarily lead to reduced binding affinity. One such example is triple mutant L10I-L63P-A71V (Table 7). While single and double mutants L10I, L63P, A71V, and L10I-L63P have reduced binding affinities for SQV (Ki = 7.4, 21, 3.2, and 39.7 nM, respectively), the combination of the mutations in the triple mutant shows enhanced binding (Ki = 2.71 nM). Another two cases include I84V-L90M, which presents more efficient SQV binding (Ki = 5.5 nM) than the individual single mutant strains I84V (Ki = 20.6 nM) and L90M (Ki = 18.2 nM), and G48V-L90M. Importantly, these observations correlate with the above findings regarding structure stabilization in these complexes (L10I-L63P-A71V and I84V-L90M, Figures 2, 4, S1, and S7 and Table 1) and clearly suggest that different mutations do not always have an additive effect in drug binding.

19 20

WT L63P V82A G48V A71V G73S I84V L10I L90M L10I-G48V G48V-V82A L10I-L63P A71V-G73S G48V-L90M I84V-L90M L10I-V82A L10I-L63P-A 71V G73S-I84VL90M L10I-G48VV82A L10I-L63PA71VG73S-I84VL90M

Ki (nM)

Ki/Kj

1.04 21 5.1 14.5 3.2 4.6 20.6 7.4 18.2 10.7 31.4 39.7 31.4 10.2 5.5 16.3 2.71

1 20.19 4.90 13.94 3.077 4.42 19.81 7.12 17.50 relative to G48V: 0.74/relative to L10I: 1.45 G48V: 2.17/V82A: 6.16 L10I: 5.36/L63P: 1.89 A71V: 9.81/G73S: 6.83 G48V: 0.70/L90M: 0.56 I84V: 0.27/L90M: 0.30 L10I: 2.20/V82A: 3.20 10: 0.37/63: 0.13/71: 0.85/10−63: 0.07

17.4

73: 3.78/84: 0.84/84−90: 3.16

120

10: 16.22/48: 8.3/82: 23.5/10−48: 11.21/48− 82: 3.82/10−82: 7.4 10: 7.89/63: 2.78/71: 18.25/73: 12.70/84: 2.8/ 90: 3.2/10−63: 1.5/71−73: 1.86/84−90: 10.62/10−63−71: 21.5/73−84− 90: 3.36

58.4

a

Kj is: Ki,WT (= 1.04) for calculations 1−9; Ki of individual single mutations (2−9) for calculations 10−16; Ki of single and double mutations (2−16) for calculations 17−19; Ki of single, double, and triple mutations (2−19) for calculation 20. bIndex: If ratio > 1 → weaker binding; if ratio < 1 → stronger binding.

The results of Tables 5 and 6 indicate that despite the fact that SQV binding in all mutant (and WT) complexes is mostly enthalpically driven (ΔH ≈ −35 kcal/mol), the entropic component is still very important, as also noticed in previous studies.27,29,50 In Figure S8, we present the enthalpy and entropy plots vs time for all complexes. Furthermore, it was observed that the total enthalpy loss upon mutations was mostly due to the loss in enthalpy around the region of the active site (Table S7). The decomposition of the binding energy into individual parts showed that contributions are fairly 9546



Article

Table 8. ΔΔGvdW (= ΔGvdWMut − ΔGvdWWT) and ΔΔGele (= ΔGeleMut − ΔGeleWT) for the Entire Protease, for Active-Site Residues (21/21′−31/31′), and for Flap Residues (44/44′−55/55′) in HIV-1 PR Complexesa

a

complex

ΔΔGvdW

ΔΔGele

ΣΔΔGvdW, active site

ΣΔΔGvdW, flaps

ΣΔΔGele, active site

ΣΔΔGele, flaps

G48V-V82A L10I-G48V G48V-L90M I84V-L90M L10I-V82A L10I-L63P A71V-G73S L10I-L63P-A71V L10I-G48V-V82A G73S-I84V-L90M L10I-L63P-A71V-G73S-I84V-L90M

−0.26 0.75 −0.19 2.15 3.84 2.43 −1.02 −0.30 −2.14 −0.99 −3.57

−0.45 1.01 1.85 1.27 −1.67 1.65 0.33 0.37 4.13 1.14 6.76

3.32 3.30 3.17 4.86 4.77 2.39 2.09 2.80 2.25 3.58 3.25

2.33 1.50 2.66 2.22 2.66 3.07 2.80 1.85 2.75 2.15 2.62

5.48 2.16 1.66 −0.01 0.10 −0.43 −0.83 0.06 −0.22 0.56 1.19

1.47 2.65 2.93 −1.66 −1.04 −0.09 1.51 −0.35 2.95 0.58 2.39

Per residue ΔΔG (the active site and flaps) for van der Waals and electrostatic energies is shown in Tables S3−S6.

Figure 6. (a, b) The structure of Gly49-Ile50-Gly51−water−Gly49′-Ile50′-Gly51′ (GWG). The geometries, which resulted from MD simulations, are associated with complexes (a) WT−SQV and (b) L10I-L63P-A71V(SQV). N, C, O and H are shown in blue, gray, red and white, respectively. (c, d) Superposition of GWG geometries: (c) The blue sticks depict the optimized (B3LYP/6-31G*) geometry, while the red sticks represent the structure derived by the MD run for WT. (d) MD structures of WT−SQV (red sticks) and L10I-L63P-A71V(SQV) (cyan sticks).

interactions between WT and each mutant may vary from favorable to unfavorable, an interesting feature was observed upon further inspection (Table 8): in all mutant complexes, a significant loss in van der Waals interactions was noticed at the active site and in the flap regions. This occurred consistently regardless of any total loss or gain in van der Waals interactions upon mutation, thus implicating the active site and flap regions

similar among complexes. The van der Waals interactions have the highest positive impact on the binding energy, followed by the nonpolar contribution to solvation (Tables 5 and 6). Similar observations (mostly regarding high van der Waals contributions) have been reported for different HIV-1 PR complexes, suggesting that this is a common trend for the protease.29,45,53,93 Even though changes in van der Waals 9547



Article

Table 9. Interaction Energy Decomposition Analysis of Selected Fragmentsa GWGb

Eel

Eex

Erep

Epl

Edisp

ΔE

dimers of GWG 1/2c 1/3c 2/3c 1/2d 1/3d 2/3d GWGe GWGg GWGh SQV-Arg8i

−4.90 −8.12 −3.21 −15.18 −8.23 −7.30 −15.96 (−12.46)f −30.17 (−25.21)f −20.48 −9.67

−15.20 −5.64 −1.97 −19.67 −8.27 −6.91 −22.78 (−13.11) −34.76 (−20.59) −31.26 −6.19

25.93 9.74 3.30 32.77 14.25 11.90 38.86 (24.72) 58.64 (38.02) 54.45 10.78

−3.53 −1.52 −0.55 −4.31 −1.83 −1.55 −4.67 (−9.49) −6.78 (−12.94) −6.36 −3.0

−7.03 −0.94 −0.71 −5.77 −1.37 −1.26 −8.59 (−14.08) −8.22 (−15.40) −5.30 −2.18

−4.73 −6.48 −3.14 −12.16 −5.46 −5.14 −13.13 (−24.42) −21.28 (−36.13) −8.95 −10.26

a The MP2/6-31G* method was used. Energy values (kcal/mol) are BSSE-corrected. bGWG: Gly49-Ile50-Gly51−H2O−Gly49′-Ile50′-Gly51′. cThe geometry of GWG has been taken from the MD simulation (Figure 6a). ΔE analysis of the three dimers is presented. dThe geometry of GWG has been optimized with B3LYP/6-31G*. eThe geometry of GWG was obtained from the MD simulation of WT−SQV. fNon-BSSE-corrected value. g The geometry of GWG has been optimized with B3LYP/6-31G*. hThe geometry of the GWG fragment is the average structure of those derived from the MD simulation of the L10I-L63P-A71V(SQV) complex. iThe geometry of Arg8−SQV (WT) was obtained from the MD simulation.

Substantial differences were also observed in most ΔE components, associated with the two geometries (e.g., Ees: −30.17/−15.96 kcal/mol). Similarly, when employing the KM energy decomposition scheme, the geometry also had a great effect on Eel and Exr (Table S12). The optimized geometry had a lower ΔE value (compared to that resulting from the MD simulation). This is due to the gas-phase-optimized geometry, where monomer 2 is closer to the water molecule (monomer 3), and this enhances the HB effect (N3−H3−O2, Figure 6). The results denote the impact of H2O on ΔE (−21.28/−12.16 kcal/mol). The sum of the interaction energies (ΔE) of dimers 1/3 and 2/3 is −9.62 kcal/mol at the geometry computed by the MD simulation. This sum slightly differs (−10.6 kcal/mol) when the geometry is optimized (B3LYP/6-31G*). The rather small change in the Epl term implies small changes in the shape of the orbitals of the dimers. III.7.3. L10I-L63P-A71V(SQV). The employed geometry is the average structure of those derived from the MD run of the L10I-L63P-A71V(SQV) complex. In comparison to ΔE of WT−SQV (−13.23 kcal/mol), a remarkable decrease (32%) was observed for this mutant (−8.95 kcal/mol). The superposition of the two associated MD geometries is shown in Figure 6d. In L10I-L63P-A71V(SQV), the H2O molecule is closer to the residues compared to the WT−SQV complex. This leads to an increase in the negative terms (Eel + Eex + Epl + Edisp) for L10I-L63P-A71V(SQV), which favor binding. The corresponding sums are −63.4 and −52.0 kcal/mol for the mutant and WT−SQV, respectively. However, the increase in the repulsion term, Erep, in L10I-L63P-A71V(SQV) is larger (54.5 kcal/mol) than in the WT−SQV complex (38.9 kcal/ mol). It was therefore suggested that WT−SQV formation is favored (ΔE = −8.95 kcal/mol for mutant/−13.1 kcal/mol for WT). III.7.4. Arg8−SQV. The average structure of those derived from the MD simulation of the WT−SQV complex has been used for the calculation. A remarkably high ΔE (−10.26 kcal/ mol) has been associated with favorable interactions between hydrogen atoms of the two −NH2 groups of arginine and two oxygen atoms of SQV (Figure S9).

as significant van der Waals stabilizing factors (Table 8). The electrostatic energies at the active site and the flap areas, however, deviate from this pattern by displaying both favorable and unfavorable contributions in the set of mutants (Tables S5 and S6). Further MM−PBSA calculations for the second half (20−40 ns) of the 40 ns simulations mentioned above (L10I-V82A and A71V-G73S, Figure S3) showed that enthalpy differences between the same complexes are ∼3−5 kcal/mol. This indicates the reproducibility of our results, since these values are within the error (standard deviation) range shown in Tables S8 and S9. The comparison (including all energy terms) is shown in Table S10. III.7. Ab Initio Analysis of Intermolecular Interactions. Employing the methods of Su and Li74 and Kitaura and Morokuma (KM),75 we have computed the interaction energy, ΔE, of some molecular fragments at the MP2/6-31G* level. These computations demonstrate the effect of geometry on ΔE, as well as the components, which contribute to ΔE. III.7.1. Analysis of the Interactions in Fragment Gly49Ile50-Gly51−water−Gly49′-Ile50′-Gly51′. Energy decomposition analysis (EDA) of the above system (briefly denoted as GWG, Figure 6a) has been performed with the method of Su and Li.74 GWG is a flap fragment of HIV-1 PR. We used two conformations of GWG: one has been taken from the WT− SQV complex, and the other, from the L10I-L63P-A71V(SQV) complex (Figure 6b). Each glycine (49, 49′, 51, 51′) has been modeled as a terminal residue, with −COOH and −NH2 end groups, to produce a closed-shell system. III.7.2. WT−SQV. Two different geometries of GWG were used in order to evaluate the interaction energy and its components. The first one is the average structure of those derived from the MD simulation of the WT−SQV complex (Figure 6a). The second one was optimized, in the gas phase, using the B3LYP/6-31G* theory (Figure 6c). We have also performed EDA computations for the three dimers (1/2, 1/3 or 2/3, see Figure 6a) of the above system. The results (Table 9) have been corrected for the basis set superposition error,81 the effect of which is significant. The employed geometries produced significantly different ΔE values (−13.13 kcal/mol for MD structure, −21.28 kcal/mol for optimized structure), although these geometries do not differ greatly (Figure 6c, Table S11). This denotes the major effect of geometry on ΔE.

IV. CONCLUSIONS A systematic approach to drug-resistance-associated mutations G48V-V82A, L10I-G48V, G48V-L90M, I84V-L90M, L10I9548



Article

V82A, L10I-L63P, A71V-G73S, L10I-L63P-A71V, L10I-G48VV82A, G73S-I84V-L90M, and L10I-L63P-A71V-G73S-I84VL90M in HIV-1 PR has been attempted by means of molecular dynamics (MD), molecular mechanics Poisson−Boltzmann surface area (MM−PBSA), and ab initio computational methodologies. The conformational features and binding modes of the first marketed, potent protease inhibitor, saquinavir (SQV), have been compared between the wild type (WT) and the 11 mutated forms of the protease to clarify the mechanism of resistance at the molecular level. A comparison of previous extensive studies from our group on the HIV-1 PR binding modes and resistance50,54−56 with available experimental results has demonstrated the ability of the above computational methodologies to provide very reliable information. According to conformational, hydrogen bond, and energy analyses of WT and mutated HIV-1 PR complexes, the SQV resistance mechanism in HIV-1 PR includes the following elements: (1) The flexibility of SQV may affect the binding significantly. Previous work from our group suggested that conformational changes of the drug may not reduce the binding, provided that its structure eventually stabilizes.50 Here, saquinavir’s conformations in all double and multiple mutants appear quite stable, but it was also observed that bent SQV structures in L10I-G48V-V82A, L10I-L63P-A71V-G73S-I84VL90M, L10I-L63P, and G48V-V82A deviate from their structures in WT and thus contribute to the reduction of the binding affinity for the protease. (2) The accumulation of mutations does not necessarily produce incremental effects. Smaller conformational changes in the structures of HIV-1 PR, SQV, and the flaps were observed for double (I84V-L90M, G48V-L90M) and triple (L10I-L63P-A71V) mutants than for individual single mutations (L10I, G48V, L63P, I84V, L90M). Also, L10I-L63P and L10I-L63P-A71V complexes presented an enhanced (compared to single mutants) HB network that stabilized SQV in the binding cavity. Interestingly, the addition of A71V to the double mutant facilitated a water-mediated HB interaction between SQV and the active site of the protease, thus resulting in improved binding. (3) Drugs designed to retain HB interactions with several residues at or near the active site (e.g., Arg8, Gly27/27′, Asp25/25′, Asp29/29′, Asp30/30′) as well as with residues at the flaps (e.g., Ile50/50′, Gly(Val)48/ 48′) of the protease lead to effective inhibition. Importantly, HBs between SQV and the active site of HIV-1 PR (especially Asp25/25′) are necessary for effective binding. Thus, mutants where SQV−Asp25/25′ interactions are prevalent (L10I-L63PA71V and G48V-L90M) displayed the closest binding affinities to the affinity for the WT. In mutant complexes L10I-G48VV82A, L10I-L63P, G48V-V82A, and the sextuple, a structural rearrangement of the binding cavity induced the displacement of SQV in a way that interactions with Asp25/25′ are not favored (Figure 7). This resulted in a significant binding energy loss (ΔΔG > 2 kcal/mol). Of major importance was the observation that the loss of enthalpy in the active-site region accounted almost entirely for the total enthalpy loss in all mutants. (4) A closed conformation of the flaps does not necessarily imply efficient binding. In agreement with our previous studies (including single mutants and several ligands),50,56 flap closing occurs in the double, triple, and sextuple complexes, regardless the binding strength. (5) Complexes which contain the L10I, L63P, I84V, or L90M single mutations will stabilize their structure and present enhanced HB interactions and binding strength upon mutation

Figure 7. Structural rearrangement of the binding cavity upon L10IG48V-V82A mutation. SQV in WT forms hydrogen bonds with Asp25/25′ (purple), while these interactions vanish in the mutated form (green).

accumulation. Complexes I84V-L90M and L10I-L63P-A71V yielded the highest binding affinities among mutants. This denotes the role of other mutations, such as A71V, as possible stabilization factors. (6) In all mutant complexes, a loss in van der Waals interactions was observed at the active site and in the flap regions. This occurred consistently regardless of any total loss or gain in van der Waals interactions upon mutation. (7) Small changes in geometry may have a great effect on the binding energy. Ab initio computations of the interaction energy (ΔE) of selected SQV−HIV-1 PR fragments showed that geometries, which do not differ greatly, yielded significantly different ΔE values (ΔΔE ≈ 8 kcal/mol). The major contributions to the interaction energy are due to electrostatic and exchange components. Finally, it has been observed that H2O has a significant effect on ΔE.

■

ASSOCIATED CONTENT

S Supporting Information *

Detailed description of the methods employed in this article. Detailed RMSD graphs for SQV and the protease for all mutated complexes. Representative conformations of SQV in the different single mutants. Depictions of ΔH and −TΔS over time for all complexes. Data on atomic fluctuations, hydrogen bonds, ΔΔE values for every mutated complex, selected bond lengths, and interaction energy decomposition analysis. This material is available free of charge via the Internet at http:// pubs.acs.org.

■

AUTHOR INFORMATION

Corresponding Authors

*(H.T.) Tel: +30-210-727-3894. E-mail: haralambostz@gmail. com. *(G.L.) Tel: +30-210-727-3894. E-mail: [email protected]. *(M.G.P.) Tel: +30-210-727-3892. E-mail: [email protected]. Notes

The authors declare no competing financial interest.

■

ACKNOWLEDGMENTS This research has been cofinanced by the European Union (European Social Fund − ESF) and Greek national funds through the Operational Program “Education and Lifelong Learning” of the National Strategic Reference Framework (NSRF) − Research Funding Program: Heracleitus II− 9549



Article

(16) Molla, A.; Korneyeva, M.; Gao, Q.; Vasavanonda, S.; Schipper, P. J.; Mo, H. M.; Markowitz, M.; Chernyavskiy, T.; Niu, P.; Lyons, N.; et al. Ordered Accumulation of Mutations in Hiv Protease Confers Resistance to Ritonavir. Nat. Med. 1996, 2, 760−766. (17) Condra, J. H.; Schleif, W. A.; Blahy, O. M.; Gabryelski, L. J.; Graham, D. J.; Quintero, J. C.; Rhodes, A.; Robbins, H. L.; Roth, E.; Shivaprakash, M.; et al. In Vivo Emergence of Hiv-1 Variants Resistant to Multiple Protease Inhibitors. Nature 1995, 374, 569−571. (18) Dunn, B. M. Anatomy and Pathology of Hiv-1 Peptidase. Essays Biochem. 2002, 38, 113−127. (19) Baldwin, E. T.; Bhat, T. N.; Liu, B.; Pattabiraman, N.; Erickson, J. W. Structural Basis of Drug Resistance for the V82a Mutant of Hiv-1 Proteinase. Nat. Struct. Biol. 1995, 2, 244−249. (20) Erickson, J. W.; Burt, S. K. Structural Mechanisms of Hiv Drug Resistance. Annu. Rev. Pharmacol. Toxicol. 1996, 36, 545−571. (21) Otto, M. J.; Garber, S.; Winslow, D. L.; Reid, C. D.; Aldrich, P.; Jadhav, P. K.; Patterson, C. E.; Hodge, C. N.; Cheng, Y. S. In Vitro Isolation and Identification of Human Immunodeficiency Virus (Hiv) Variants with Reduced Sensitivity to C-2 Symmetrical Inhibitors of Hiv Type 1 Protease. Proc. Natl. Acad. Sci. U.S.A. 1993, 90, 7543− 7547. (22) el-Farrash, M. A.; Kuroda, M. J.; Kitazaki, T.; Masuda, T.; Kato, K.; Hatanaka, M.; Harada, S. Generation and Characterization of a Human Immunodeficiency Virus Type 1 (Hiv-1) Mutant Resistant to an Hiv-1 Protease Inhibitor. J. Virol. 1994, 68, 233−239. (23) Genoni, A.; Morra, G.; Merz, K. M., Jr.; Colombo, G. Computational Study of the Resistance Shown by the Subtype B/Hiv1 Protease to Currently Known Inhibitors. Biochemistry 2010, 49, 4283−4295. (24) Dandache, S.; Sevigny, G.; Yelle, J.; Stranix, B. R.; Parkin, N.; Schapiro, J. M.; Wainberg, M. A.; Wu, J. J. In Vitro Antiviral Activity and Cross-Resistance Profile of Pl-100, a Novel Protease Inhibitor of Human Immunodeficiency Virus Type 1. Antimicrob. Agents Chemother. 2007, 51, 4036−4043. (25) King, J. R.; Wynn, H.; Brundage, R.; Acosta, E. P. Pharmacokinetic Enhancement of Protease Inhibitor Therapy. Clin. Pharmacokinet. 2004, 43, 291−310. (26) Bergersen, B. M.; Sandvik, L.; Dunlop, O.; Birkeland, K.; Bruun, J. N. Prevalence of Hypertension in Hiv-Positive Patients on Highly Active Retroviral Therapy (Haart) Compared with Haart-Naive and Hiv-Negative Controls: Results from a Norwegian Study of 721 Patients. Eur. J. Clin. Microbiol. Infect. Dis. 2003, 22, 731−763. (27) Ode, H.; Neya, S.; Hata, M.; Sugiura, W.; Hoshino, T. Computational Simulations of Hiv-1 Proteases–Multi-Drug Resistance Due to Nonactive Site Mutation L90m. J. Am. Chem. Soc. 2006, 128, 7887−7895. (28) Toth, G.; Borics, A. Flap Opening Mechanism of Hiv-1 Protease. J. Mol. Graphics Modell 2006, 24, 465−474. (29) Hou, T.; McLaughlin, W. A.; Wang, W. Evaluating the Potency of Hiv-1 Protease Drugs to Combat Resistance. Proteins 2008, 71, 1163−1174. (30) Kovalevsky, A. Y.; Liu, F.; Leshchenko, S.; Ghosh, A. K.; Louis, J. M.; Harrison, R. W.; Weber, I. T. Ultra-High Resolution Crystal Structure of Hiv-1 Protease Mutant Reveals Two Binding Sites for Clinical Inhibitor Tmc114. J. Mol. Biol. 2006, 363, 161−173. (31) Kovalevsky, A. Y.; Tie, Y.; Liu, F.; Boross, P. I.; Wang, Y. F.; Leshchenko, S.; Ghosh, A. K.; Harrison, R. W.; Weber, I. T. Effectiveness of Nonpeptide Clinical Inhibitor Tmc-114 on Hiv-1 Protease with Highly Drug Resistant Mutations D30n, I50v, and L90m. J. Med. Chem. 2006, 49, 1379−1387. (32) Tie, Y.; Kovalevsky, A. Y.; Boross, P.; Wang, Y. F.; Ghosh, A. K.; Tozser, J.; Harrison, R. W.; Weber, I. T. Atomic Resolution Crystal Structures of Hiv-1 Protease and Mutants V82a and I84v with Saquinavir. Proteins 2007, 67, 232−242. (33) Lee, T.; Laco, G. S.; Torbett, B. E.; Fox, H. S.; Lerner, D. L.; Elder, J. H.; Wong, C. H. Analysis of the S3 and S3′ Subsite Specificities of Feline Immunodeficiency Virus (Fiv) Protease: Development of a Broad-Based Protease Inhibitor Efficacious against

Investing in Knowledge Society through the European Social Fund. Funding provided by the European Commission for FP7-REGPOT-2009-1 project “ARCADE” (grant agreement no. 245866), for FP7-PEOPLE-2011-IRSES project “NanoBRIDGES” (grant agreement no. PIRSES-GA-2011-295128), and for project “NanoPUZZLES” (grant agreement no. NMPSL-2012-309837) is also acknowledged.

■

REFERENCES

(1) Croteau, G.; Doyon, L.; Thibeault, D.; McKercher, G.; Pilote, L.; Lamarre, D. Impaired Fitness of Human Immunodeficiency Virus Type 1 Variants with High-Level Resistance to Protease Inhibitors. J. Virol. 1997, 71, 1089−1096. (2) Wlodawer, A.; Erickson, J. W. Structure-Based Inhibitors of Hiv-1 Protease. Annu. Rev. Biochem. 1993, 62, 543−585. (3) Carpenter, C. C.; Fischl, M. A.; Hammer, S. M.; Hirsch, M. S.; Jacobsen, D. M.; Katzenstein, D. A.; Montaner, J. S.; Richman, D. D.; Saag, M. S.; Schooley, R. T.; et al. Antiretroviral Therapy for Hiv Infection in 1997. Updated Recommendations of the International Aids Society-USA Panel. J. Am. Med. Assoc. 1997, 277, 1962−1969. (4) FDA http://www.fda.gov/ForConsumers/ByAudience/ ForWomen/ucm118597.htm, accessed March 12, 2013. (5) Mastrolorenzo, A.; Rusconi, S.; Scozzafava, A.; Barbaro, G.; Supuran, C. T. Inhibitors of Hiv-1 Protease: Current State of the Art 10 Years after Their Introduction. From Antiretroviral Drugs to Antifungal, Antibacterial and Antitumor Agents Based on Aspartic Protease Inhibitors. Curr. Med. Chem. 2007, 14, 2734−2748. (6) Hammer, S. M.; Squires, K. E.; Hughes, M. D.; Grimes, J. M.; Demeter, L. M.; Currier, J. S.; Eron, J. J., Jr.; Feinberg, J. E.; Balfour, H. H., Jr.; Deyton, L. R.; et al. A Controlled Trial of Two Nucleoside Analogues Plus Indinavir in Persons with Human Immunodeficiency Virus Infection and Cd4 Cell Counts of 200 Per Cubic Millimeter or Less. Aids Clinical Trials Group 320 Study Team. N. Engl. J. Med. 1997, 337, 725−733. (7) Coffin, J. M. Hiv Population Dynamics in Vivo: Implications for Genetic Variation, Pathogenesis, and Therapy. Science 1995, 267, 483−489. (8) Imamichi, T. Anti-Hiv Drugs and Resistance: Inhibitors of Hiv-1 Reverse Transcriptase and Protease. Front. Med. Chem. 2006, 3, 23− 44. (9) Nicholson, L. K.; Yamazaki, T.; Torchia, D. A.; Grzesiek, S.; Bax, A.; Stahl, S. J.; Kaufman, J. D.; Wingfield, P. T.; Lam, P. Y.; Jadhav, P. K.; et al. Flexibility and Function in Hiv-1 Protease. Nat. Struct. Biol. 1995, 2, 274−280. (10) Freedberg, D. I.; Ishima, R.; Jacob, J.; Wang, Y. X.; Kustanovich, I.; Louis, J. M.; Torchia, D. A. Rapid Structural Fluctuations of the Free Hiv Protease Flaps in Solution: Relationship to Crystal Structures and Comparison with Predictions of Dynamics Calculations. Protein Sci. 2002, 11, 221−232. (11) Lapatto, R.; Blundell, T.; Hemmings, A.; Overington, J.; Wilderspin, A.; Wood, S.; Merson, J. R.; Whittle, P. J.; Danley, D. E.; Geoghegan, K. F.; et al. X-Ray Analysis of Hiv-1 Proteinase at 2.7 a Resolution Confirms Structural Homology among Retroviral Enzymes. Nature 1989, 342, 299−302. (12) Hornak, V.; Okur, A.; Rizzo, R. C.; Simmerling, C. Hiv-1 Protease Flaps Spontaneously Close to the Correct Structure in Simulations Following Manual Placement of an Inhibitor into the Open State. J. Am. Chem. Soc. 2006, 128, 2812−2813. (13) Hornak, V.; Okur, A.; Rizzo, R. C.; Simmerling, C. Hiv-1 Protease Flaps Spontaneously Open and Reclose in Molecular Dynamics Simulations. Proc. Natl. Acad. Sci. U.S.A. 2006, 103, 915− 920. (14) Vondrasek, J.; Wlodawer, A. Hivdb: A Database of the Structures of Human Immunodeficiency Virus Protease. Proteins 2002, 49, 429−431. (15) Flexner, C. Hiv-Protease Inhibitors. N. Engl. J. Med. 1998, 338, 1281−1292. 9550



Article

Fiv, Siv, and Hiv in Vitro and Ex Vivo. Proc. Natl. Acad. Sci. U.S.A. 1998, 95, 939−944. (34) Mo, H.; Lu, L.; Dekhtyar, T.; Stewart, K. D.; Sun, E.; Kempf, D. J.; Molla, A. Characterization of Resistant Hiv Variants Generated by in Vitro Passage with Lopinavir/Ritonavir. Antiviral Res. 2003, 59, 173−180. (35) Johnson, V. A.; Brun-Vezinet, F.; Clotet, B.; Gunthard, H. F.; Kuritzkes, D. R.; Pillay, D.; Schapiro, J. M.; Richman, D. D. Update of the Drug Resistance Mutations in Hiv-1: December 2009. Top HIV Med. 2009, 17, 138−145. (36) Johnson, J. A.; Li, J. F.; Wei, X.; Lipscomb, J.; Irlbeck, D.; Craig, C.; Smith, A.; Bennett, D. E.; Monsour, M.; Sandstrom, P.; et al. Minority Hiv-1 Drug Resistance Mutations Are Present in Antiretroviral Treatment-Naive Populations and Associate with Reduced Treatment Efficacy. PLoS Med. 2008, 5, e158. (37) Brik, A.; Wong, C. H. Hiv-1 Protease: Mechanism and Drug Discovery. Org. Biomol Chem. 2003, 1, 5−14. (38) Coates, L.; Tuan, H. F.; Tomanicek, S.; Kovalevsky, A.; Mustyakimov, M.; Erskine, P.; Cooper, J. The Catalytic Mechanism of an Aspartic Proteinase Explored with Neutron and X-Ray Diffraction. J. Am. Chem. Soc. 2008, 130, 7235−7237. (39) Hong, L.; Zhang, X. C.; Hartsuck, J. A.; Tang, J. Crystal Structure of an in Vivo Hiv-1 Protease Mutant in Complex with Saquinavir: Insights into the Mechanisms of Drug Resistance. Protein Sci. 2000, 9, 1898−1904. (40) Hou, T.; Yu, R. Molecular Dynamics and Free Energy Studies on the Wild-Type and Double Mutant Hiv-1 Protease Complexed with Amprenavir and Two Amprenavir-Related Inhibitors: Mechanism for Binding and Drug Resistance. J. Med. Chem. 2007, 50, 1177−1188. (41) Nijhuis, M.; van Maarseveen, N. M.; Lastere, S.; Schipper, P.; Coakley, E.; Glass, B.; Rovenska, M.; de Jong, D.; Chappey, C.; Goedegebuure, I. W.; et al. A Novel Substrate-Based Hiv-1 Protease Inhibitor Drug Resistance Mechanism. PLoS Med. 2007, 4, e36. (42) Saen-oon, S.; Aruksakunwong, O.; Wittayanarakul, K.; Sompornpisut, P.; Hannongbua, S. Insight into Analysis of Interactions of Saquinavir with Hiv-1 Protease in Comparison between the WildType and G48v and G48v/L90m Mutants Based on Qm and Qm/Mm Calculations. J. Mol. Graphics Modell. 2007, 26, 720−727. (43) Mammano, F.; Trouplin, V.; Zennou, V.; Clavel, F. Retracing the Evolutionary Pathways of Human Immunodeficiency Virus Type 1 Resistance to Protease Inhibitors: Virus Fitness in the Absence and in the Presence of Drug. J. Virol. 2000, 74, 8524−8531. (44) Lefebvre, E.; Schiffer, C. A. Resilience to Resistance of Hiv-1 Protease Inhibitors: Profile of Darunavir. AIDS Rev. 2008, 10, 131− 142. (45) Alcaro, S.; Artese, A.; Ceccherini-Silberstein, F.; Ortuso, F.; Perno, C. F.; Sing, T.; Svicher, V. Molecular Dynamics and Free Energy Studies on the Wild-Type and Mutated Hiv-1 Protease Complexed with Four Approved Drugs: Mechanism of Binding and Drug Resistance. J. Chem. Inf. Model. 2009, 49, 1751−1761. (46) Prabu-Jeyabalan, M.; Nalivaika, E. A.; King, N. M.; Schiffer, C. A. Viability of a Drug-Resistant Human Immunodeficiency Virus Type 1 Protease Variant: Structural Insights for Better Antiviral Therapy. J. Virol 2003, 77, 1306−1315. (47) Surleraux, D. L.; Tahri, A.; Verschueren, W. G.; Pille, G. M.; de Kock, H. A.; Jonckers, T. H.; Peeters, A.; De Meyer, S.; Azijn, H.; Pauwels, R.; et al. Discovery and Selection of Tmc114, a Next Generation Hiv-1 Protease Inhibitor. J. Med. Chem. 2005, 48, 1813− 1822. (48) Tie, Y.; Boross, P. I.; Wang, Y. F.; Gaddis, L.; Hussain, A. K.; Leshchenko, S.; Ghosh, A. K.; Louis, J. M.; Harrison, R. W.; Weber, I. T. High Resolution Crystal Structures of Hiv-1 Protease with a Potent Non-Peptide Inhibitor (Uic-94017) Active against Multi-DrugResistant Clinical Strains. J. Mol. Biol. 2004, 338, 341−352. (49) Kovalevsky, A. Y.; Ghosh, A. K.; Weber, I. T. Solution Kinetics Measurements Suggest Hiv-1 Protease Has Two Binding Sites for Darunavir and Amprenavir. J. Med. Chem. 2008, 51, 6599−6603. (50) Tzoupis, H.; Leonis, G.; Mavromoustakos, T.; Papadopoulos, M. G. A Comparative Molecular Dynamics, Mm-Pbsa and

Thermodynamic Integration Study of Saquinavir Complexes with Wild-Type Hiv-1 Pr and L10i, G48v, L63p, A71v, G73s, V82a and I84v Single Mutants. J. Chem. Theory Comput. 2013, 9, 1754−1764. (51) Ghosh, A. K.; Chapsal, B. D.; Weber, I. T.; Mitsuya, H. Design of Hiv Protease Inhibitors Targeting Protein Backbone: An Effective Strategy for Combating Drug Resistance. Acc. Chem. Res. 2008, 41, 78−86. (52) Ghosh, A. K.; Chapsal, B. D.; Steffey, M.; Agniswamy, J.; Wang, Y. F.; Amano, M.; Weber, I. T.; Mitsuya, H. Substituent Effects on P2Cyclopentyltetrahydrofuranyl Urethanes: Design, Synthesis, and XRay Studies of Potent Hiv-1 Protease Inhibitors. Bioorg. Med. Chem. Lett. 2012, 22, 2308−2311. (53) Liu, F.; Kovalevsky, A. Y.; Tie, Y.; Ghosh, A. K.; Harrison, R. W.; Weber, I. T. Effect of Flap Mutations on Structure of Hiv-1 Protease and Inhibition by Saquinavir and Darunavir. J. Mol. Biol. 2008, 381, 102−115. (54) Tzoupis, H.; Leonis, G.; Megariotis, G.; Supuran, C. T.; Mavromoustakos, T.; Papadopoulos, M. G. Dual Inhibitors for Aspartic Proteases Hiv-1 Pr and Renin: Advancements in AidsHypertension-Diabetes Linkage Via Molecular Dynamics, Inhibition Assays, and Binding Free Energy Calculations. J. Med. Chem. 2012, 55, 5784−5796. (55) Leonis, G.; Steinbrecher, T.; Papadopoulos, M. G. A Contribution to the Drug Resistance Mechanism of Darunavir, Amprenavir, Indinavir, and Saquinavir Complexes with Hiv-1 Protease Due to Flap Mutation I50V: A Systematic MM−PBSA and Thermodynamic Integration Study. J. Chem. Inf. Model. 2013, 53, 2141−2153. (56) Leonis, G.; Czyznikowska, Z.; Megariotis, G.; Reis, H.; Papadopoulos, M. G. Computational Studies of Darunavir into Hiv1 Protease and Dmpc Bilayer: Necessary Conditions for Effective Binding and the Role of the Flaps. J. Chem. Inf. Model. 2012, 52, 1542−1558. (57) Politi, A.; Leonis, G.; Tzoupis, H.; Ntountaniotis, D.; Papadopoulos, M. G.; Grdadolnik, S. G.; Mavromoustakos, T. Conformational Properties and Energetic Analysis of Aliskiren in Solution and Receptor Site. Mol. Inform. 2011, 30, 973−985. (58) Tzoupis, H.; Leonis, G.; Durdagi, S.; Mouchlis, V.; Mavromoustakos, T.; Papadopoulos, M. G. Binding of Novel Fullerene Inhibitors to Hiv-1 Protease: Insight through Molecular Dynamics and Molecular Mechanics Poisson-Boltzmann Surface Area Calculations. J. Comput.-Aided Mol. Des 2011, 25, 959−976. (59) Case, D. A.; Darden, T. A.; Cheatham, I. T. E.; Simmerling, C.; Wang, J.; Duke, R. E.; Luo, R.; Walker, R. C.; Zhang, W.; Merz, K. M. et al. AMBER 11; University of California: San Francisco, CA, 2010. (60) Hornak, V.; Abel, R.; Okur, A.; Strockbine, B.; Roitberg, A.; Simmerling, C. Comparison of Multiple Amber Force Fields and Development of Improved Protein Backbone Parameters. Proteins 2006, 65, 712−725. (61) Pettersen, E. F.; Goddard, T. D.; Huang, C. C.; Couch, G. S.; Greenblatt, D. M.; Meng, E. C.; Ferrin, T. E. Ucsf Chimera–a Visualization System for Exploratory Research and Analysis. J. Comput. Chem. 2004, 25, 1605−1612. (62) Wang, J.; Wolf, R. M.; Caldwell, J. W.; Kollman, P. A.; Case, D. A. Development and Testing of a General Amber Force Field. J. Comput. Chem. 2004, 25, 1157−1174. (63) Jakalian, A.; Jack, D. B.; Bayly, C. I. Fast, Efficient Generation of High-Quality Atomic Charges. Am1-Bcc Model: Ii. Parameterization and Validation. J. Comput. Chem. 2002, 23, 1623−1641. (64) Jorgensen, W. L.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulationg Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (65) Izaguirre, J. A.; Catarello, D. P.; Wozniak, J. M.; Skeel, R. D. Langevin Stabilization of Molecular Dynamics. J. Chem. Phys. 2001, 114, 2090−2098. (66) Ryckaert, J. P.; Ciccotti, G.; Berendsen, H. J. C. Numerical Integration of the Cartesian Equations of Motion of a System with Constraints: Molecular Dynamics of N-Alkanes. J. Comput. Phys. 1977, 23, 327−341. 9551



Article

(67) Darden, T.; York, D.; Pedersen, L. Particle Mesh Ewald: An N· Log(N) Method for Ewald Sums in Large Systems. J. Chem. Phys. 1993, 98, 10089−10092. (68) Shao, J. Y.; Tanner, S. W.; Thompson, N.; Cheatham, T. E. Clustering Molecular Dynamics Trajectories: 1. Characterizing the Performance of Different Clustering Algorithms. J. Chem. Theory Comput 2007, 3, 2312−2334. (69) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong, L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; et al. Calculating Structures and Free Energies of Complex Molecules: Combining Molecular Mechanics and Continuum Models. Acc. Chem. Res. 2000, 33, 889−897. (70) Srinivasan, J.; Cheatham, I. T. E.; Cieplak, P.; Kollman, P. A.; Case, D. A. Continuum Solvent Studies of the Stability of DNA, Rna, and Phosphoramidate−DNA Helices. J. Am. Chem. Soc. 1998, 120, 9401−9409. (71) Wang, W.; Donini, O.; Reyes, C. M.; Kollman, P. A. Biomolecular Simulations: Recent Developments in Force Fields, Simulations of Enzyme Catalysis, Protein-Ligand, Protein-Protein, and Protein-Nucleic Acid Noncovalent Interactions. Annu. Rev. Biophys. Biomol. Struct. 2001, 30, 211−243. (72) Liu, H.; Bao, W.; Ding, H.; Jang, J.; Zou, G. Binding Modes of Flavones to Human Serum Albumin: Insights from Experimental and Computational Studies. J. Phys. Chem. B 2010, 114, 12938−12947. (73) Ersmark, K.; Nervall, M.; Hamelink, E.; Janka, L. K.; Clemente, J. C.; Dunn, B. M.; Blackman, M. J.; Samuelsson, B.; Aqvist, J.; Hallberg, A. Synthesis of Malarial Plasmepsin Inhibitors and Prediction of Binding Modes by Molecular Dynamics Simulations. J. Med. Chem. 2005, 48, 6090−6106. (74) Su, P.; Li, H. Energy Decomposition Analysis of Covalent Bonds and Intermolecular Interactions. J. Chem. Phys. 2009, 131, 014102. (75) Kitaura, K.; Morokuma, K. New Energy Decomposition Scheme for Molecular-Interactions within Hartree-Fock Approximation. Int. J. Quantum Chem. 1976, 10, 325−340. (76) Morokuma, K.; Kitaura, K. Chemical Applications of Atomic and Molecular Electrostatic Potentials; Plenum Press: New York, 1981. (77) Ziegler, T.; Rauk, A. Co, Cs, N2, Pf3, and Cnch3 as SigmaDonors and Pi-Acceptors - Theoretical-Study by the Hartree-FockSlater Transition-State Method. Inorg. Chem. 1979, 18, 1755−1759. (78) Hayes, I. C.; Stone, A. J. An Intermolecular Perturbation-Theory for the Region of Moderate Overlap. Mol. Phys. 1984, 53, 83−105. (79) Chen, W.; Gordon, M. S. Energy Decomposition Analyses for Many-Body Interaction and Applications to Water Complexes. J. Phys. Chem. 1996, 100, 14316−14328. (80) Schmidt, M. W.; Baldridge, K. K.; Boatz, J. A.; Elbert, S. T.; Gordon, M. S.; Jensen, J. H.; Koseki, S.; Matsunaga, N.; Nguyen, K. A.; Su, S. J.; et al. General Atomic and Molecular Electronic-Structure System. J. Comput. Chem. 1993, 14, 1347−1363. (81) Boys, S. F.; Bernardi, F. Calculation of Small Molecular Interactions by Differences of Separate Total Energies - Some Procedures with Reduced Errors. Mol. Phys. 1970, 19, 553−&. (82) Meng, E. C.; Pettersen, E. F.; Couch, G. S.; Huang, C. C.; Ferrin, T. E. Tools for Integrated Sequence-Structure Analysis with Ucsf Chimera. BMC Bioinform. 2006, 7, 339. (83) Tomasi, J.; Mennucci, B.; Cammi, R. Quantum Mechanical Continuum Solvation Models. Chem. Rev. 2005, 105, 2999−3093. (84) Barone, V.; Cossi, M. Quantum Calculation of Molecular Energies and Energy Gradients in Solution by a Conductor Solvent Model. J. Phys. Chem. A 1998, 102, 1995−2001. (85) Bostrom, J.; Norrby, P. O.; Liljefors, T. Conformational Energy Penalties of Protein-Bound Ligands. J. Comput. Aided Mol. Des. 1998, 12, 383−396. (86) Perola, E.; Charifson, P. S. Conformational Analysis of DrugLike Molecules Bound to Proteins: An Extensive Study of Ligand Reorganization Upon Binding. J. Med. Chem. 2004, 47, 2499−2510. (87) Siebel, C. W.; Kollman, P. A.; Pergamon: Oxford, 1990; Vol. 4.

(88) Ohtaka, H.; Schon, A.; Freire, E. Multidrug Resistance to Hiv-1 Protease Inhibition Requires Cooperative Coupling between Distal Mutations. Biochemistry 2003, 42, 13659−13666. (89) Perryman, A. L.; Lin, J. H.; McCammon, J. A. Hiv-1 Protease Molecular Dynamics of a Wild-Type and of the V82f/I84v Mutant: Possible Contributions to Drug Resistance and a Potential New Target Site for Drugs. Protein Sci. 2004, 13, 1108−1123. (90) Zoete, V.; Michielin, O.; Karplus, M. Relation between Sequence and Structure of Hiv-1 Protease Inhibitor Complexes: A Model System for the Analysis of Protein Flexibility. J. Mol. Biol. 2002, 315, 21−52. (91) Altman, M. D.; Ali, A.; Reddy, G. S.; Nalam, M. N.; Anjum, S. G.; Cao, H.; Chellappan, S.; Kairys, V.; Fernandes, M. X.; Gilson, M. K.; et al. Hiv-1 Protease Inhibitors from Inverse Design in the Substrate Envelope Exhibit Subnanomolar Binding to Drug-Resistant Variants. J. Am. Chem. Soc. 2008, 130, 6099−6113. (92) Ohtaka, H.; Velazquez-Campoy, A.; Xie, D.; Freire, E. Overcoming Drug Resistance in Hiv-1 Chemotherapy: The Binding Thermodynamics of Amprenavir and Tmc-126 to Wild-Type and Drug-Resistant Mutants of the Hiv-1 Protease. Protein Sci. 2002, 11, 1908−1916. (93) Stoica, I.; Sadiq, S. K.; Coveney, P. V. Rapid and Accurate Prediction of Binding Free Energies for Saquinavir-Bound Hiv-1 Proteases. J. Am. Chem. Soc. 2008, 130, 2639−2648.

9552


Systematic Molecular Dynamics, MMâPBSA, and Ab Initio Approaches

Recommend Documents

Systematic Molecular Dynamics, MMâPBSA, and Ab Initio Approaches