Binding Isotope Effects as a Tool for Distinguishing Hydrophobic and

Aug 16, 2014 - different inhibitors of HIV-1 RT, one of the three key proteins responsible ...... This work was supported by the Spanish Ministerio de...
0 downloads 0 Views 5MB Size
Article pubs.acs.org/JPCB

Binding Isotope Effects as a Tool for Distinguishing Hydrophobic and Hydrophilic Binding Sites of HIV‑1 RT ,†,§ ́ Agnieszka Krzemińska,† Piotr Paneth,† Vicent Moliner,*,‡ and Katarzyna Swiderek* †

Institute of Applied Radiation Chemistry, Lodz University of Technology, 90-924 Lodz, Poland Departament de Química Física i Analítica, Universitat Jaume I, 12071 Castelló, Spain § Departament de Química Física, Universitat de València, 46100 Burjassot, Spain ‡

S Supporting Information *

ABSTRACT: The current treatment for HIV-1 infected patients consists of a cocktail of inhibitors, in an attempt to improve the potency of the drugs by adding the possible effects of each supplied compound. In this contribution, nine different inhibitors of HIV-1 RT, one of the three key proteins responsible for the virus replication, have been selected to develop and test a computational protocol that allows getting a deep insight into the inhibitors’ binding mechanism. The interaction between the inhibitors and the protein have been quantified by computing binding free energies through FEP calculations, while a more detailed characterization of the kind of inhibitor−protein interactions is based on frequency analysis of the ligands in the initial and final state, i.e. in solution and binding the protein. QM/MM calculation of heavy atoms (13C, 15N, and 18O) binding isotope effects (BIE) have been used to identify the binding sites of the different inhibitors. Specific interactions between the isotopically labeled atoms of the inhibitors and polar residues and magnesium cations on the hydrophilic pocket of the protein are responsible for the frequencies shifting that can be detected when comparing the IR spectra of the compounds in solution and in the protein. On the contrary, it seems that changes in vdW interactions from solution to the final state when the ligand is interacting with residues of the hydrophobic cavity, does not influence frequency modes and then no BIE are observed. Our results suggest that a proper computational protocol can be a valuable tool which in turn can be used to increase the efficiency of anti AIDS drugs.

1. INTRODUCTION HIV-1 RT is one of the enzymes, together with protease (PR) and integrase (IN), involved in the human immunodeficiency virus (HIV) replication, which explains the fact that 13 of the 26 individual anti-AIDS approved antiretroviral drugs target reverse transcriptase (RT).1 HIV-1 RT is a heterodimer composed of two different chains that is responsible for copying the single-stranded viral RNA genome into doublestranded DNA.2 HIV-1 RT contains two main domains with two enzymatic activities: a DNA polymerase domain and a ribonuclease H (RNase H) domain (see Figure 1). The former is responsible for copying either an RNA or DNA, while the function of the RNase H domain is to cleave and degrade the template RNA after DNA synthesis so that the new DNA can generate a second DNA strand. The RNase H domain is also responsible for the integration of the duplex DNA into the host cell chromosome.3 Nucleoside HIV-1 RT inhibitors (NRTIs) bind the polymerase active site and, despite having received much attention in the drugs development,4 the fact that they also compete with natural nucleoside substrates of human host polymerases induces toxicity that causes serious adverse effects ranging from lactic acidosis to hepatic steatosis during longterm use as anti-AIDS drug.5,6 During the last years, and © XXXX American Chemical Society

probably due to these side effects, there has been considerable interest in developing non-nucleoside HIV-1 RT inhibitors (NNRTIs) that specifically inhibit the RNase H activity of RT by binding to a hydrophobic pocket in the palm subdomain adjacent to the base of the thumb subdomain. This excludes any binding interaction with other polymerases.7 As observed in Figure 1, both sites are clearly identified and located in different domains of the protein. The hydrophilic polymerase site is characterized by the presence of polar residues and two magnesium cations, and it is exposed to the solvent water molecules. The hydrophobic pocket is located in the palm of the p66 subunit of the enzyme close to the polymerase active site (∼10 Å). On the basis of the close proximity of the two binding sites new compounds have been proposed joining the two kind of drugs through a linker.8,9 Nevertheless, together with the problem of undesired side-effects of the long-term treatments, patients develop HIV-1 RT mutations even after a Special Issue: William L. Jorgensen Festschrift Received: June 19, 2014 Revised: August 11, 2014

A

dx.doi.org/10.1021/jp506119h | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 1. Representation of HIV-RT with bound DNA:RNA helix (in gold) and the residues forming the two binding sites, hydrophobic and hydrophilic, highlighted in cyan balls. The two binding cavities are amplified and docking two inhibitors, 3-cyclopentyl-1,4-dihydroxy-1,8naphthyridin-2(1H)-one and Nevirapine, which are displayed in licorice. Water molecules around the hydrophilic cavity are also displayed in licorice.

single dose, provoking drugs to become inactive. Consequently, discovery of new inhibitors without side effects and effective against existing variant strains of the virus is an active field of research. Estimation of binding affinities of different NNRTIs inhibitors has been the subject of previous theoretical studies by Jorgensen and co-workers based on Monte Carlo simulations.10−12 Interaction energies to the different residues of the allosteric binding pocket have been identified by MMPBSA combined with molecular docking,13 by quantum methods14 and by structure energy optimizations with ONIOM methods.15,16 More recently, Warshel and co-workers have evaluated the performance of PDLD/S-LRA/b method by computing the absolute binding free energies of a set of NNRT inhibitors.17 A contribution to the field has been done at our laboratory by a recent theoretical study of inhibition of HIV-1 RT by five different NRTIs based on the alchemical free energy perturbation method, FEP, and by a pathway method, in which the ligands were physically pulled away from the binding site.18 Application of FEP methods to biological processes, and in particular to studies of proteins, was first introduced by Warshel and co-workers,19 who obtained absolute binding free energies by means of thermodynamic cycles.20 According to our own experience with HIV-1 RT18 a simplified thermodynamic cycle can be used to compute enzyme−ligands binding free energies from alchemy free energy perturbation methods. These methods, while not providing detailed information about the full process of binding, render reliable values of binding free energy at a reasonable computational cost. In this study we are presenting a new additional tool based on calculation of heavy atoms binding isotope effects (BIEs) of NRTIs and NNRTIs that can be used to get a deep insight into the HIV-1 RT inhibitory process that allows distinguishing the site of the protein the inhibitors are binding to. Isotopic

