Spatial Arrangement of the Potential anti-HIV Drug, 1-[2-(2-Benzoylphenoxy)ethyl]-6-methyluracil, in Crystal and Ligand-Receptor Complex

1. Introduction

Atomistic simulations of binding between anti-HIV agents and various biological systems containing complex biomacromolecules supplies a theoretical basis for the investigation of the microscopic level processes during the antiretroviral therapy. X-Ray crystallographic data about ligand–receptor complexes between proteins or nucleic acids receptors and drug molecules often become a starting model for these simulations. Alternatively, the starting models can be obtained from molecular docking of small molecules of drugs and the receptor molecules which requires a high-resolution 3D structure of the receptor. We believe that descriptors based only on analysis of structural data from experimental or computational studies cannot supply a detail description of ligand-receptor binding. The methods that use electron density distribution function or wavefunction are much more informative for the aforementioned purpose. Indeed, a number of studies of such ligand–receptor complexes in terms of methods employing electron density distribution function [1,2], NCI [3,4] and source function [5,6] were published recently. As a result, the detailed picture of ligand-receptor binding was recovered including nature and energy of all possible ligand – receptor interatomic interactions.

In the present article we discuss the intermolecular interactions of 1-[2-(2-benzoylphenoxy)ethyl]-6-methyluracil (VMU-2012-05, Scheme 1) that was patented in Russia as new prospective nonnucleoside inhibitor of reverse transcriptase in its single crystal and complexes with protein and nucleobases to evaluate the conformational stability of VMU-2012-05 and its ability to bind with biomolecules (Patent RU2648998C1 [7]). According to published preliminary studies, VMU-2012-05 is characterized by rather low toxicity on rats [8] and demonstrated high activity for wild and mutant HIV-1 that 2.5 times exceeds that for nevirapine [7]. Thus, our study will be essential for a comprehensive understanding of the effect of VMU-2012-05 on molecular scale. The electron density distribution, the type and energies of intermolecular interaction for crystalline VMU-2012-05 could be obtained experimentally as a reference. The quantum crystallography methods were applied to molecular models of varying complexity from dimers of VMU-2012-05 with nucleosides to a ligand-receptor complex as the reliable computational scheme to assess the role of the intermolecular interactions using the theoretical electron density distribution function. The study of the ligand-protein receptor binding can be carried out with QM/MM calculations, but the level of information value is in direct relation with the size of QM subsystem. The description of ligand to DNA/RNA binding is more complex task due to lack of accurate crystallographic data. In our previous study the model with ~ 500 atoms was used for the study of binding of ferrocene derivative to the fragment of DNA [9] but the size of the model QM system required for correct description of DNA/RNA binding is unclear yet.

In order to evaluate effects of intermolecular interactions on molecular conformation of a flexible molecule in terms of energy, we performed conformational analysis of an isolated molecule of VMU-2015-05, compared stable conformations of an isolated molecule with that experimentally observed in crystal, and those in a model drug–nucleoside or a drug–aminoacid dimers as well as the ligand–protein complex. Herein we demonstrated that the main features of the stacking interactions and the H-bonding between the molecule of anti-HIV drugs and DNA/RNA can be sufficiently described using rather small models of stacking dimers. Vice versa, dimers of the drug and an aminoacid do not represent the ligand-protein binding; exhaustive study of binding situation with the methods based on electron density function requires all aminoacids surrounding the molecule of drug to be a part of the QM subsystem. The applied methodology paves the way to estimate the ability of a drug molecule to bind with biomolecules provided that the 3D structure of the receptor is unknown.

2. Materials and Methods

2.1. Dataset Collection, Reduction and Refinement

Single-crystals were grown by slow evaporation from a water-ethanol mixture at room temperature. A good quality single-crystal was selected from the precipitate and mounted on a MiTeGen® MiroMount® plastic loop. X-ray diffraction data were collected at 100(2) K on the Agilent Super Nova diffractometer equipped with the Oxford Cryostream cooling device and the Supernova microfocus X-ray source with the molybdenum anode (λ = 0.71073 Å). Omega and phi scans were carried out at several detector positions utilizing various exposure times to reach completeness and maintain sufficient redundancy at high diffraction angles. The measured intensities were integrated and corrected for absorption with the Rigaku OD CrysAlis Pro software [10].

