What do docking and QSAR tell us about the design of HIV-1 reverse transcriptase nonnucleoside inhibitors?

Despite vigorous studies, effective nonnucleoside inhibitors of HIV-1 reverse transcriptase (NNRTIs) are still in demand, not only due to toxicity and detrimental side effects of currently used drugs but also because of the emergence of multidrug-resistant viral strains. In this contribution, we present results of docking of 47 inhibitors to 107 allosteric centers of HIV-1 reverse transcriptase. Based on the average binding scores, we have constructed QSAR equations to elucidate directions of further developments in the inhibitor design that come from this structural data. Electronic supplementary material The online version of this article (10.1007/s00894-017-3489-3) contains supplementary material, which is available to authorized users.


Introduction
According to the latest estimates from The Joint United Nations Programme on HIV and AIDS, there were 36.7 million people living with HIV in 2015, up from 33.3 million in 2010. Among them, 1.1 million died from AIDS-related illnesses [1,2]. The currently licensed antiviral medications include drugs falling into five main classes targeting different steps in the HIV life cycle: reverse transcription (nucleoside/ nucleotide and non-nucleoside reverse transcriptase inhibitors), assembly, budding, and maturation (protease inhibitors), viral entry (fusion inhibitors and co-receptor antagonists), and integration (integrase inhibitors). Standard first-line antiretroviral therapy (ART) for adults consists of three reverse transcriptase inhibitors (two nucleosides plus nonnucleoside) or two nucleoside reverse transcriptase inhibitors in combination with an integrase inhibitor [3,4]. Although ART has its limitations, it has proved largely effective, particularly in the early stages of the disease, reducing the mortality from AIDS and turning the disease from lethal to chronic [5]. Because of the rapid emergence of multidrug-resistant viral strains, toxicity, and detrimental side effects caused by long-term drug treatment [6], however, the discovery of new antiviral agents, in particular those targeting HIV-1 reverse transcriptase, is a research priority.
HIV-1 reverse transcriptase (RT) is an asymmetric heterodimeric enzyme composed of a catalytically active 66-kDa subunit with three domains: polymerase (residues 1-319), connection (residues 320-440) and RNase H (residues 441-560), and a catalytically inactive 51-kDa subunit containing non-functional polymerase and connection domains [7]. Nearly half of approved antiretroviral drugs target the polymerase active site of RT or the allosteric cavity located 10 Å away [8]. Five of them are nonnucleoside RT inhibitors that bind to the allosteric site, called the non-nucleoside inhibitor binding pocket (NNIBP), in a noncompetitive manner, resulting in limitation of conformational flexibility required for efficient DNA synthesis by RT [9]. Nevirapine (NVP), delavirdine (DLV), and efavirenz (EFV) are first-generation non-nucleoside reverse transcriptase inhibitors (NNRTIs) with This paper belongs to Topical Collection P. Politzer 80th Birthday Festschrift Electronic supplementary material The online version of this article (https://doi.org/10.1007/s00894-017-3489-3) contains supplementary material, which is available to authorized users. potent anti-HIV-1 activity, but their adverse effects, particularly those affecting the central nervous system and  hepatotoxicity, poor resistance profile, and low genetic  barriers for viral resistance, especially single mutants  K103N, Y181C, and double mutant K103N/Y181C, lim-ited their clinical application [10]. The second-generation NNRTIs etravirine (ETR) and rilpivirine (RPV) have been approved by the Food and Drug Administration and European Union in 2008 and 2011, respectively. Both display a broad spectrum of activity against clinically relevant HIV-1 mutant and wild-type (wt) HIV-1 strains at low nanomolar concentrations. However, etravirine must be given twice daily and is not recommended for initial treatment of HIV infection because of insufficient data in treatment-naive individuals whilst rilpivirine has suboptimal efficacy in patients with viral loads greater than 100,000 copies/ml or CD4 cell counts below 200 copies/ml at baseline because of higher virologic failure rate. Moreover, their use in some settings is limited by hypersensitivity reactions or other adverse effects. In addition, although these drugs have higher genetic barriers to resistance than first-generation NNRTIs, drug resistance still emerges, and K103N and E138K are the most common m utations selected by ETV and RPV. Importantly, when resistance to RPV is selected after virologic failure, cross-resistance to all nonnucleoside reverse transcriptase inhibitors is commonly observed [11][12][13][14][15][16][17][18][19][20][21]. Therefore, there is a considerable need for new drugs with anti-HIV-1 reverse transcriptase activity. To address this problem, in this contribution we present results of docking of 47 inhibitors bound to 107 allosteric centers of HIV-1 reverse transcriptase, i.e., all data available in the Protein Data Bank. Based on the average binding scores, we have constructed QSAR equations to elucidate directions of further developments in the inhibitor design that come from this structural data.
In the above 107 structures, 48 ligands are present. They represent a number of different classes of chemicals. By far the most frequently studied (25 entries) ligand is NVP. Codes and structures of the remaining ligands are collected in Table 1. Out of compounds currently in clinical use, only DLV is not present in the set, as it was crystallized in the pocket of HIV-2 reverse transcriptase, which differs structurally from the allosteric pocket of HIV-1 enzyme [23]. Furthermore, we have excluded from the studies QO9 ligand, as it contains silicon atoms for which there is no parametrization. These ligands have been used in docking and subsequently in QSAR analysis.
Docking was performed using the FlexX program [24,25], as implemented in the LeadIT software package [26]. The receptor was prepared using graphical interface of the package. Both protein chains were selected, and the binding site was defined to include residues within a 6.5-Å radius around the native ligand. The library of ligands was imported from the .mol2 file. Protonation states corresponding to aqueous solution were used. Soft docking (allowing for volume overlap up to 100 Å 3 ) was performed. For ligand base placement, the default hybrid Enthalpy and Entropy strategy was used. The Clash Factor was set to 0.6. Other parameters were kept at default. Scores of the top-ranked poses of each ligand were averaged for a given enzyme structure.
QSAR analysis was performed in two different modes. In both cases, the training set included 47 structures; the endpoint was the average docking score obtained from FlexX calculations after exclusion of results for QO9, which as the only one yielded positive values of the score function. The size of the data set was too small to allow splitting it into training and validation sets, so the internal leave-one-out cross-validation approach was chosen instead.
In the first approach, based on classical descriptors, several models were created using the Complete Topological QSAR method and regression equation, created by feature selection with Enhanced Replacement Method [ERM] as implemented in SCIGRESS Suite software [27]. In the second approach, QSAR [28] was based on molecular fragments contribution. The common substructures were extracted from a training set, yielding a set of 96 substructure-count descriptors by means of an algorithm based on RECAP using ADMEWORKS ModelBuilder [29,30]. The substructures, present in less than three of training structures, were removed by means of zero-test with a threshold of 6%, leaving 39 substructure count descriptors. Particle Swarm Optimization algorithm [31] was employed for feature selection with a target of selecting 15 descriptors. After approximately 7000 iterations of 10,000 model population, the process was manually interrupted. Eighteen of the most often used descriptors were selected. The final model was created using leaps-and-bounds multiple linear regression model, a variation of backward stepwise regression.

Results and discussion
All 47 ligands present in the PDB that are bound in the allosteric cavity have been docked to all 107 structures and averaged scores for a given ligand were obtained separately for wild-type (wt) and for mutated enzyme (individual data provided in Table S1 in Supporting information). The obtained poses have been inspected for correct orientation within the allosteric cavity (for examples of overlap with the native ligand see Figs. S1 and S2 in the Supporting Information). The results are collected in Table 2. Averaged binding scores have been compared for wt and mutated enzymes. The results are illustrated graphically by Fig. 1. The strong linear correlation obtained indicates that there is no significant difference between binding in either form of the enzyme. Furthermore, as illustrated by Fig. 2, a slight preference for binding in the allosteric pocket of either wt enzyme or its mutated form is random and does not correlate with the energy of binding. The difference is symmetrically distributed between positive and negative values showing practically no systematic preference of binding to either wild-type or one of the mutated forms of the enzyme. Similarly, we have found no correlation between the standard deviation of the average binding score and the binding energy. This observation indicates that activity against mutated HIV-1 RT forms is not governed by the strength of binding. Allosteric ligands impair enzyme action by a wedge mechanism, hindering domain mobility toward opening and closing the access to the active site. However, final allosteric site architecture is achieved upon ligand binding. In order to account for this flexibility and possible clash between the protein and a ligand, we have used large overlap volume (100 Å 3 ). Lack of systematic difference between binding to wt and mutated enzyme seems thus to indicate that activity against mutants is connected with the structural features of the ligand rather than their binding energy. Interactions within the allosteric site are mostly associated with van der Waals forces and to a lesser extend to hydrogen bonding [32]. As illustrated by the most suited for mutant enzymes ligand, EFZ, its success seems to come from hydrogen bonding to lysine 101 rather than lysine 103, which is the most frequent mutation (see left panel of Fig. S1). Based on the results obtained from docking, we have performed QSAR studies. In order to understand the very general structural features affecting the binding affinity, we attempted two approaches to QSAR modeling, aiming for models containing descriptors that have clear "chemical interpretation" and can be easily interpreted to give suggestions for new compound design. In the first approach, we have used "traditional" QSAR that involves a library of descriptors as implemented in the SCIGRESS Suite software.
Several models gave acceptable correlation coefficient, while when only easily interpretable descriptors were used the following equation turned out to be statistically the best: FlexX score ¼ −2:035 * double bond count−3:268 * sec−amine count−2:570 * cyano countþ 6:83817e−05 * nonpolar area 2 −16:375 ð1Þ with r 2 equal to 0.8358; the degrees-of-freedom adjusted r 2 equal to 0.8202; and the leave-one-out cross-validated r 2 of 0.7819. The standard deviation in the error predicted by leave-one-out cross-validation is 2.5150. The F-ratio is 53.4646. The probability that a greater F-ratio can be obtained by chance alone is 0.0000. Thus, since the probability is less than 0.05, there is at least one significant descriptor in the model. Partial F values (listed below) indicate relative importance of each descriptorthe larger the value the more important is the descriptor. Topological results are illustrated by Fig. 3.
The ranges of values for the descriptors in the training set and obtained corresponding partial-F values (in parenthesis) were: Since the objective is to have compounds with the lowest (most negative) FlexX score, the model given by Eq. (1) suggests that molecules should contain nitrile and secondary amine groups, and the area of the molecule incapable of hydrogen bonding (either as a donor or an acceptor) should be as small as possible.
The second attempt aimed at creating QSAR using fragment contribution approach using common substructures present in the training set using ADMEWORKS ModelBuilder. Due to size of the training set, the set of six descriptors was chosen. As illustrated by Fig. 4, this is the lowest number of descriptors that yields acceptable statistically significant results. The set contained X-H (hydrogen attached to any atom) substructure count descriptor. For simpler mechanistic interpretation, the descriptor was manually replaced with C-H count (hydrogens attached to carbon) to calculate the final model. The obtained results are presented in Fig. 5, while the final statistical parameters of this model are collected in Table 3.
The model cross-validation r 2 of less than 70% does not encourage its use for direct prediction of unknown compounds. However, the sign of the linear regression equation's weight vector coefficients is a measure of the influence of a given substructure contribution to activity. In particular, negative values indicate improvement in binding, while positive values suggest that the corresponding substructures should be eliminated or their presence minimized. The obtained results are summarized in Table 4, where substructures with positive contribution to binding are presented in bold, while those which should be avoided are distinguished by italics.
As can be seen, the interpretation of results of the fragmentbased QSAR leads to similar conclusions, as the interpretation of the topological QSARthe presence of sec-amino and cyano groups is improving activity, while parts of purely hydrocarbon-nature should be avoided.

Conclusions
Two main conclusions are drawn from the research presented in this study: 1. HIV-1 reverse transcriptase (RT) has relatively low fidelityits error rate is evaluated to be in the range of 10 −5 per nucleotide addition [33,34]  main source of the activity toward mutants rather than the binding energy, 2. The interpretation of results of fragment-based QSAR leads to similar conclusions as interpretation of the topological QSAR: the presence of sec-amino and cyano groups is improving binding, while hydrocarbon fragments should be avoided.