Theoretical studies of energetics and binding isotope effects of binding a triazole-based inhibitor to HIV-1 reverse transcriptase

Understanding of protein-ligand interactions is crucial for rational drug design. Binding isotope effects, BIEs, can provide intimate details of specific interactions between individual atoms of an inhibitor and the binding pocket. We have applied multi-scale QM/MM simulations to evaluate binding energetics of a novel triazole-based non-nucleoside inhibitor of HIV-1 reverse transcriptase and to calculate associated BIEs. The binding sites can be distinguished based on the O-BIE.


Introduction
Human immunodeficiency virus (HIV-1) is a cytopathic human retrovirus that primarily infects CD+ helper T lymphocytes leading to the loss of their number and function.The progressive failure of the immune system associated with the loss of predominant T lymphocytes and other cells critical to this system, e.g.macrophages or dendritic cells, leads to acquired immune deficiency syndrome (AIDS).The progressive disappearance of immunity weakens the body, allows opportunistic infections to evolve, leads to cancers and ultimately to premature death. 1 Nowadays, treatment of HIV-1 infected patients relies on cocktails of different drugs.Due to undesired side-effects, especially occurring during long-term treatments, there has been considerable interest in developing new inhibitors targeting three HIV-1 enzymes: integrase (IN), protease (PR) and reverse transcriptase (RT). 2 Herein we focus on an inhibitor of the latter enzyme.
HIV-1 RT is an asymmetric heterodimer, composed of p66 and p51 subunits, which catalyses the transcription of singlestranded viral RNA into viral DNA double helices. 3The enzyme contains one DNA polymerase active site and one ribonuclease H (RNase H) active site, both of which reside in the p66 subunit in spatially distinct regions while the smaller p51 chain is responsible for structural scaffolding.The hydrophilic polymerase active site is located between palm (light blue), finger (blue) and thumb (cyan) subdomains (see Fig. 1) and is made up of Asp110, Asp185 and Asp186 amino acids that can coordinate two Mg 2+ or two Mn 2+ cations. 4The allosteric cavity, composed of Lys103, Tyr181, Tyr188, Phe227, Trp229, Leu234 and Tyr318 residues 5 (highlighted in Fig. 1), lies approximately 10 Å from the polymerase active site.The RNase H active site comprising Asp443, Asp498 and Asp549 amino acids and two bound Mg 2+ cations 6 is located 60 Å from the allosteric pocket.The two active sites and the allosteric cavity constitute a key target of antiviral research.Nucleoside reverse transcriptase inhibitors (NRTIs) inhibit the polymerase active site and need to be phosphorylated by cellular kinases to active triphosphate form, which competes with natural dNTPs and finally becomes incorporated into the growing primer.Nucleotide reverse transcriptase inhibitors (NtRTIs), in contrast to NRTIs, are structures already equipped with the phosphonate group and after the phosphorylation step are able to terminate DNA synthesis also by the competitive mechanism. 7Non-nucleoside RT inhibitors (NNRTIs), on the other hand, follow the non-competitive mechanism of inhibition.They are considered to be highly specific compounds that bind tightly in the allosteric pocket, acting on the enzyme subunits as a molecular wedge. 8,9Simultaneous inhibition by combining nucleoside and non-nucleoside inhibitors is also currently explored. 10,11To date eight NRTIs, two NtRTIs and five NNRTIs have been approved by FDA. 12 In the present study we focus on a ligand that can bind either in the allosteric cavity or the RNase H active site of HIV-1 RT and compare it to ligands designed specifically for one of these sites (structural formulas of all ligands are presented in ESI, † Fig. S1).
An accurate description of binding affinity of an inhibitor to an enzyme is still a challenging task. 13,146][17][18] Binding energies of NNRTIs have also been studied by MM-PBSA 19 and quantum calculations 20 including the ONIOM method. 21Warshel and co-workers obtained the absolute binding free energies of NNRTIs by the PDLD/S-LRA/b method. 22Our own experience, first validated on the FDA approved drugs, indicates that binding free energies of HIV-1 RT complexes with NNRTIs and NRTIs can be efficiently obtained using a thermodynamic cycle from alchemical free energy perturbation (FEP) calculations. 23,24ecently, Schramm and co-workers showed that binding isotope effects, BIEs, 25 exhibit different values upon binding an inhibitor, para-aminobenzoic acid, to an enzyme (dihydropteroate synthase) from different sources (Staphylococcus aureus and Plasmodium falciparum).In this communication we extend these findings to the identification of the binding site within a single enzyme using BIEs on the example of a novel potent HIV-1 RT triazole-based inhibitor.