2.2. Independent Atom Refinement

Crystal structure was solved using the SHELXT [11] software and refined against F²_hkl in anisotropic approximation for non-hydrogen atoms with SHELXL [12] and OLEX2 [13] programs. Riding hydrogens were used to calculate and refine hydrogen atomic positions. The obtained structure is shown in Figure 1.

2.3. Multipole Refinement

The charge distribution in VMU-2012-05 was obtained using multipole formalism [14] as implemented in the MoPro package [15] with core and valence electron density derived from wave functions fitted to a relativistic Dirac-Fock solution. In the first step, the scale factor was refined from all data. After that, the high-order refinement (sin(θ) / λ > 0.7 Å) of atomic positions and atomic displacement parameters incorporating all non-hydrogen atoms was employed, followed by refinement of hydrogen atoms positions, with C–H, N–H, and O–H distances fixed at values obtained from neutron diffraction [16]. The ADPs for hydrogen atoms were estimated using SHADE3 software [17]. All non-hydrogen atoms were described at the octupolar level. For the H atoms, only the monopole population and the dipole populations were refined, while κ and κ’ were fixed to H atom theoretical values [18]. Separate κ parameters were used for all non-hydrogen atoms. During the final stage of the refinement, the multipole parameters, positions, and thermal parameters of all non-hydrogen atoms and monopoles were refined. All bonded pairs of atoms satisfied the Hirshfeld criterion. The experimental and refinement parameters are listed in Table S1 (Electronic Supporting Information, ESI). Maps of electron deformation density, of residual electron density, analysis of the residual density according to Meindl & Henn [19,20] and the DRK-plot obtained via the WinGX suite [21] are also given and discussed in the Supporting Information.

2.4. Computational studies

Quantum chemical calculations were carried out by ORCA quantum chemistry program package [22]. Initial geometry for quantum chemical calculations was taken from the X-ray study and then used to perform subsequent conformational analysis. Mercury program [23] was applied to generate the initial structures for conformational search using molecular mechanics. The structures with the lowest energy were then optimized with B3LYP/def2-TZVPP method/basis combination. The results of conformational studies were used to construct the stacked dimers and to research the stacking interactions with nucleobases. To construct the initial models for the staked dimers the amino groups of the nucleobases were positioned opposite to the plane of monosubstitued phenyl cycle of VMU-2012-05 at the distance of ~ 4 Å. Energies of the interactions between the nucleobase and VMU-2012-05 were calculated as the difference between the total energy of the stacked dimer and its components. To account the dispersion interactions and to obtain more exact geometry on the hybrid DFT level the D3BJ correction [24] was applied. Second derivatives matrices were calculated for all optimized structures to check the absence of positive eigenvalues. Harmonic frequencies and thermodynamic properties were calculated to define the values of zero-point energies and to correct previously calculated interaction energies in stacked dimers. To compare the values of the interaction energies in the stacked dimers obtained on hybrid DFT level, the additional DLPNO-CCDC(T)/6-311G(d,p) calculations were carried out.

2.5. Molecular Docking

Molecular docking was carried out with the Cambridge Structural Database Gold software [25] using the genetic algorithm and the 1S1V structure from the Protein Data Bank (PDB) [26] as a starting point. Prior to the docking procedure, the molecules of external water were removed from the 1S1V structure. After that the hydrogen atoms were added to the protein residues. To score the solutions, the GOLD score function was applied, and the results were then additionally rescored with the Chemscore function.

2.6. QM/MM Calculations