substitutions of atoms involved in chemical reactions are traditionally used as a potent tool in many fields of science and technology.21 The change in the mass alters the vibrational frequencies while the underlying potential energy surface (PES) remains unaltered according to the Born−Oppenheimer approximation. Thus, accurate measurements of the relative reaction rate constants involving reactants that only differ in the masses of selected atoms, kinetic isotope effects (KIEs), can provide information about the nature of the chemical transition state. Although less studied, application of this physical concept to equilibrium steps can provide information on the relative population of two stable states in thermodynamic equilibrium.22 These equilibrium isotope effects (EIEs), when applied to the binding step of a substrate or ligand to an enzyme or a nonbiological host system, are known as binding isotope effects (BIE). Changes in normal-mode force constants surrounding the isotopic substitution between free (usually dissolved in water) and bound substrates give rise to a binding preference for light or heavy isotope in the bound state and hence a measurable BIE.23 Considering the process of binding of an inhibitor, from solution to the formation of the inhibitor− enzyme complex, BIE would be defined as the ratio between the binding equilibrium constant of the light (L), or natural abundant isotope species, and the heavy (H) or isotopically substituted one, BIE =

KL KH

(1)

BIEs can also be expressed in terms of the binding free energy of the light and heavy species, considering the relationship between the equilibrium constant and the binding free energy: BIE = e−(ΔG L −ΔG H)/RT B

(2)

dx.doi.org/10.1021/jp506119h | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

In eq 2, ΔGL and ΔGH are the binding free energies for the light and heavy istopomers, respectively, R is the universal gas constant and T the temperature. Then, as depicted in Figure 2,

protocol; four of them are NRTIs binding the hydrophilic site, while the remaining five inhibitors are the NNRTIs that the US Food and Drug Administration (FDA) has approved as drugs, which bind the hydrophobic active site of RNaseH. The interaction between the ligands and the protein will be quantified by means of FEP calculations while a more detailed characterization of the kind of ligand−protein interactions will be based on frequency analysis of the ligands in the initial and final state, i.e. in solution and in the active site of the protein.