System setup
Novel triazole derivative N-(2-chloro-4-sulphamoylphenyl)-2-((4-(2,4-dimethyl-phenyl)-5-(thiophen-2-yl)-4H-1,2,4-triazol-3yl)sulphanyl)-acetamide, L-1, 26 which was found to be a potent NNRTI, has been docked to the allosteric cavity and the RNase H active site of HIV-1 RT based on the crystallographic structure from the Protein Data Bank with PDB ID 2RKI 27 using Glide package as implemented in the Schro ¨dinger program. 28Both complexes are depicted in Fig. 1.The 2RKI structure was crystalized with 4-benzyl-3-[(2-chlorobenzyl)sulfanyl]-5-thiophen-2-yl-4H-1,2,4-triazole mole-molecule bound in the allosteric cavity.1][32] Assuming a physiological pH value, hydrogen atoms were added using the tLEAP module of AMBER program. 33eutralization of the L-1-allosteric cavity and L-1-RNase H active site complexes was completed by adding 15 and 9 Cl À counterions, respectively.Finally, both complexes were placed in orthorhombic boxes of TIP3P 34 water molecules of 144 Â 160 Â 144 Å 3 .Parameters for L-1 were obtained from GAFF 35 as implemented in AMBER Tools.Afterwards, several optimization and dynamics simulations were performed using AMBER ff03 36 implemented in NAMD 37 program with a time step of 1 fs and periodic boundary conditions using the particle mesh Ewald method. 38The cut-off for nonbonding interactions was applied using a smooth switching function with a radius range from 14.5 to 16 Å.First energy optimizations were carried out by means of the conjugated gradient algorithm.The systems were heated by increasing temperature from 0 to 300 K with 0.001 K increment.The equilibration of the systems was achieved during 300 ps of Langevin-Verlet dynamics at 300 K. Then 2 ns of NVT MM MD were carried out.Subsequently, 200 ps for the L-1-cavity and 1 ns for L-1-RNase H complexes were simulated using QM/MM MD simulations at the AM1 39 /AMBER:TIP3P theory level as implemented in fDynamo. 40All atoms beyond 20 Å from L-1 were kept frozen.In order to describe the binding energy of L-1 to HIV-1 RT, a reference system of L-1 in an aqueous solution of the equivalent size (144 Â 160 Â 144 Å 3 ) was simulated during 200 ps of AM1/TIP3P dynamics.

Alchemical free energy perturbation (FEP) method
The enzyme-ligand binding affinity is a main goal of drug design.Among the different methods developed for this purpose, the alchemical free energy perturbation (FEP) combined with QM/MM molecular dynamics, while not providing detailed information about the full binding process, delivers reliable values of binding free energy, DDG bind , at a reasonable computational cost.The FEP method by means of the thermodynamic cycle was first introduced by Warshel and co-workers. 41Herein, a series of QM/MM MD have been carried out at 300 K in the NVT ensemble with two parameters l and g, which correspond to electrostatic and van der Waals interactions, respectively: Both l and g were smoothly changed from 0 to 1 with 0.02 increments.According to our own experience with HIV-1 RT FEP calculations in each window 5 ps of relaxation followed by 100 ps of QM/MM MD need to be carried out and a simplified thermodynamic cycle, depicted in Fig. 2, can be used to obtain DDG bind based on eqn (2).
The potential energy

In order to describe individual interactions between allosteric cavity amino acids and the triazole molecule, the analysis of
This journal is © the Owner Societies 2016 averaged potential energy has been performed through QM/ MM MD based on eqn (3) where E INT QM=MM is composed of electrostatic, E elect QM=MM , and van der Waals, E vdw QM=MM , interaction terms as a function of interaction energy between QM and MM parts.The polarized QM subsystem includes the coulombic interaction of the QM nuclei (Z QM ) and the electrostatic interaction of the polarized electronic wave-function (C) with enzyme amino acid charges (q MM ), while the MM part brings non-polarizable potential contribution.

Molecular electrostatic potential (MEP)
The molecular electrostatic potential, MEP, describes the interaction energy between the charge distribution of a molecule and a unit positive charge.This method allows predicting nuclear and electronic charge distribution of a given molecule.Despite the fact that polarization is not included, charge distribution remains unperturbed by the external test charge, and the MEP is still effective for the interpretation and prediction of chemical reactivity. 42Three-dimensional MEP surfaces for L-1 extracted from the allosteric cavity and the RNase H active site were obtained at the B3LYP/6-31+G(d,p) theory level using Gaussian09. 43Negative electronic potential corresponds to the most nucleophilic regions (red), while positive electronic potential indicates the most electrophilic regions (blue).