To reveal the strength and the nature of the interactions between VMU-2012-05 and a protein molecule, a simplified model was calculated. The structure of the HIV-1 reverse transcriptase complex with the VMU-2012-05 described above was used as the starting point. The coordinates of 906 atoms from the VMU-2012-05 molecule, the surrounding peptide chains, and the water molecules were marked as an active region in the corresponding pdb file. After that, the QM region of 224 atoms was selected from the active region containing VMU-2012-05 and its closest environment involved in intermolecular interactions. The number of the amino acid residues (18) exceeds typical molecular coordination number and should be sufficient to describe all interactions between VMU-2012-05 and receptor molecules in the selected binding site. The next step was related to the structure optimization. The following combination was applied to optimize the geometry of the models: B3LYP/6-31G(d,p) with D3BJ correction for QM part and CHARMM force field for the rest of the amino acids′ residues. The water molecules in the MM part were described using the TIP3P model [27]. The CHARMM parameters were used for VMU-2012-05 moieties [28,29]. After the optimization, the size of the QM region was expanded up to 448 atoms to include the residues that more distant from the VMU-2012-05 molecule to mimic the distribution of the electron density in the region of the ligand-receptor interactions. The geometry optimization, the hessian and the charge density calculations of the drug–receptor complex were carried out using the ORCA quantum chemistry program package [22]. The charge density functions for the drug–receptor complex were analyzed in terms of QTAIM theory [30], and the NCI method [31] using Multiwfn software [32].

3. Results and Discussion

3.1. X-ray Study and Experimental Electron Density Analysis

Single crystals contain one molecule of VMU-2012-05 in the asymmetric unit in the folded conformation (Figure 1). The molecule remains neutral upon the crystallization with the hydrogen atom solely localized on the nitrogen atom of the heterocycle. This hydrogen atom takes part in the N2-H2...O2 bonding (N...O and H...O distances are 2.847(3) and 1.84 Å, respectively, ∟NOH = 170.3º) to assemble the VMU-2012-05 molecules into the centrosymmetric H-bonded dimers. In addition, there are several weak C-H...O bonds. Analysis of the crystal packing also reveals the stacking interaction between the phenyl group and the uracil moiety of neighboring molecules with slightly inclined disposition of rings (the value of interplanar angle is 5.81(7)°, the distance between centroids of corresponding cycles is equal to 3.398(3) Å).

A characteristic feature of this molecule is that it has five acceptor groups per the only donor group. Note that most of the nucleoside inhibitors of HIV contain several groups that can form strong hydrogen bonds, especially NH, OH, C=O. Therefore, among all non-covalent interactions the strong hydrogen bonds play a key role for the stabilization of the crystal packing and the realization of the ligand-receptor complexes [33,34]. Instead, non-nucleoside HIV inhibitors usually do not possess potential donor groups suitable for appearance of the strong hydrogen bonds in the crystals. As a result, their crystalline packing is stabilized mainly by weak C–H...O bonds and van-der-Waals H...H interactions, while the contribution of strong hydrogen bonds to the lattice energy here is small. Recently, we demonstrated that the contribution of the hydrophobic interactions to the binding energy in the complex between a protein and lamivudine is comparable to that for the lattice energy of the tetragonal polymorph II of free base lamivudine despite the discrepancies in the molecular conformation, while the energies of the hydrophilic interactions decrease [35]. For antitumor agent bicalutamide, both contributions and the total interaction energies were found to be nearly the same both for the experimental charge density studies and the QM/MM calculations [36].

We performed the detailed study of the intermolecular interactions in the crystal of VMU-2012-05 within the formalism of the Quantum Theory “Atoms in Molecules” (QTAIM) formulated by R. Bader, as well as using the Non-Covalent Interactions (NCI) method, using the results of the multipole refinement based on the high-resolution X-ray diffraction data. To evaluate the energy of the weak interactions we used the Espinosa-Mollins-Lecomte correlation formula that allows one to calculate required values using the potential energy density V(r) in bond critical points (bcps) as E_int = 0.5V(r) [37]. In addition, another formula E_int = 0.429G(r) has been used for the evaluation of the hydrogen bonds and the weak interaction energies via kinetic energy density G(r) [38]. The corresponding values are given in the ESI.