2. COMPUTATIONAL METHODS System Setup. The four selected structures of NRTIs correspond to 2,7-dihydroxy-4-1(methylethyl)-2,4,6-cycloheptatrien-1-one or β-thujaplicinol (LNA-1), 3-cyclopentyl-1,4dihydroxy-1,8-naphthyridin-2(1H)-one (LNA-2), (ethyl 1,4dihydroxy-2-oxo-1,2-dihydro-1,8-naphthyridine-3-carboxylate (LNA-3) and 3-[4-(diethylamino)phenoxy]-6-(ethoxycarbonyl)5,8-dihydroxy-7-oxo-7,8-dihydro-1,8-naphthyridin-1-ium (LNA4). The corresponding protein−inhibitor complexes were prepared as described in our previous paper.18 The structures of HIV-1 RT with NNRTI were obtained from Protein Data Bank as structures with PDB IDs: 3V8125− LNN-1 (nevirapine), 1FK926− LNN-2 (efavirenz), 2ZD127− LNN-3 (rilpivirine), 1S9E28− LNN-4 (etravirine), and 1KLM29− LNN-5 (delaviridine). For all five structures the pKa for titratable amino acids was calculated using PROPKA software.30−32 In order to neutralize the systems Cl− contraions were added (22 for LNN1, LNN-4 and LNN-3, 20 for LNN-2 and 17 for LNN-5) and then systems were soak in 140 × 100 × 100 Å3 box of water molecules using NAMD33 program. All missing parameters for inhibitors, in order to run long MM MD simulations were obtained at AM1 level using Atechamber software implemented in AMBERTools.34 In order to adjust the system to the AMBER force field, to relax possible steric clashes, and to distribute water molecules, energy minimizations were carried out by means of a conjugate gradient algorithm. Subsequently, the systems were heated by increasing temperature from 0 to 300 K with temperature increment of 0.001 K. Then 2 ns Langevin−Verlet (NVT) MM MD simulation using AMBER force field implemented in NAMD was done. During the MM MD simulations all atoms were free to move, facilitating the structural changing of enzyme induced by influence of ligand. Periodic boundary conditions with cut-offs for the nonbonding interactions were applied using a switching scheme, within a range radius from 14.5 to 16 Å. After long MM MD relaxation, 200 ps of QM/MM MD simulations (at 300 K) at AM135/MM level were done. During these simulations, the atoms of the ligands were selected to be treated by QM. The rest of the system (protein plus water molecules) were described using the AMBER36−38 and TIP3P39 force fields, respectively, implemented in the fDYNAMO library.40,41 In order to reduce the cost of the calculations, the protein atoms and all water molecules 20 Å apart from any of the atoms of the initial inhibitor were kept frozen in the remaining calculations with the same switching function described above to treat the nonbonding interactions. Binding Isotope Effects (BIEs). Eleven structures of each inhibitor (five NNRTIs and four NRTIs) in aqueous solution and binding the corresponding protein site were extracted from the last 20 ps of QM/MM MD simulations to compute semiclassical binding isotope effects, BIEs, without scaling the vibrational frequencies, at AM1/MM (for all tested inhibitors) and B3LYP/6-31+G(d,p)/MM levels (for LNA-1 and LNN-1).

Figure 2. Schematic representation of the origin of inverse BIEs (stronger confining well, higher force constants in the enzyme than in solution) due to contribution of zero-point vibrational energies in a binding process of an inhibitor from aqueous solution to an enzyme cavity. The molecular potential energy curves have been approximated by parabolas, according to the harmonic approximation.

assuming that only zero-point vibrational energies (ZPVE) are isotope sensitive and the harmonic approximation, inverse BIEs (values lower than unity), indicate stronger interactions of the heavier isotopically substituted atoms in the enzyme cavity (higher force constants and then steeper parabola) than in solution, i.e. higher (more negative) values of the free energy of binding, then BIEs can be used to characterize the increase or decrease in binding upon substitution of an atom with its heavier isotope,23 or in other words, they can render information on the protein−ligand interaction. A combination of experimental measurements and theoretical computation of BIE can be used to identify binding sites on a protein. BIEs were already successfully used as a tool in analyzing interactions between ligand and active sites of enzymes, binding site of carrying proteins, or cavities of nonbiological host systems.22 Since increased atomic mass alters the bond vibrational environment of a molecule, the binding equilibrium of a ligand between aqueous solution and the protein cavity is sensitive to isotopic substitution, and then, the magnitude of the obtained BIE can render an estimation of the interactions established in both media.23 This can be of special interest in those proteins presenting more than one possible receptor site. Recently, BIE based approaches, together with other experimental techniques, have been used to address the question of tautomerism in the context of an RNA aptamer−ligand complex.24 Nevertheless, the fact that calculations were performed in gas phase or within a continuum model precludes the effect of specific interactions to be incorporated in the simulations, and thus, a correlation with the almost indistinguishable BIEs obtained for the different tautomers was very difficult. In the present study, BIEs are going to be computed with hybrid Quantum Mechanics/ Molecular Mechanics (QM/MM) methods, thus including the effect of the environment (protein and solvent water molecules) in an explicit way. Keeping in mind that the current treatment for HIV infected patients consists of a cocktail of inhibitors, and the recent efforts in developing bifunctional RT inhibitors, a detailed knowledge of the binding mechanism can be a valuable tool to increase the efficiency of the drugs by designing more selective compounds. Nine different inhibitors have been selected to develop and test our C

dx.doi.org/10.1021/jp506119h | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

⎛ EQM/MM(λ , γ ) = ⟨Ψ|Ĥ 0|Ψ⟩ + λ⎜⎜∑ ⎝

The subset of atoms used to define the Hessian for these BIE calculations were those of the QM region consistent with the “cut-off rule” and the local nature of isotope effects. From the definition of the free energy of a state, Gi, as a function of the internal energy Ui, the total partition function Qi, and the zero point vibrational energy, ZPEi, Gi = Ui − RT ln Q i + ZPE i

+

qMM re,MM

Ψ

ZQMqMM ⎞ ⎟ + γE vdW QM/MM + E MM rQM,MM ⎟⎠

(6)

According to Scheme 1, the binding free energy of a particular ligand would be expressed as

(3)

the free energy of binding of a inhibitor from solution “w” to the cavity of the enzyme “e” can be expressed according to eq 4,

Scheme 1. Thermodynamic Cycle to Compute Enzyme− Inhibitor Binding Free Energies from Alchemy Free Energy Perturbation Methodsa

ΔG = Ge − Gw = (Ue − Uw ) + ( − RT ln Q e + RT ln Q w) + (ZPEe − ZPE w ) = ΔU + RT ln

∑∑

Ψ

Qw Qe (4)

+ ΔZPE

Then, the binding isotope effect of substituting a light atom “L” by a heavier isotope “H” can be computed by substituting eq 4 into eq 2.;

( ) BIE = ( ) Qe

Qw

L

Qe

Qw

H

e−1/ RT(ΔZPEL −ΔZPEH) (5)

In eq 5, the total partition function, Q, was computed as the product of the translational, rotational, and vibrational partition functions for the isotopologs inhibitors in solution and in the enzyme. The Born−Oppenheimer, rigid-rotor, and harmonic oscillator approximations were considered to independently compute the different contributions. The full 3N × 3N Hessians have been subjected to a projection procedure to eliminate translational and rotational components, which give rise to small nonzero frequencies, as previously described.42 Thus, it has been assumed that the 3N − 6 vibrational degrees of freedom are separable from the 6 translational and rotational degrees of freedom of the inhibitor in both the binding site and in water. It is important to point out that, as shown by Williams,43 an alternative procedure where all 3N degrees of freedom of the inhibitors are considered as vibrations, thereby including the (possible) isotopic sensitivity of the librational motions of the inhibitor relative to its environment could involve smaller errors. These errors can be significant when a reduced number of atoms (cutoff model) is used to compute the Hessians. Hessians have been also used to compute IR spectra. Finally, 121 BIEs values obtained with all possible combinations of structures in water and in the enzyme (11 × 11) were averaged for each isotope substitution. Alchemical Free Energy Perturbation (FEP) Methods. Ligand−protein interaction free energies for the NNRTIs were computed as described in ref 18 for NRTIs. Thus, series of QM/MM MD simulations have been carried out in the hydrophobic cavity of the protein and in a box of solvent water molecules, introducing two parameters; λ and γ parameters in the electrostatic and van der Waals QM/MM interaction terms, respectively:

a

E is the enzyme with inhibitor (I) in its binding site, E′ is the apo form of enzyme, (I)0 is the ligand in gas phase. elec vdW ΔΔG bind = ΔGw − ΔG E = Δ(GQM/MM + GQM/MM )w elec vdW − Δ(GQM/MM + GQM/MM )E

(7)

where the two terms of the right side of eq 7 are computed by applying eqs 8 and 9, in water and in the enzyme, respectively. During simulations the λ smoothly changes between values of 1, corresponding to a full MM charge-wave function interaction, and 0; whereas the electrostatic interaction with the force field is not considered. Once the electrostatic charges of the ligand (QM region) are annihilated, subsequently γ smoothly changes between values of 1, corresponding to a full QM/MM van der Waals interaction, and 0 where there is no van der Waals interaction with the force field. The calculation of the free energy difference of two consecutive windows of the two stages, annihilation of charges and van der Waals parameters, is performed by means of FEP methods, and then the total free energy change is evaluated as the sum of all the windows covering the full transformation from the initial to the final state. In practice, this is done in two steps (eqs 8 and 9), once the charges are annihilated then the van der Waals: 1 elec ΔGQM/MM = − [∑ ln⟨e−β[E(λi+1) − E(λi)⟩λi ]γ = 1 β

(8)

1 vdW ΔGQM/MM = − [∑ ln⟨e−β[E(λi+1) − E(λi)⟩γi ]λ = 0 β

(9)

We have used a total of 50 windows to evaluate the electrostatic interaction term, from λ = 0 (no electrostatic D

dx.doi.org/10.1021/jp506119h | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Table 1. Electrostatic (ΔGe) and van der Waals (ΔGvdW) Terms of the Water−Ligand and Protein−Ligand Interaction Free Energy, and Electrostatic, van der Waals, and Total Free Energy Differences of Interaction between Water and the Binding Site of the Protein (ΔΔGe, ΔΔGvdW, and ΔGbind)a ΔGe‑w

ΔGe‑HIV

ΔΔGe

ΔGvdW‑w

ΔGvdW‑HIV

ΔΔGvdW

ΔGbind

−16.8 −23.7 −23.1 −36.5

−19.6 −28.6 −23.8 −39.2

−2.8 −4.9 −0.7 −2.7

−27.9 −17.8 −30.2 −19.1

−22.1 −21.3 −29.5 −30.3 −37.1

−38.4 −37.5 −50.9 −53.8 −62.0

−16.3 −16.2 −21.4 −23.5 −24.9

−20.7 −19.6 −29.2 −27.6 −25.8

NRTIs LNA-1 LNA-2 LNA-3 LNA-4

−5.0 −12.0 −8.2 −13.4

−30.1 −24.9 −37.7 −29.8

−25.1 −12.9 −29.5 −16.4

LNN-1 LNN-2 LNN-3 LNN-4 LNN-5

−7.2 −8.1 −7.1 −8.1 −28.9

−11.6 −11.5 −14.9 −12.2 −29.8

−4.4 −3.4 −7.8 −4.1 −0.9

NNRTIs

a

All values are reported in kcal·mol−1. Values of NRTIs come from ref 18.

Figure 3. Contributions of electrostatic (gray bars) and vdW (blue bars) interactions to the difference of total free energy of binding of the different inhibitors to the corresponding cavity (ΔΔGw − ΔΔGe): (A) hydrophilic RNaseH active site and (B) hydrophobic cavity. Values are in percentages.

mol−1. On the contrary, the analysis of the values obtained for the NNRTIs shows that the electrostatic term of the protein− inhibitor free energy of interaction are quite close to the electrostatic term of the waters−inhibitor free energy of interaction in aqueous solution. Values obtained for LNN-1, LNN-2, LNN-3, and LNN-4 range from −7 to −8 kcal·mol−1 in solution and from −12 to −15 kcal·mol−1 in the enzyme. These values render differences between enzyme and water that range from 3 to 8 kcal·mol−1. Interestingly, the van der Waals terms of the protein−ligand free energy of interaction are markedly higher than the corresponding values obtained for the waters− ligand interactions in solution. In this case, not only the absolute values are higher but also the differences between both media, ranging from −16 to −24 kcal·mol−1. The comparative analysis can be clearly observed in Figure 3, where differences of contributions of electrostatic and vdW interactions to the total free energy of binding of all studied inhibitors between water and the corresponding cavity of the protein (ΔΔGw − ΔΔGe) are presented as a bars plot. Clearly, electrostatic interactions must be the driving force of the docking process for the NRTIs, while binding of NNRTIs to the hydrophobic pocket must be due to nonelectrostatic interactions. LNN-5 presents a peculiar behavior. As observed in Table 1, the values of the van der Waals, and especially the electrostatic terms of the free energy of interaction in the enzyme and in solution, are much higher than the rest of the NNRTIs. In fact this inhibitor presents values in the range of the ones obtained for the NRTIs. A closer inspection of the interactions exhibited in both media reveals that the methylsulfonyl amino group of this

interaction) to λ = 1 (full interaction), which turns into a δλ of 0.02. The same number of windows has been used to calculate the van der Waals interaction term, from β = 0 (no interaction) to γ =1 (full interaction), which turns into a δγ of 0.02. The procedure can be applied both forward (i → i + 1) and backward (i + 1 → i) as long as the comparison provides information about the convergence of the process. In each window a total of 5 ps of relaxation followed by 50 ps of production QM/MM MD have been performed using the NVT ensemble at the reference temperature of 300 K.

3. RESULTS AND DISCUSSION The ligand−protein binding free energies and their decomposition into van der Waals and electrostatic terms computed from FEP methods are reported in Table 1. As expected, binding of NRTIs to the hydrophilic cavity is basically due to favorable electrostatic interactions in the enzyme with respect to the aqueous media, while the van der Waals interactions are mostly responsible for the binding of NNRTIs in the hydrophobic cavity (see ΔΔGe and ΔΔGvdW columns in Table 1). Thus, the electrostatic term of the free energy of the NRTIs−RT interaction ranges from −30 to −32 kcal·mol−1, values much higher than the corresponding values in aqueous solution that range from −5 to −13 kcal·mol−1. The resulting electrostatic contribution to the binding process appears to be quite significant, with differences between aqueous and protein environment ranging from −13 to −30 kcal·mol−1. van der Waals terms in solution and in the enzyme cavity, however, are quite similar, providing differences between −1 and −5 kcal· E

dx.doi.org/10.1021/jp506119h | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 4. Computed heavy atoms BIEs for nucleoside analog RT inhibitors bonded to the hydrophilic cavity, top, and non-nucleoside RT inhibitors bonded to the hydrophobic cavity of HIV-1 RT. Inverse BIEs are reported in red, normal BIEs in blue, and isotopic substitutions rendering no effects are in black.

water and protein, the values are in agreement with the rest of studied NNRTIs (see differences of interaction energies in Table 1 and the graphical representation in Figure 3). Keeping in mind the two types of inhibitors under study, the results of the contributions of electrostatic and vdW interactions to the difference of total free energy of binding of the NRTIs and NNRTIs to their corresponding cavity, as shown in Figure 3, appear to be quite reasonable. Binding of

butterfly-like shape class of compound is strongly interacting with water molecules not only (and obviously) when it is free solvated in solution but also when it is bonded to the hydrophobic pocket, where this group of the inhibitor is exposed to the solvent. A representation of this and the rest of NNRTIs inhibitors bonded to the hydrophobic pocket of the RT protein is depicted in Supporting Information (SI) (Figure S1). Nevertheless, when doing a comparative analysis between F

dx.doi.org/10.1021/jp506119h | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 5. (a) Schematic representation of the interactions of key oxygen atoms of LNA-1 with water molecules of the aqueous solution and with residues of the hydrophilic and hydrophobic sites of HIV-1 RT, respectively. (b) IR spectra of LNA-1 in water and in the enzyme (top and middle panels) and superposition of IR spectra in water (blue line) and in the enzyme (black line) of LNN-1 (bottom panel). All figures and data are based on B3LYP/MM calculations.

of binding free energies reported in Table 1 for the same inhibitors (LNN-4, LNN-1, and LNN-2): −27.6, −20.7, and −19.6 kcal·mol−1, respectively. Once the binding free energies are analyzed, the next step is the study of the BIEs. The results of heavy atoms BIEs (13C, 15 N, and 18O) corresponding to the equilibrium of the four NRTIs and the five NNRTIs established between aqueous solution and the hydrophilic and hydrophobic cavities of HIV-1 RT, respectively, are reported in Figure 4 (see the decomposition of the different terms contributing to the BIEs for some of the most significant isotope substitutions in LNA-1 in Table S1 in the SI). According to the BIEs results, the first conclusion that can be obtained is the lack of BIEs for most of the substitutions on the NNRTIs. Interestingly, it is noticeable that the normal 18O BIEs observed in the carbonyl oxygen of LNN-1 (1.011 ± 0.002) is a result that can be rationalized on the basis of stronger interactions with the water molecules in solution than when bonded to the hydrophobic site. Other 18O substitutions on other oxygen atoms of NNRTIs, ether or carbonyl oxygen atoms, do not render BIEs within the statistical deviations. This is due to the fact that, while there are interactions that are established with water molecules in solution, the labeled atom is also exposed to solvent water molecules when bonded into the cavity and to the N−H bond of the backbone of residues such as Lys101 (in the case of carbonyl oxygen atom in LNN-2) or Lys103 (in the case of carbonyl oxygen atom in LNN-5). In contrast, many non-

NRTIs would be governed by electrostatic interactions with the protein, while binding of NNRTIs to the hydrophobic pocket of the protein would be based on van der Waals interactions. The validation of the absolute binding free energies of the NNRTIs to the hydrophobic pocket, by comparison with experimental data, is more problematic. As mentioned in our previous paper,18 in the context of the analysis of our predicted values for NRTIs, the use of IC50 (value which expresses the concentration of an inhibitor required to produce 50% inhibition of an enzymatic reaction at a specific substrate concentration) would not be adequate since it measures other steps of the full inhibition processes. The inhibitory binding constant of a competitive inhibitor to the enzyme, or the reciprocal of the dissociation constant of the enzyme−inhibitor complex (KI), should be a more appropriate value that could be straightforward compared with binding free energy values.44 Nevertheless, there is no single paper where a systematic study of KIs of the different NNRTIs, measured at the same experimental conditions, was reported. Keeping these limitations in mind, the trend of the computed binding free energies of the nine inhibitors appear to be in reasonable qualitative agreement with experimental evidence: secondgeneration NNRTIs show higher binding energies than the first-generation compounds. Moreover, in a recent study of Agneskog et al.,45 the measurements of IC50 for Etravirine, Nevirapine, and Efavirenz on wild-type HIV-1 RT recovered from plasma of HIV-1 infected patients were 0.9, 1.1, and 3.6 μM, respectively. These values are in agreement with the trend G

dx.doi.org/10.1021/jp506119h | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

interactions in the protein cavity than with water molecules in solution. By comparison with LNA-2, a compound that has the same ester group, the posing of both compounds in the active site of the protein is significantly different. Thus, while this moiety of LNA-4 is more exposed to the water molecules of the solvent, the ether oxygen atom of LNA-2 is surrounded by the backbone atoms of residues His539 and Asp549. This is probably the origin of a less inverse 18O BIE (0.998 vs 0.989 in LNA-2 and in LNA-4, respectively). A stronger intramolecular hydrogen bond between the carbonyl oxygen of the ester group in LNA-2 would explain its inverse value, 0.990, in contrast with the normal 18O BIE obtained in LNA-4 (1.021). This hydrogenbond interaction in LNA-4 also explains the different values of 18 O BIE obtained for the involved hydroxyl group. A stronger hydrogen bond makes a less stiff bonded oxygen atom of the hydroxyl group, rendering a normal effect: 1.004 in LNA-2 and 0.994 in LNA-4. Arguments similar to the ones used to analyze 18O BIEs can be used to explain the 15N BIEs obtained in LNA-2, LNA-3, and LNA-4. In these inhibitors, lone electron pairs of nitrogen atoms appear to be interacting with the magnesium cations in the protein site, thus in agreement with the inverse effect. Regarding the 13C BIEs, normal effects are obtained in those atoms that are confined in the protein site without any specific interaction, while in solution they can interact, despite weakly, with water molecules through weak C−H···O hydrogen bonds. As observed, higher values are obtained in methyl groups since an additive effect of the three hydrogen atoms can take place (LNA-1, LNA-2, and LNA-4) and in carbon atoms of the nonaromatic five-membered ring of LNA-3. Other heavy atom BIEs cannot be directly analyzed due to more complex interactions. Moreover, in most of them the effects are not so dramatic and consequently not so easily experimentally detected, a priori. The IR spectra of the different compounds, and in particular the peaks corresponding to frequencies associated with oxygen atoms, can complement the study of 18O BIEs based on analysis of geometries. The comparison of the predicted IR spectra of LNA-1 in aqueous solution and in the hydrophilic site of HIV-1 RT, presented in Figure 5, reveals a noticeable shift in frequencies of those oxygen atoms that rendered measurable 18 O BIEs. It can be observed how the OC stretching of the carbonyl group and the C−O stretching of the hydroxyl groups shift from 1589.1 and 1321.8 cm−1, in water, to 1226.5 cm−1 in the protein (both values in the protein are indistinguishable). The O−H stretching shifts from 3633.3 cm−1 in water to 2727.8 cm−1 in the protein. According to the harmonic approximation higher values of frequencies should be associated with higher bond force constants and consequently stiffer vibrational modes in solution than in the enzyme. Nevertheless, in order to understand the shifts on the frequencies, we have to consider all interactions established with each functional group. Thus, the carbonyl group and hydroxyl groups of the LNA-1 interact with the magnesium cations, Asp549 and Glu478 when docked in the active site of the protein. These interactions provoke a displacement of the electron density from the bond to the species the oxygen atom is interacting with and, consequently, weaker CO, C−O and O−H interactions, thus rendering smaller frequencies. This can be supported by the analysis of the distances. Thus, CO distance is 1.27 and 1.30 Å in solution and in the protein, respectively, while change in O−H distance from solution to the enzyme is more dramatic; from 0.98 Å in solution to 1.01 and 1.04 Å in the enzyme (see

negligible values of BIEs have been obtained for the NRTIs. In order to analyze the results, it must be kept in mind that, within the harmonic approximation, inverse BIEs (lower than unity) mean that the isotopically substituted atom experiences a stronger interaction in the binding site of the protein than when solvated in aqueous solution (see Figure 2). The steeper the sides of the potential, the greater the force constant. This implies, in molecular terms, stiffer bonds or a stronger confining atom. Then, the fact that no significant values of BIEs are obtained in the NNRTIs is in agreement with no specific interactions between the ligands and the protein in the hydrophobic pocket, as remarked above. The molecular structure of these inhibitors shows almost no single strong polar group capable of discriminatorily interacting with ionic residues of the protein. The five NNRTIs present interactions with aromatic rings belonging to the π-box; Tyr181, Tyr188, Phe227, or Trp229, with hydrophobic residues such as Pro59, Leu100, Val106, Leu234, or Pro236, and weak interactions with some hydrophilic residues such as Lys101, Lys103, Ser105, Asp132, or Glu224. The intensity of these interactions with hydrophilic residues appear to be at the same order of magnitude as the interactions that the corresponding groups of the inhibitor exhibit with water molecules in solution. Then, the BIEs appear to be in agreement with this analysis. Quite the reverse, the computed values for some of the heavy atom BIEs obtained for the NRTIs suggest that they could be experimentally measured. Therefore, almost all hydroxyl groups of these compounds show inverse 18O BIEs ranging from 0.991 to 0.956. The highest 18O BIEs of the hydroxyl groups (0.966 in LNA-1, 0.988 in LNA-2, 0.956 in LNA-3, and 0.982 in LNA-4) are obtained in those cases when the oxygen atom of the hydroxyl group is strongly interacting with the magnesium cations of the active site. Figure 5 shows a schematic representation of the interactions of a fragment of the LNA-1 with water molecules in solution and with the magnesium cations and some residues in the hydrophilic site (representation of the rest of the NRTIs are reported in Figure S2 in the SI). The interactions with the divalent magnesium cations appear to be much stronger than the ones established with water molecules when the ligands are in solution. The other hydroxyl group of LNA-1 shows a less inverse 18O BIE, 0.993 (see Figure 4), which can be rationalized by taking into account that the hydrogen atom of the group is interacting with Asp549 when the ligand is bonded in the active site. This interaction polarizes the O−H bond and makes the interactions established with the oxygen atom less strong or, in other words, weaker confined wells of its molecular potential energy curve and, consequently, lower BIE. The hydrogen atom of the first hydroxyl group of LNA-1 is also interacting with a residue when bonded to the active site, in this case with Glu478, but as shown in Figure 5, not so strongly. Other hydroxyl groups in compounds LNA-2, LNA-3, and LNA-4 do not show 18O BIEs, within the deviation uncertainty, in agreement with the fact that no strong interactions are present in the active site, by comparison with the interactions in aqueous solution (see Figure S2 in the SI). The study of the rest of 18O BIEs, those on carbonyl oxygen atoms and the ether group of LNA-4, can be done following the same analysis. As observed, only a significant value is obtained in the case of LNA-1 (0.962), consistent with the strong interaction established with the magnesium cations in the protein cavity. The ether oxygen atom of LNA-4 presents a 18O BIE of 0.989. This is consistent, again, with stronger H

dx.doi.org/10.1021/jp506119h | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

Figure 6. Comparison of AM1/MM (black squares) and B3LYP/MM (blue circles) heavy atom BIEs of LNA-1 when binding to the hydrophilic site of HIV-1 RT (a) and of LNN-1 when binding to the hydrophobic site and of HIV-1 RT (b).

a couple of cases, DFT/MM calculations render higher values of BIEs (more inverse or more normal) than AM1/MM calculations which means that an easier analysis could be done. On the basis of this comparative analysis, our predictions based on AM1/MM calculations can be considered as robust.

Figure 5). Nevertheless, despite these interactions becoming weaker in the enzyme, the oxygen atoms are stiffer since they are constrained by the strong interatomic interactions established with the mentioned residues, much stronger than the possible interactions with water molecules in solution. Then, the measured inverse BIE can be explained. The shifts on the C−H frequencies are not so visible but, in any case, it is possible to observe the opposite trend; an increase from 3109.0 cm−1 in solution to 3126.8 cm−1 in the enzyme. In this case, the behavior is due to the lack of polar interaction in the hydrophilic pocket, by comparison with the weak interactions that can be established with the water molecules in solution. The different values of 18O BIEs of the oxygen atoms of carbonyl groups and hydroxyl groups of the different tested compounds can be rationalized on the basis of the frequency modes of the IR spectra, keeping in mind that a reduction in the force constant of the original bonds are compensated with additional bonds established in the protein site. Larger shifts to smaller values are associated with stronger protein−ligand interactions and, consequently, more inverse BIEs. Thus, inverse 18O BIE can be associated with displacements of vibrational frequencies (from the IR obtained in water to the one in the protein) to lower values. As commented above regarding the analysis of some heavy atom BIEs, analysis of variations in frequencies associated with other atoms such as carbon or nitrogen is not so straightforward. When analyzing the corresponding IR spectra of the NNRTIs in solution and in the hydrophobic cavity, no shifts are detected (see IR spectra of LNN-1 in Figure 5b). According to the previous arguments, this observation is in agreement with the fact that no BIEs are observed for the binding of this class of RT inhibitors. Finally, keeping in mind the possible limitations of the semiempirical AM1 Hamiltonian in describing frequencies, we have computed frequencies and heavy atom BIEs for LNA-1 and LNN-1 using DFT/MM methods, with the standard B3LYP functional. The 13C, 15N, and 18O BIEs computed with both methods are shown in Figure 6. We must take into account that AM1/MM BIEs are reported as average values within the standard deviations, while B3LYP/MM values come from only one structure in solution and one in the protein. Nevertheless, as observed in the figure, both methods render the same trend of BIEs for both ligands when binding their corresponding protein sites. It is important to stress that there is no single position of the compounds where AM1 and B3LYP gave different kind of BIEs (normal vs inverse). Moreover, except for

4. CONCLUSION Heavy atom (13C, 15N, and 18O) BIEs of the two classes of HIV-1 RT inhibitors, NRTIs and NNRTIs, have been computed by means of hybrid QM/MM methods for the binding process from aqueous solution to the corresponding site on the protein, hydrophilic and hydrophobic sites, respectively. Structures of the tested inhibitors bonded to the protein have been obtained after long QM/MM MD simulations, which allows optimizing the posing of the compounds in their respective protein site. The interaction between the inhibitors and the protein has been quantified by computing the free energy of binding from aqueous solution to the active site of the protein based on FEP calculations. The analysis of the results has confirmed the different nature of the protein−ligand interactions. Thus, while the electrostatic term is the most important contribution to the stabilization of the nucleoside analogue inhibitors in their corresponding protein site, the preferential stabilization of the non-nucleoside analogue inhibitors for the hydrophobic cavity remains in the van der Waals protein−ligand interactions. The results of heavy atom BIEs reveal that these magnitudes can be used to distinguish the cavity where the putative inhibitors are docking to. Thus, it seems that changes in vdW interactions from the initial state, when interacting with water molecules in solution, to the final state, when NNRTIs are interacting with residues of the hydrophobic cavity, do not influence frequency modes, and then no BIE are observed. On the contrary, when computing heavy atom BIEs for the NRTIs, and in particular 18O BIEs, significant values are obtained. Specific interactions between the isotopically labeled oxygen atoms of these inhibitors and polar residues and magnesium cations on the protein are responsible for the frequency shifts that can be detected when comparing the IR spectra of the compounds in solution and in the protein. Interestingly, if only the IR spectra are compared, an a priori contradictory result is observed since peaks corresponding to frequencies associated with functional groups such as carbonyl or hydroxyl are shifted to lower values. In order to correlate the IR spectra with the computed heavy atom BIEs, a detailed analysis of the total amount of interactions established with the functional group I

dx.doi.org/10.1021/jp506119h | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

nucleoside Inhibition Mechanism. Nat. Struct. Mol. Biol. 2012, 19, 253−259. (8) Smerdon, S. J.; Jager, J.; Wang, J.; Kohlstaedt, L. A.; Chirino, A. J.; Friedman, J. M.; Rice, P. A.; Steitz, T. A. Structure of the Binding Site for Non-nucleoside Inhibitors of the Reverse Transcriptase of Human Immunodeficiency Virus Type 1. Proc. Natl. Acad. Sci. U.S.A. 1994, 91, 3911−3915. (9) Bailey, C. M.; Sullivan, T. J.; Iyidogan, P.; Tirado-Rives, J.; Chung, R.; Ruiz-Caro, J.; Mohamed, E.; Jorgensen, W. L.; Hunter, R.; Anderson, K. S. Bifunctional Inhibition of Human Immunodeficiency Virus Type 1 Reverse Transcriptase: Mechanism and Proof-ofConcept as a Novel Therapeutic Design Strategy. J. Med. Chem. 2013, 56, 3959−3968. (10) Smith, R. H., Jr.; Jorgensen, W. L.; Tirado-Rives, J.; Lamb, M. L.; Janssen, P. A.; Michejda, C. J.; Kroeger Smith, M. B. Prediction of Binding Affinities for TIBO Inhibitors of HIV-1 Reverse Transcriptase Using Monte Carlo Simulations in a Linear Response Method. J. Med. Chem. 1998, 41, 5272−5286. (11) Rizzo, R. C.; Tirado-Rives, J.; Jorgensen, W. L. Estimation of Binding Affinities for HEPT and Nevirapine Analogs with HIV-1 Reverse Transcriptase via Monte Carlo Simulations. J. Med. Chem. 2001, 44, 145−154. (12) Jorgensen, W. L.; Bollini, M.; Thakur, V. V.; Domaoal, R. A.; Spasov, K. A.; Anderson, K. A. Efficient Discovery of Potent Anti-HIV Agents Especially Targeting the Tyr181Cys Variant of HIV Reverse Transcriptase. J. Am. Chem. Soc. 2011, 133, 15686−15696. (13) Wang, J.; Morin, P.; Wang, W.; Kollman, P. A. Use of MMPBSA in Reproducing the Binding Free Energies to HIV1 RT of TIBO Derivatives and Predicting the Binding Mode to HIV1 RT of Efavirenz by Docking and MM-PBSA. J. Am. Chem. Soc. 2001, 123, 5221−5230. (14) He, X.; Mei, Y.; Xiang, Y.; Zhang, D. W.; Zhang, J. Z. H. Quantum Computational Analysis for Drug Resistance of HIV-1 Reverse Transcriptase to Nevirapine through Point Mutations. Proteins: Struct, Funct. Bioinf. 2005, 61, 423−432. (15) Saen-oon, S.; Kuno, M.; Hannongbua, S. Binding Energy Analysis for Wild-type and Y181C Mutant HIV-1 RT/8-Cl TIBO Complex Structures: Quantum Chemical Calculations Based on the ONIOM Method. Proteins, Struct. Funct. Bioinf. 2005, 61, 859−869. (16) Nunrium, P.; Kuno, M.; Saen-oon, S.; Hannongbua, S. Particular Interaction between Efavirenz and the HIV-1 Reverse Transcriptase Binding Site as Explained by the ONIOM2 Method. Chem. Phys. Lett. 2005, 405, 198−202. (17) Singh, N.; Warshel, A. Absolute Binding Free Energy Calculations: On the Accuracy of Computational Scoring of Protein−ligand Interactions. Proteins: Struct, Funct. Bioinf. 2010, 78, 1705−1723. ́ (18) Swiderek, K.; Martí, S.; Moliner, V. Theoretical Studies of HIV1 Reverse Transcriptase Inhibition. Phys. Chem. Chem. Phys. 2012, 14, 12614−12624. (19) Warshel, A.; Sussman, F.; King, G. Free Energy of Charges in Solvated Proteins: Microscopic Calculations Using a Reversible Charging Process. Biochemistry. 1986, 25, 8368−8372. (20) Lee, F. S.; Chu, Z. T.; Bolger, M. B.; Warshel, A. Calculations of Antibody-Antigen Interactions: Microscopic and Semi-microscopic Evaluation of the Free Energies of Binding of Phosphorylcholine Analogs to McPC603. Protein Eng. 1992, 5, 215−228. (21) Kohen, A.; Limbach, H. H., Eds. Kinetic Isotope Effects as Probes for Hydrogen Tunneling in Enzyme Catalysis. Isotope Effects in Chemistry and Biology; CRC Press: Boca Raton, FL, 2006; Vol. 28, pp 743−764. ́ (22) Swiderek, K.; Paneth, P. Binding Isotope Effects. Chem. Rev. 2013, 113, 7851−7879. (23) Schramm, V. L. Binding Isotope Effects: Boon and Bane. Curr. Opin. Chem. Biol. 2007, 11, 529−536. (24) Singh, V.; Peng, C. S.; Li, D.; Mitra, K.; Silvestre, K. J.; Tokmakoff, A.; Essigmann, J. M. Direct Observation of Multiple Tautomers of Oxythiamine and their Recognition by the Thiamine Pyrophosphate Riboswitch. ACS Chem. Biol. 2014, 9, 227−236.

the atoms belong to (by bond polarization, rotational constraints, or hydrogen bonding) has to be considered. Thus, the stronger the total interaction is observed, the steeper the walls of the molecular potential energy curves and more inverse BIEs would be expected. The obtained results can be also important when KIEs are studied. In a more general sense, our results suggest that intrinsic KIE would be different from observed KIE when BIEs exist, which can be the case if ligands are binding into hydrophilic cavities. Keeping in mind that the current treatment for HIV-1 infected patients consists of a cocktail of inhibitors (most of them against HIV PR and HIV RT), a detailed knowledge of the binding mechanism and which cavity is being occupied by each inhibitor, can be a valuable tool to increase the efficiency of the anti-AIDS therapy by designing new compounds.



ASSOCIATED CONTENT

S Supporting Information *

Figure S1; representative snapshots of NNRTIs bonded to the hydrophobic pocket of RT protein. Figure S2; schematic representation of key interactions between LNA-2, LNA-3, and LNA-4 with the magnesium cations and some residues in the hydrophilic site. Table S1; contribution of different terms to the BIEs for key substitutions in LNA-1. Full reference 34. This material is available free of charge via the Internet at http:// pubs.acs.org.



AUTHOR INFORMATION

Corresponding Authors

*E-mail: [email protected]. *E-mail: [email protected]. Notes

The authors declare no competing financial interest.



ACKNOWLEDGMENTS This work was supported by the Spanish Ministerio de Economiá y Competitividad (ref CTQ2012-36253-C03), Universitat Jaume I (ref P1·1B2011-23) and the Polish National Science Center (NCN) for grant 2011/02/A/ST4/00246. The authors acknowledge the Servei d’Informatica, Universitat Jaume I, for generous allotment of computer time.



REFERENCES

(1) Das, K.; Arnold, E. HIV-1 Reverse Transcriptase and Antiviral Drug Resistance. Curr. Opin. Virol. 2013, 3, 111−118. (2) Telesnitsky, A.; Goff, S. P. In Retroviruses; Coffin, J. M., Hughes, S. H., Varmus, H. E., Eds.; Cold Spring Harbor Laboratory Press: Plainview, NY, 1997. (3) Safarianos, S. G.; Das, K.; Tantillo, C.; Clark, A. D.; Ding, J.; Whitcomb, J. M.; Boyer, P. L.; Hughes, S. H.; Arnold, E. Crystal Structure of HIV-1 Reverse Transcriptase in Complex with a Polypurine Tract RNA:DNA. EMBO J. 2001, 20, 1449−1461. (4) Molina, J. M. Efficacy and Safety of Once-Daily Regimens in the Treatment of HIV Infection. Drugs 2008, 68, 567−578. (5) Koczor, C. A.; Lewis, W. Nucleoside Reverse Transcriptase Inhibitor Toxicity and Mitochondrial DNA. Expert Opin. Drug Metab. Toxicol. 2010, 6, 1493−1504. (6) Desai, V. G.; Lee, T.; Delongchamp, R. R.; Leakey, J. E.; Lewis, S. M.; Lee, F.; Moland, C. L.; Branham, W. S.; Fuscoe, J. C. Nucleoside Reverse Transcriptase Inhibitors (NRTIs)-Induced 181-195. Expression Profile of Mitochondria-Related Genes in the Mouse Liver. Mitochondrion 2008, 8, 181−195. (7) Das, K.; Martinez, S. E.; Bauman, J. D.; Arnold, E. HIV-1 Reverse Transcriptase Complex with DNA and Nevirapine Reveals NonJ

dx.doi.org/10.1021/jp506119h | J. Phys. Chem. B XXXX, XXX, XXX−XXX

The Journal of Physical Chemistry B

Article

(25) Das, K.; Martinez, S. E.; Bauman, J. D.; Arnold, E. HIV-1 Reverse Transcriptase Complex with DNA and Nevirapine Reveals Non-nucleoside Inhibition Mechanism. Nat. Struct. Mol. Biol. 2012, 19, 253−259. (26) Ren, J.; Milton, J.; Weaver, K. L.; Short, S. A.; Stuart, D. I.; Stammers, D. K. Structural Basis for the Resilience of Efavirenz (DMP266) to Drug Resistance Mutations in HIV-1 Reverse Transcriptase. Struct. Fold. Des. 2000, 8, 1089−1094. (27) Das, K.; Bauman, J. D.; Clark, A. D.; Frenkel, Y. V.; Lewi, P. J.; Shatkin, A. J.; Hughes, S. H.; Arnold, E. High-Resolution Structures of HIV-1 Reverse Transcriptase/TMC278 Complexes: Strategic Flexibility Explains Potency Against Resistance Mutations. Proc. Natl. Acad. Sci. U.S.A. 2008, 105, 1466−1471. (28) Das, K.; et al. Roles of Conformational and Positional Adaptability in Structure-Based Design of TMC125-R165335 (Etravirine) and Related Non-nucleoside Reverse Transcriptase Inhibitors That Are Highly Potent and Effective against Wild-Type and Drug-Resistant HIV-1 Variants. J. Med. Chem. 2004, 47, 2550− 2560. (29) Esnouf, R. M.; Ren, J.; Hopkins, A. L.; Ross, C. K.; Jones, E. Y.; Stammers, D. K.; Stuart, D. I. Unique Features in the Structure of the Complex between HIV-1 Reverse Transcriptase and the Bis(heteroaryl)piperazine (BHAP) U-90152 Explain Resistance Mutations for this Nonnucleoside Inhibitor. Proc. Natl. Acad. Sci. U.S.A. 1997, 94, 3984−3989. (30) Li, H.; Robertson, A. D.; Jensen, J. H. Very Fast Empirical Prediction and Interpretation of Protein pKa Values. Proteins: Struct. Funct. Bioinf. 2005, 61, 704−721. (31) Bas, D. C.; Rogers, D. M.; Jensen, J. H. Very Fast Prediction and Rationalization of pKa Values for Protein-Ligand Complexes. Proteins: Struct. Funct. Bioinf. 2008, 73, 765−783. (32) Olsson, M. H. M.; Søndergard, C. R.; Rostkowski, M.; Jensen, J. H. PROPKA3: Consistent Treatment of Internal and Surface Residues in Empirical pKa predictions. J. Chem. Theory Comput. 2011, 7, 525− 537. (33) Phillips, J. C.; Braun, R.; Wang, W.; Gumbart, J.; Tajkhorshid, E.; Villa, E.; Chipot, C.; Skeel, R. D.; Kale, L.; Schulten, K. Scalable Molecular Dynamics with NAMD. J. Comput. Chem. 2005, 26, 1781− 1802. (34) Case, D. A.; et al. AMBER 9; University of California: San Francisco, 2006. (35) Dewar, M. J. S.; Zoebisch, E. G.; Healy, E. F.; Stewart, J. J. P. Development and Use of Quantum Mechanical Molecular Models. 76. AM1: a New General Purpose Quantum Mechanical Molecular Model. J. Am. Chem. Soc. 1985, 107, 3902−3909. (36) Weiner, S. J.; Kollman, P. A.; Nguyen, D. T.; Case, D. A. An All Atom Force Field for Simulations of Proteins and Nucleic Acids. J. Comput. Chem. 1986, 7, 230−252. (37) Weiner, S. J.; Kollman, P. A.; Case, D. A.; Singh, U. C.; Ghio, C.; Alagona, G.; Profeta, S., Jr.; Weiner, P. K. A New Force Field for Molecular Mechanical Simulation of Nucleic Acids and Proteins. J. Am. Chem. Soc. 1984, 106, 765−784. (38) Weiner, P. K.; Kollman, P. A. AMBER: Assisted Model Building with Energy Refinement. A General Program for Modeling Molecules and Their Interactions. J. Comput. Chem. 1981, 2, 287−303. (39) Jorgensen, W. L.; Chandrasekhar, J.; Madura, J. D.; Impey, R. W.; Klein, M. L. Comparison of Simple Potential Functions for Simulating Liquid Water. J. Chem. Phys. 1983, 79, 926−935. (40) Field, M. J. A Practical Introduction to the Simulation of Molecular Systems; Cambridge University Press: Cambridge, U.K., 1999. (41) Field, M. J.; Albe, M.; Bret, C.; Proust-de Martin, F.; Thomas, A. The Dynamo Library for Molecular Simulations Using Hybrid Quantum Mechanical and Molecular Mechanical Potentials. J. Comput. Chem. 2000, 21, 1088−1100. (42) Ruggiero, G. D.; Guy, S. J.; Martí, S.; Moliner, V.; Williams, I. H. Vibrational Analysis of the Chorismate Rearrangement: Relaxed Force Constants, Isotope Effects and Activation Entropies Calculated for Reaction in Vacuum, Water and the Active Site of Chorismate Mutase. J. Phys. Org. Chem. 2004, 17, 592−601.

(43) Williams, I. H. Kinetic Isotope Effects from QM/MM Subset Hessians: “Cut-Off” Analysis for SN2 Methyl Transfer in Solution. J. Chem. Theory Comput. 2012, 8, 542−553. (44) Cheng, Y.; Prusoff, W. H. Relationship between the Inhibition Constant (K1) and the Concentration of Inhibitor which causes 50% Inhibition (I50) of an Enzymatic Reaction. Biochem. Pharmacol. 1973, 22, 3099−3108. (45) Agneskog, E.; Nowak, P.; Källander, C. F. R.; Sönnerborg, A. Evaluation of Etravirine Resistance in Clinical Samples by a Simple Phenotypic Test. J. Med. Virol. 2013, 85, 703−708.

K

dx.doi.org/10.1021/jp506119h | J. Phys. Chem. B XXXX, XXX, XXX−XXX