Binding isotope effects (BIEs)
Binding isotope effects can provide information on conformational changes and ligand-enzyme interactions. 44In combination with the experimental values they can also be used to determine the actual binding site.We have already successfully used BIEs to analyse interactions between ligands and HIV-1 RT. 24 BIEs can be expressed by the binding free energy of the light, L, and heavy, H, isotopologs considering the relationship between the binding free energy and the equilibrium constant: where the total partition function Q e , for the ligand in the enzyme-ligand complex, and Q w , for the ligand in the ligandwater complex, are computed as products of the translational, rotational and vibrational partition functions.The different contributions were separately computed using the Born-Oppenheimer, rigid rotor and harmonic oscillator approximations.Due to the fact that in each case the ligand is in a condensed phase, rotational and translational motions are replaced by six librational motions of the ligand within its external environment.However, we decided to subject the full 3N Â 3N Hessians for the ligand to a projection procedure to eliminate translational and rotational components, which gives rise to small non-zero frequencies, as previously described. 44,45Thus, 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 substrate (in contrast to ref.

Results and discussion
The obtained total free energies of binding inhibitors to HIV-1 RT binding sites are collected in Table 1.A more extensive comparison is presented in Fig. 3.There are some difficulties in the validation of DDG bind using experimental data.Our experience shows that IC 50 measurements are not suitable for this purpose 23,24 in contrast to the reciprocal of the dissociation constant (K I ) of the protein-ligand complex. 47Unfortunately, no systematic study of K I values of different NNRTIs measured under the same experimental conditions can be found in the literature.The evaluation of our FEP calculations was made by using NRTI 23 and NNRTI 24   In this case the magnitude of the total energy of binding of L-1 to the allosteric cavity is thus solely due to the van der Waals term (see the green bars in Fig. 3), which is preferable in the hydrophobic pocket.Usually, the van der Waals term is much smaller than the electrostatic term, hence the relative DDG vdW value obtained for L-1 is the highest among NNRTIs.This makes L-1 one of the most potent NNRTIs.Binding of L-1 to the RNase H active site is due to preferable electrostatic forces with DDG bind equal to À36.5 kcal mol À1 .There is a meaningful difference between L-1 interactions in aqueous solution and the enzyme environment in both electrostatic and van der Waals terms.NRTIs show the same behaviour in the case of electrostatic interactions, but opposite in van der Waals interactions.Consequently, binding L-1 to the RNase H site is accompanied by both strong electrostatic and van der Waals terms, and indicates that L-1 can successfully bind into the RNase H active site.The potential energy of electrostatic and van der Waals interactions between L-1 and allosteric cavity amino acids has been obtained along the QM/MM dynamics; the averaged results are displayed in Fig. 4. Negative values correspond to favourable interactions and stabilize the enzyme-ligand complex, while positive values represent unfavourable interactions.The strength of interactions depends Table 1 The total free energy of binding (kcal mol À1 ), DDG bind , inhibitors to HIV-1 RT binding sites obtained from the FEP method.DG elect and DG vdW components correspond to the electrostatic and van der Waals terms, respectively, aq indicates the ligand in aqueous solution while HIV in an enzymatic environment Fig. 3 Contributions of electrostatic and van der Waals terms to total binding free energy obtained from FEP calculations.

Inhibitor
Fig. 4 Analysis of potential energy of interactions between L-1 and allosteric cavity amino acids.
on the nature of the ligand (size, stereochemistry, number and type of functional groups) and differs among NNRTIs.The crystallographic structures 5,48 indicate that Leu100, Lys101, Lys103, Val106, Thr107, Val108, Val179, Tyr181, Tyr188, Val189, Gly190, Phe227, Trp229, Leu234, and Tyr318 all from chain A and Glu138 from chain B are usually responsible for inhibitorenzyme interactions.The L-1-allosteric cavity complex has meaningful interactions with Leu100, Lys103, Lys104, Val106 and Pro236 as shown in Fig. 4. L-1 interacts also with Phe227, His235 and Tyr318.Residues Lys101, Thr107, Val108, Val189, Gly190, Trp229, Leu234, and Tyr318 are too far, while Val179, Tyr181, Tyr188, and Glu138 from chain B exercise negligible bending affinity toward L-1.The main contribution of favourable electrostatic interactions, shown in Fig. 4 and 5, is observed between 1,2-triazole nitrogen atoms and Lys103.However, these interactions do not enhance binding affinities, because interactions of L-1 with hydrogen atoms of solvent molecules fully compensate them.Significantly weaker are favourable electrostatic interactions with Lys104 and Tyr318, wherein carbonyl oxygen from Lys104 interaction with the hydrogen atom from the aromatic ring is further stabilized by the p-stacking with Tyr318.These positive electrostatic interactions are compensated by the electrostatic repulsion of Asp192.The determinant factor of L-1-RT attraction comes from the van der Waals interactions, not present in water, especially with Val106, Leu100, Pro236, and His235.Residues Tyr181 and Tyr188 are stabilizing one aromatic ring of L-1 with two methyl groups and Phe227 interacts by p-stacking with another aromatic ring neighbouring the sulphate group.The electrostatic positive and negative terms are equal, therefore the binding process is a result of van der Waals interactions between L-1 and the allosteric cavity; the potential energy of van der Waals interactions is sufficient for incorporating and maintaining L-1 in this hydrophobic pocket.
In contrast to the L-1-allosteric cavity complex, the electrostatic interactions between L-1 and RNase H active sites are crucial to hold the ligand in the active site centre.The RNase H active site is made up of Asp443, Glu478, Asp498 and Asp549 amino acids, which coordinate two Mg 2+ .The magnesium cation tends to create six coordination bonds but in the active site usually only five are occupied.Hence, two cations can interact with ligands by two relatively strong coordination bonds (Fig. 6).
Small ligands, like L NA -1, 23 interact with both Mg 2+ without reorganizing active sites.An opposite situation is observed for bigger molecules such as L NA -4 and L-1.Both alter the RNase H active site in the same way, causing an increase of the distance between magnesium cations and forcing coordination of Val552 and Ser553 by one Mg 2+ .The open and flexible character of the RNase H region supports the reorganization of the active site.This dynamic effect, observed along MM MD, supports binding of a bigger ligand such as L-1.
Recently, effective medical treatment of HIV infection has been focused on the composite approach combining different targets (enzymes) of the virus as well as different sites within a single enzyme.It is assumed that inhibitors of different types bind to different sites of the enzyme.However, recent reports on finding several alternative binding sites within HIV-1 RT shed new light on the actual binding sites of particular inhibitors. 49In order to analyse the specificity of the triazole-based inhibitors we have performed calculations of the binding energetics of L-1 in two binding spots: the allosteric centre and the RNase H active site.After QM/MM MD simulations we have extracted the last structure of L-1 from the L-1-cavity and L-1-RNaseH complexes, and calculated three-dimensional MEP surfaces.The obtained MEP surfaces are shown in Fig. 7, where the most nucleophilic regions are in red (À0.1 a.u.) and the most electrophilic regions are in blue (+0.1 a.u.).Both structures have similar distribution of nuclear and electronic charges.As expected, the triazole ring and the sulphate group exhibit negative electronic potential.The sulphate group is dangling almost outside of the hydrophobic pocket, thus it can interact with water molecules.In the L-1-RNase H complex, the active site is open to the aqueous environment and also in this case the sulphate group can interact with solvent molecules.The L-1 nitrogen atoms are deep inside the allosteric cavity and in the presence of Lys103 they act as nucleophiles, while in the RNase H active site they interact electrostatically with Asp498, Gln500 and Tyr501.Consequently, there is no big difference between the charge distribution of L-1 extracted from an allosteric cavity and from the RNase H active site.These results indicate that L-1 can bind to both binding sites.Another analysis of interactions of L-1 in complexes with the allosteric cavity and the RNase H active site has been based on heavy atom ( 13 C, 15 N, 18 O, 37 Cl and 34 S) BIEs.The calculated values are collected in Fig. 7. BIEs can provide information on differences in interactions between L-1 dissolved in water and bound in RT binding sites.Despite the dissimilarity in the molecular character of the hydrophobic cavity and the hydrophilic RNase H active site, in both RT binding sites lack of BIEs for most of the atoms has been observed.The largest isotope effects were obtained for 18 O of the sulphate group and apparently this isotope effect can be experimentally used to distinguish the place in which L-1 is bound; a normal 18 O-BIE (1.016 AE 0.002) is expected when L-1 binds in the allosteric cavity, where oxygen atom interacts (see Fig. 4) with the hydrogen atom from Val106 (the distance is 4.23 Å) and is also surrounded by three water molecules.On the other hand an inverse BIE (0.988 AE 0.002) for sulphate group 18 O atoms is expected for the complex of L-1 with RNase H, originating from strong interactions with one of the active site Mg 2+ cations (the distance is 2.41 Å) and additionally from stabilization by the hydrogen atom (the distance is also 2.41 Å) from the methyl group of Ala445 and one water molecule.Isotope effect differentiation of the other oxygen of the sulphonyl group is not possible as it is only slightly inverse (0.997 AE 0.001) when L-1 is bound in the RNase H active site in contrast to the lack of the isotope effect of this oxygen atom in the allosteric cavity.This slightly inverse value of BIE originates from the interaction with the second Mg 2+ cation (the distance is 2.56 Å).Although electrostatic interaction between the nitrogen atom of the triazole ring and the hydrogen atom from Lys103 (the distance is 2.52 Å) is unfavourable in the allosteric cavity, while in the RNase H active site electrostatic interactions of this atom with Gln500 and Tyr501 hydrogen atoms tend to stabilize the complex (the distances are 2.55 Å and 2.10 Å, respectively), there is no reflection of these differences in the isotopic binding patterns with both 15 N-BIEs being slightly inverse (0.997 AE 0.001) and equal.Finally, a negligible normal BIE (1.002 AE 0.002) was calculated for the 37 Cl-BIE when L-1 binds in the RNase H active site, and a similar value was obtained when L-1 is complexed in the allosteric cavity.Although it should be kept in mind that equilibrium isotope effects, such as BIEs, are generally small and hard to measure experimentally, the above analysis illustrates that it is possible to distinguish binding sites based on their values in favourable cases.

Conclusions
The novel inhibitor of HIV-1 RT, L-1, has been proposed and has been examined based on theoretical simulations.The L-1 molecule was docked to both RT binding sites: allosteric cavity and RNase H active site.Final complexes of the L-1-binding site have been obtained after long QM/MM dynamics.Once the stability of both systems was confirmed, the interactions between L-1 and RT binding sites have been quantified by free energy of binding computed using the FEP method.The analysis of this affinity has been also enriched by potential energy calculations, three-dimensional MEP isosurfaces, structural analysis and heavy atom ( 13 C, 15 N, 18 O, 37 Cl and 34 S) BIEs.The obtained results confirmed the differences between the allosteric cavity and the RNase H binding sites, as well as dissimilarity in the ligand binding pathways.In the allosteric cavity the inhibitor binds by van der Waals affinities, which are preferable in this hydrophobic pocket, while the electrostatic interactions are unfavourable.The opposite situation was observed for binding the ligand to the RNase H active site, where binding occurs by electrostatic forces.
Inhibition is accompanied by the release of standard free energy, DDG bind .The values of DDG bind were obtained by the alchemical FEP method.Comparison of different NNRTIs, and NRTIs, L-1 bound in the allosteric cavity and the RNase H active site shows that L-1 has the largest negative standard free energy of binding, indicating that L-1 is able to successfully inhibit both the allosteric cavity and RNase H sites. Also the three-dimensional MEP isosurface display shows that there is no large difference between L-1 charge distribution in the allosteric cavity and the RNase H active site.
BIE values are almost identical for both L-1-enzyme complexes, with one exception.In contrast to the L-1-allosteric cavity complex, where a normal (larger than unity) 18 O-BIE is expected, in the L-1-RNase H active site an inverse BIE for two oxygens from the sulphate group is predicted.This occurs due to the strong interactions with two Mg 2+ originally belonging to the RNase H active site.Thus BIEs allow us to distinguish to which site L-1 binds.Experimental determination of this 18 O-BIE is being currently carried out in our laboratory.
Science and Higher Education.The authors also acknowledge the Servei d'Informatica, Universitat Jaume I for generous allotment of the computer time.

Fig. 1
Fig. 1 Representation of HIV-1 RT with two binding sites, allosteric cavity (upper right box) and RNase H active site (lower right box) with a bound L-1 inhibitor.
46).Although the environment (protein or water) surrounding the ligand affects the values of each Hessian, there is no coupling between the ligand and its environment.Eleven snapshots from last 20 ps of QM/MM MD calculations were extracted for L-1-water, the L-1-allosteric cavity and L-1-RNase H active site complexes.The combination of 121 BIE individual values obtained with eleven enzyme-ligand and eleven ligand-water (11 Â 11) structures allowed for the calculations of BIEs with their uncertainties.All formalisms used for BIE calculations yielded practically the same values (see TablesS1 and S2in the ESI †).

Fig. 2
Fig.2Thermodynamic cycle to compute enzyme-ligand binding free energies from the FEP method, where (E) is the enzyme with ligand (L) in its binding site, E 0 is the apo form of enzyme, and (L) 0 is the ligand in the gas phase.

Fig. 5
Fig. 5 Arrangement of L-1 at the allosteric cavity with the average of key distances in Å.

Fig. 6
Fig. 6 Arrangement of L-1 at the RNase H active site with average key distances in Å. Magnesium cations are highlighted as green balls.