The analysis of the experimental charge density function in the crystal of VMU-2012-05 reveals bcps related to all covalent bonds and the weak intra-/intermolecular interactions. All covalent bonds are characterized by the negative values of Laplacian and local energy density that corresponds to the shared type of interactions (detailed information can be found in ESI, Table S1 and Figs. S1 – S5). At the same time, weak intra- and intermolecular interactions correspond to the closed-shell type of interactions with the positive values of Laplacian and local energy density in respective bcps. QTAIM study reveals three intramolecular interactions, namely, the C–H...O bond between the methyl group bonded to the uracil fragment and bridging oxygen atom (C20–H202...O3), as well as two unusual contacts, O4...C4 and C19...C20. The first contact can be described as a peak-hole interaction, while the second one most probably belongs to weak dispersion interactions. These findings are in good agreement with the evaluated energies of the interaction (-3.7 and -1.5 kJ/mol). The hydrogen bond N2–H2...O2 was found to be the strongest intermolecular interaction, (-29.5 kJ/mol) while the energies of C–H...O and C–H...N bonds are much weaker (-2.3 – -7.6 kJ/mol and -2.2 – -3.7 kJ/mol). At the same time, only one bcp (C16...N1) can be associated with the stacking interaction, so its strength was estimated as -2.7 kJ/mol. The sum of all intermolecular interactions, i.e. lattice energy (E_latt), was found to be -214.6 kJ/mol. The sum of the partial contributions of hydrophilic C–H...O and C–H...N weak interactions gives more than 50% (53.1 and 5.5%) of E_latt. Thus, the energy and the characteristics of the intermolecular interactions found in the crystal under study are typical for organic compounds. Comparison of the molecular conformations in crystal and isolated molecules, and analysis of the energies of rotational barriers and intermolecular interactions could allow revealing which of interactions mainly affect the conformation in a condensed matter.

3.2. Molecular Flexibility and Conformational Analysis

Conformation of small flexible molecules in ligand-receptor complexes with biomacromolecules typically differs from that in their solid forms (polymorphs, solvates, cocrystals and salts). We can claim with a high degree of accuracy that an extended or linear conformation is preferable for ligand-receptor complexes, while in crystals of small molecules folded or globular conformations are favorable. Besides, prevalence of particular type of interactions can stabilize an extended or a folded conformation in both crystals and ligand-receptor complexes [39]. The VMU-2012-05 molecule contains three rigid aromatic groups connected by several single bonds that allow to change the mutual orientation of the former groups. As a result, the conformational energy can be described as a function of six dihedral angles (Scheme 1). Accordingly, the conformational analysis of the isolated VMU-2012-05 molecule reveals four stable conformers (Table 1). The energy of the least favorable conformation exceeds the energy of the most favorable one by 14.23 kJ/mol. The majority of the conformers are in folded conformation except for the conformer3 realize. Moreover, for all of them the conformation significantly differs from that found in the crystal structure (Table 1, Figure 2a).

In the case of VMU-2012-05, the conformational changes can be expressed by the values of ψ1-ψ6 angles depicted in Scheme 1 (the curves of corresponding rotation potentials are presented in Figures S6–S11 in the ESI). The rotation barriers around the bonds in the oxyethyl moiety exceed 60 kJ/mol. At the same time, the rotation barrier that corresponds to the mutual displacement of two phenyl groups around two bridging C–C bonds is also too high for free rotation in solution (up to 32 kJ/mol). Upon transition of the most stable conformer4 to the less favorable conformer1 the ψ angles are changing not significantly, with the exception of ψ₄. Nevertheless, the change in the ψ₄ angle from -7.1 to -109.0° increases the full energy by less than 20 kJ/mol. On the other hand, more than 60 kJ/mol is necessary to change the ψ₃ angle by almost 360° (almost complete rotation around the O–C bond in oxyethyl moiety) and, therefore, to change significantly the mutual orientation of uracil and phenyl moieties. In the case of transition from conformer4 to another unfavorable conformer3 it is necessary to overcome several potential barriers as high as 35 kJ/mol. Thus, the conformational changes in VMU-2012-05 are not fully forbidden but the energies of potential barriers significantly exceed experimental values of individual intermolecular interactions in crystal besides the hydrogen bond N2–H2...O2. Hence, a strong hydrogen bond is able to stabilize one of unfavorable conformations in the condensed matter, but the sum of weak interactions occurred between a pair of molecules can also act in a similar manner.

3.3. Modelling the Interactions of VMU-2012-05 in Dimers with Nucleobases

According to the results of macromolecular crystallographic studies, for instance [40,41,42], the molecules of nucleoside inhibitors in its phosphate forms can bind with the fragments of DNA or RNA due to the stacking interactions between π-systems of corresponding nucleobases and the aromatic heterocycle of nucleoside fragment. Here the heterocycle is a chemically modified nucleobase fragment usually situated on top of the RNA or DNA fragment, so its’ role can be identified as the key factor that prevents the process of RNA or DNA transcription. The existence of such interactions in the case of non-nucleoside inhibitors of HIV is an open question because in the most of the complexes with the HIV inhibitors the corresponding molecules were found in the region far from any DNA or RNA fragments. According to ref. [43] the role of the non-nucleoside inhibitor is to lock the HIV polymerase in an inactive conformation by means of weak intermolecular interactions. On the other hand, numerous examples of stacking interactions of the uracil fragment with nucleosides are observed in the Cambridge Structural Database. In order to evaluate stability of stacking interactions between VMU-2015-05 and DNA/RNA we optimized dimers constructed by the stacking interaction between the 5-methyluracil fragment of one molecule of VMU-2012-05 with one molecule of adenine, cytosine, guanine, thymine or uracil.

Initial geometry of VMU-2015-05 molecule in the stacked dimers was the most stable conformation4. After optimization of geometry the conformation of VMU-2012-05 has changed due to rotation of phenyl group around C–C bonds that correspond to changes of ψ₅ and ψ₆ torsion angles (Figure 2b). However, the dimers remain to be connected by stacking interactions between the nucleoside cycles (Figure 3). The analysis of the potential energy surface reveals two stable stacked dimers for adenosine and only one for the other nucleobases. Two stacked dimers with the adenosine molecule differ by the nature of the stacking interactions. In the first stacked dimer, VMU-2012-05·with adenine (dimer1), we found the stacking interaction between the uracil moiety of VMU-2012-05 and the pyrimidine cycle of adenosine. In the second stacked dimer, VMU-2012-05· adenine (dimer2) the 1,3-pyrazol cycle of adenosine participates in stacking interaction instead of pyrimidine one.

In general, the mutual orientations of the uracil moiety and the heterocyclic fragment of nucleobases are non-parallel. The characteristic angle φ between the planes of the cycles involved in stacking interaction, varies in the wide range (1.968 – 12.923°). At the same time, the values of the distances between the centroids of the involved cycles and the inter-centroid shifts fall in the ranges 3.51 – 3.67 and 0.41 – 1.80 Å, respectively, that is indicative for rather high extent of overlapping of the corresponding π-systems. Energies of dissociation of VMU-2012-05·adenine (dimer2), VMU-2012-05·cytosine and VMU-2012-05·guanine dimers are significantly lower than those of other stacking dimers (Table 2). The QTAIM analysis of stacked dimers demonstrates that the additional stabilization is provided by the hydrogen bond between the amino group bonded to the pyrimidine cycle of nucleobase and the carbonyl group of VMU-2012-05(Tables S4–S6, ESI). The effect of additional stabilization due to hydrogen bonding is assessed differently by B3LYP and DLPNO-CCDC(T) methods in the case of dimer1 and dimer2. In the case of the other stacked dimers, both methods gave nearly the same results.

Thus, the results of the calculations of the isolated stacked dimers of VMU-2012-05 allowed us to assume that the binding of this nonnucleoside anti-HIV agent with DNA or RNA can be rather strong even at absence of H-bonding. However, to afford the possibility of such interaction, the molecule of VMU-2012-05 should noticeably change conformation.

3.4. Modelling the Interactions of VMU-2012-05 with Receptor

The majority of previously reported reverse transcriptase (RT) inhibitors, both based on uracyl derivatives [44,45,46,47,48,49], and others [50,51,52] were reported to occupy the same binding site forming a hydrogen bond with Lys101 and a number of hydrophobic interactions with Tyr181, Tyr188 and Trp229 residues. Taking into account that the VMU-2012-05 belongs to 1,6-disubstituted uracyles, we can assume that the VMU-2012-05 molecule can occupy the same binding site as TNK-651 [45], HEF [46], KRV [47], MKC [48] and GCA [49] ligands. To mimic the interaction of VMU-2012-05 with a protein receptor, the 1S1V structure of L100I mutant HIV-1 reverse transcriptase as receptor with its potential inhibitor TNK-651 (6-benzyl-1-benzyloxymethyl-5-isopropyluracyl) [44] was taken from PDB as a starting model.

The overall conformation of VMU-2012-05 molecule after optimization is close to those in the most unfavorable conformer1 with an exception of rotation along the N1–C5 bond (Table 1, Figure 2). The complex resulting from QM/MM calculations possesses intermolecular distances between the receptor and VMU-2012-05 typical for molecular crystals of organic compounds. Thus, the geometry of the calculated model is suitable for the subsequent analysis of the interactions between VMU-2012-05 and HIV-1 reverse transcriptase in terms of QTAIM theory and NCI method. We found three strong hydrogen bonds between the uracil moiety of ligand and the receptor residues (Figure 4). These are N-H...O bonds where VMU-2012-05 acts as a donor of N-H or acceptor of O=C groups for H-bonding with the LYS103 (instead of expected LYS101), VAL179 and VAL106 residues. Interesting that in X-rayed RT complexes with inhibitors, the amide groups of VAL179, VAL106 and SER105 residues are often do not take part in intraprotein H-bonds or take part in H-bonding with water molecules. Thus, geometry optimization affected not only ligand, but also macromolecular 3D conformations. The other interactions between VMU-2012-05 and the peptide chains are hydrophobic as expected (Figure 4, Table S16, ESI). All H...O distances corresponding to the C-H...O interactions between ligand and receptor vary in the range 2.39-2.51 Å. At that, the most of the H...H distances between ligand and receptor exceed 2.14 Å, only one short H...H distance (2.022 Å) is observed that can correspond to the interaction between methyl group of isoleucine residue ILE96 and the methylene bridge of VMU-2012-05. No stacking interactions were found in the complex.

To calculate the electron density function, the QM region was expanded thus to cover all possible interactions of VMU-2012-05 with the surrounding aminoacid residues. Subsequent study of the calculated electron density in terms of QTAIM theory has revealed that the strongest interactions are the hydrogen N-H...O bonds between uracil moiety of VMU-2012-05 and LYS103/VAL179/VAL106 residues, respectively (-23.7, -15.6 and -13.8 kJ/mol, Table S16, ESI). Note, that these interactions in previously reported complexes All interactions of other types (C-H...O, C-H...N, C-H...π, H...H, O...O, O...N) are noticeably weaker (-0.7 – -12.5 kJ/mol). The sum of energies of all interatomic interactions, i.e. ligand-receptor binding energy, is equal to -252.5 kJ/mol, which is 37.9 kJ/mol lower than the value of lattice energy. This discrepancy comes mainly from the additional H…H bonding and O…O plus the O…N bonding observed only in the complex. As a result, total energy of the hydrophobic interactions was higher in the complex, while the total energy of the hydrophilic hydrogen bonds decreased probably due to the larger cage occupied by VMU-2012-05. Indeed, in the crystal and the complex the volume of the molecular Voronoi polyhedron was found to be 416 and 446 ų. The Voronoi tessellation also indicates increase in the contribution of hydrophobic interactions in complex as compared with crystal, although in general it underestimates the role of hydrophilic interactions (Figure 5b).

3.5. NCI Study of Intermolecular Interaction in Crystal, Stacking Dimers and the Model of Ligand-Receptor Interaction

NCI analysis was performed for pure VMU-2012-5, its′ dimers with nucleobases, and in the ligand–receptor complex. The study of weak interactions in terms of NCI and the comparison with the results of QTAIM study reveals qualitative agreement between the positions of bcps and RDG maxima. In addition, the shape and the volume of RDG maxima gives us complimentary information about the area of interaction between the atomic domains. In crystal of VMU-2012-05 the stacking interaction corresponds to the shapeless planar isosurface (Figure 6a). The most of the isosurface area is colored by light brown that is characteristic for the van-der-Waals interactions. Nevertheless, approximately one fourth of the corresponding area is light blue that reveals the attractive interactions between two atomic pairs (N2...C2 and H...C) of uracil and phenyl moieties, respectively. The intramolecular C-H...O bond can be characterized as mostly van-der-Waals type. In contrast to the stacking interaction and the C-H...O bond, the attractive nature is quite expectedly demonstrated for intermolecular N-H...O bonds (Figure 6b). In the dimers between VMU-2012-05 and nucleobases the study of RDG clearly reveals the role of stacking interaction and strong hydrogen bond in its stabilization. The stacking interaction in all studied dimers is characterized by rather large interaction area (Figures S27–S29, ESI). For instance, in dimer1 the stacking interaction can be distinguished as number of van-der-Waals interactions between the nitrogen and the carbon atoms of the adenine and the uracil moiety of VMU-2012-05 assisted by the weak C-H...N bond. In more stable dimers VMU-2012-05·adenine (dimer2) and VMU-2012-05·guanine the stacking interaction is assisted by one strong N-H...O and several weak H-bonds of other types. In addition, several intramolecular interactions in VMU-2012-05 are observed that can be attributed as the C-H…π and C-H...O interactions.

The binding of VMU-2012-05 with model HIV-1 reverse transcriptase is somewhat different as compared to the corresponding crystal structure and the stacked dimers. The VMU-2012-05 molecule participates in several H-bonds of attractive character including two C-H...O bonds. At the same time, there are a lot of interactions of van-der-Waals type that play key role in the ligand-receptor binding (Figure S30, ESI).

Similar visualization of the intermolecular surface between the interacted molecules can be obtained from the molecular Voronoi surfaces. It also reveals both stacking and H-bonded pairwise interactions both in the crystal (Figure 7) and in the ligand-receptor complex (Figure 8). The view of surfaces is similar with that obtained from the RDG function, demonstrating the potential of the molecular Voronoi surfaces for fast low-resource demanding analysis of the closest environment of molecules in any crystal, solid and associate.

4. Conclusions

In this work we demonstrate that even for molecules with large number of potential donors and acceptors of hydrogen bonding, the role of weak van der Waals interactions in the stabilization of its associates with biomacromolecules should not be underestimated. Sum of energies of these interactions can be large enough to stabilize an unfavorable conformation of a flexible molecule, and is often higher than energy of a strong, but individual interaction such as a hydrogen bond. For VMU-2012-05 molecule, the conformation analysis reveals four conformers with relative energies 0 – 14.2 kJ/mol. None of them were found in the condensed phases, either experimental (crystal), or theoretical (six dimers with nucleosides and a complex with HIV-1 reverse transcriptase). The experimental folded conformation is supported by the C–H…O, C…O and C…C intramolecular interactions with the overall energy of -13.0 kJ/mol, while all theoretical conformations are extended due to the rotation along oxyethyl moiety. Analysis of VMU-2012-05 : nucleoside associates demonstrates that overall (total) energies of C…C, N…C and C…H interactions between the nucleoside and the phenyl rings can reach -14.0 kJ/mol that is comparable with the energy of N–H…O hydrogen bonds observed in the crystal and calculated for the complex with HIV-1 reverse transcriptase (-7.5 – -29.5 kJ/mol). Counterintuitively, despite the number of H-bonds in the ligand–receptor complex exceeds that in the crystal due to the presence of numerous donors of H-bonding in the binding pocket of a macromolecule, the overall energy of N-H…O bonds in the crystal and in a complex is nearly the same. Stacking interactions in the binding pocket are also absent. Instead, the sum energy of H…H interactions which are typically ignored during discussion of structure-forming contacts increases from -0.9 to -61.7 kJ/mol instead. This fact indicates necessity of comparative analysis of full continuum of intermolecular interactions in both materials and biomolecules, that at the moment can be performed only by means of ED-based methods.