Under ONIOM Layers: Analysis of BCR-ABL Enzyme Inhibitors Through Bond-Critical Points and Natural Orbitals

TL;DR

This study compares how rebastinib and ponatinib interact with the BCR-ABL enzyme using advanced computational methods to assess their effectiveness.

Contribution

The paper introduces a detailed electronic structure analysis of BCR-ABL enzyme inhibitors using ONIOM and BCPs.

Findings

Ponatinib forms stronger interactions with key residues like Glu286, Met318, and Asp381 compared to rebastinib.

Ponatinib shows a larger HOMO-LUMO gap, suggesting greater chemical stability and potential to inhibit catalytic activity.

ONIOM calculations revealed significant differences in the electronic properties of the two inhibitors.

Abstract

Considering the relevance of hydrogen bonds and other intermolecular interactions in regulating the activity of the tyrosine kinase class of enzymes, an in-depth electronic structure study of these forces in the context of the BCR-ABL protein was performed through full optimizations using the ONIOM method. Rebastinib and ponatinib were docked to the target enzyme using AutoDock Vina to provide starting-point geometries, which were then optimized through ONIOM calculations. This study evaluated Frontier Molecular Orbitals (FMOs) and Bond-Critical Points (BCPs) located in the sites of interactions formed with accessible residues, such as Glu286, Met318, and Asp381. Ponatinib’s ONIOM-optimized structure was shown to not only form and preserve prominent interactions, which were shown to be significantly stronger than those formed by rebastinib, but also to be associated with a significant increase in the HOMO (Highest Occupied Molecular Orbital)−LUMO (Lowest Unoccupied Molecular Orbital) gap, indicating its potential to hinder catalytic activity by providing higher chemical stability when compared to rebastinib.

Figures

Full text

1. Introduction

The BCR-ABL transcript is part of the tyrosine kinase class of enzymes, being related to the development of Chronic Myeloid Leukemia (CML) and some forms of severe Acute Lymphoblastic Leukemia (ALL) [1,2,3,4,5,6,7,8,9]. This enzyme has its origin in a mutual translocation of genetic material involving the 9th and 22nd pairs of chromosomes, which produces a chimeric arrangement—the Philadelphia chromosome—that contains the mutant BCR-ABL gene [10,11,12,13]. In the context of CML, the detection of this gene through karyotyping or in situ hybridization tends to be among the best techniques for diagnosis [1,14,15,16].

The oncogenic activity begins with an essential phosphorylation of the Tyr177 residue, which promotes the recruitment of other cytoplasmic proteins to form a complex structure with multiple protein subunits [17]. The first of these proteins is the Growth Factor Receptor-Bound Protein 2 (GRB2), which joins BCR-ABL and results in a complex capable of recruiting the Son of Sevenless (SOS) factor [9]. The BCR-ABL/GRB2/SOS complex promotes activation of the Rat Sarcoma (Ras) factors through its conversion to the GTP-bound form. Thus, the final consequence is the activation of the replication cascade mediated by Mitogen-Activated Proteins (MAPs) and mitogen-activated protein kinases signals [3,18,19,20,21]. Other oncogenic mechanisms were also identified, such as the hindrance of cell adhesion, an increase in the reactive oxygen species concentration, and the inhibition of apoptosis [22,23,24,25].

The primary successful treatment strategy for CML is the use of tyrosine kinase inhibitors. Among these, imatinib was considered the gold-standard CML drug, being associated with a survival rate of over 80% at 10 years [26,27,28,29]. However, further mutations emerged and proved capable of disrupting its inhibitory and therapeutic activity [30,31].

There is a series of point mutations that result in different modulations of the structural profile of the enzyme and its response to a variety of drugs, including imatinib [30,32]. Among these, the T315I mutation is notorious for being associated with the worst cases of resistance by intensely disturbing the interaction network of first-choice inhibitors with the target enzyme [2,30]. This resistance conference is mainly associated with a strong stabilization of the active conformation of the BCR-ABL enzyme [33]. Access to the hydrophobic pocket is also significantly hampered, as isoleucine has a larger side chain than threonine [34], resulting in additional steric bumps.

To tackle these obstacles, ponatinib was designed as part of the third generation of tyrosine kinase inhibitors [35,36]. Its structure was explicitly designed to penetrate deep into the hydrophobic pocket of the target enzyme, resulting in the effective inhibition of the catalytic activity [34]. However, although ponatinib proved to be a valuable option for the treatment of CML when the T315I mutation is present, it was also shown to be significantly cardiotoxic [37,38].

Other tyrosine kinase inhibitors have been studied as potential candidates for the treatment of CML, including rebastinib, originally a tyrosine kinase inhibitor targeting the Tie2 series of protein kinases, which is associated with angiogenesis and is overexpressed in various tumors [39,40,41,42]. Although rebastinib was initially used for the inhibition of Tie2 kinase proteins, recent studies have shown its potential as an alternative drug for the treatment of CML [43,44,45], since it has been associated with effective inhibitory activity due to the presence of molecular scaffolds that interact with the ATP-binding region and the Glu282 and Arg386 residues in its structure [46]. This dual capacity is associated with an additional effect of conformational control that could assist in traditional competitive inhibition, resulting in highly potent inhibition of the T315I isoform of the BCR-ABL tyrosine kinase [46].

Given that attractive interactions, especially hydrogen bonds, are intimately related to the structure, control, regulation, and inhibition of kinase-type enzymes and other receptors [47,48,49,50,51,52,53], an in-depth study of rebastinib and ponatinib interactions with the BCR-ABL protein was performed to better understand and indicate different aspects in their electronic interactions with the target enzyme that may translate to their different therapeutic effects and inhibitory capacities. This context was achieved through full ONIOM optimizations of the entire protein complexes, as well as the Quantum Theory of Atoms in Molecules (QTAIM) and Natural Bond Orbital (NBO) analysis, which enabled intricate evaluations of large biological systems. Therefore, this study explores the electronic interaction of rebastinib and ponatinib with the BCR-ABL protein, focusing on the role of hydrogen bonding in the kinase-type enzyme structure.

2. Results and Discussion

The structures of ponatinib and rebastinib, which are deposited in PDB codes 3IK3 and 3QRJ, respectively, were used to validate the optimized B3LYP and docking structures. To assess the differences in the geometries, we calculated the root mean square deviation (RMSD) values of the docking and B3LYP geometries of ponatinib in relation to the PDB of 3IK3, which are 0.7766 and 1.0520 Å, respectively. In the case of rebastinib, the RMSD values are smaller. The RMSDs of the docking and B3LYP geometries in relation to the PDB of rebastinib in 3QRJ are 0.0098 and 0.5089 Å, respectively. The RMSD between docking and optimized B3LYP geometries for ponatinib is 0.9869, and for rebastinib, it is 0.5099. Therefore, RMSD values smaller than 2.0 Å were obtained, which is a reasonable threshold for ligand-binding conformations [54,55].

We used the optimized geometries of the complexes to evaluate the distances and angles of the main hydrogen bonds involving residues Glu286, Met318, and Asp381—more specifically, the ionic oxygen atoms from the lateral chains of the Glu286 and the Asp381 residues that interact with polar hydrogen atoms of the TKIs. These residues are described in the literature as important regulators of the active conformation of the BCR-ABL enzyme [46,56,57,58,59]. Therefore, a study of the the electronic structure of their interactions with ligands can point out important aspects for the inhibition of the target protein. Furthermore, the energetic aspects of the systems were evaluated—specifically, the energies of the obtained frontier orbitals.

Figure 1 depicts the 2D structures and residues Glu286, Met318, and Asp381. The snapshots of NCI profiles of rebastinib and ponatinib are included. The complete NCI isosurface is shown in the Supplemental Materials.

Table 1 compiles the evaluated data of frontier molecular orbitals of BCR-ABL and BCR-ABL:ligand, which can be used to predict reactivity properties [60,61]. The gap value is the difference between the energy values of LUMO and HOMO ( $[eqn]$ ). Gap values are dependent on the choice of density functional approximation (DFA), and it is well known that B3LYP underestimates the gap values [62]. However, the goal of this work is to determine the trends of frontier orbitals, which is effectively achieved with DFAs, not the gap’s absolute values [61,63,64,65,66,67].

It is important to note that the energies of the molecular orbitals, as well as the gap value, are only a pair of the indicators of stability and reactivity and that it is also necessary to consider other aspects related to the interaction of the protein with the studied ligands. The distances and angles of intermolecular interactions may indicate their intensity [68] and may also be related to the stabilization of the complex through this energy balance.

The B3LYP binding energy (BE) was calculated using the method of counterpoise correction proposed by Boys and Bernardi to evaluate the known Basis-Set Superposition Error (BSSE) [69]. The use of the BSSE is important to take into account for this small basis set used to compensate for the dispersion interaction [70]. The binding energy of rebastinib is 39.3 kcal·mol^−1^, which is smaller than the BE of ponatinib of 209.0 kcal·mol^−1^. However, in order to take into account the dispersion correction using the Grimme D3 dispersion correction, based on the B3LYP geometry, the BEs are 86.7 and 248.5 kcal·mol^−1^ for ponatinib and rebastinib, respectively. Moreover, since vibrational frequencies were used to characterize the minimum-energy DFT calculations, we also accounted for thermal corrections to the energies of these stable complexations. The Gibbs free energies for ponatinib and rebastinib are 81.5 kcal·mol^−1^ and 310.6 kcal·mol^−1^.

Table 2 shows the geometries of the key attractive interactions formed by rebastinib and ponatinib with the BCR-ABL enzyme. Many intermolecular interactions—especially hydrogen bonds—tend to achieve optimal geometries and, consequently, the highest stability at shorter distances and angles close to 180° [68].

After considering the initial information provided by the analysis of the geometry of each interaction, Table 3 yields further insight using the QTAIM and NBO data related to the interactions formed with Glu286, Met318, and Asp381. The interactions formed by ponatinib tend to be stronger. The hydrogen bond established with the Glu286 side chain surpasses both of its counterparts in the rebastinib complex, with the Laplacian values indicating that the electronic density of the interaction formed by ponatinib is closer to the atomic center.

According to Figure 1, the NCI profiles of rebastinib and ponatinib show different behaviors. It can be observed that ponatinib forms an additional and strong attractive interaction with the negatively charged side chain of Asp381, which is indicated by the dark-blue highlighted isosurface. Although rebastinib forms a hydrogen bond with Met318, it is suggested to be significantly weaker, as indicated by its isosurface. The absence of the contributing effects of an aromatic scaffold similar to that of ponatinib may be a cause for a weaker interaction [71]. As mentioned, both ligands also form attractive interactions with the Glu286 side chain.

Ponatinib has the highest gap value (Table 1) and is known as the best-performing molecule when applied to the BCR-ABL enzyme, which was expected, since its structure was designed for application to this protein [35]. Rebastinib also showed an increase in the gap value, but it was modest when compared to ponatinib while alternatively showing an increased value for the binding energy. This data corroborates its potential as an effective BCR-ABL inhibitor, which is described to form additional interactions with residues that ponatinib does not—specifically, Glu282 and Arg386 [46]. However, despite these extra interactions, it can be assumed, considering the evidence provided by the frontier molecular orbitals, that ponatinib may provide greater stabilization to the enzyme when examined through this method.

In relation to the docking structure, rebastinib improves its interaction distance (Table 2) with Glu286 through both hydrogen bonds identified along the ONIOM-optimized geometry. These interactions are depicted in Figure 1 as rebastinib:Glu286(i) and rebastinib:Glu286(ii), both involving the negatively charged carboxylate group of the Glu286 residue and hydrogen-bond donors in the form of the hydrogens bound to the nitrogens in the carbamide functional group of the rebastinib structure.

The rebastinib:Glu286(i) interaction had its length decreased by 0.119 Å while increasing the established angle by 10.46°. Alternatively, the rebastinib:Glu286(ii) interaction had its length decreased by 0.117 Å and its angle increased by 18.78°. An improvement in the geometric configuration of both interactions could be observed, indicating a strengthening of the attractive forces, given the influence of distances and angles on hydrogen-bond strengths [68]. The ONIOM-optimized lengths of the rebastinib:Glu286 interactions were shortened in relation to the docking and experimental data [46].

Ponatinib also forms an interaction with the carboxylate group of the Glu286 residue (Figure 1) using the N-H scaffold of its amide functional group as the hydrogen-bond donor. This interaction was labeled ponatinib:Glu286, and it reduces its distance to the side chain of the Glu286 residue throughout the geometry optimization by 0.041 Å while simultaneously increasing its interaction angle by 2.71°. It is worth noting that, similarly to what was observed in the interactions between the Glu286 residue and rebastinib, the distance and angle values became geometrically more favorable for strong hydrogen bonding, also indicating the formation of a stronger intermolecular interaction when compared to rebastinib due to the shorter length and larger angle. Additionally, both ligands were also shown to differ when their interactions with Asp381 were evaluated. Here, ponatinib reaches an equilibrium with an interaction distance only 0.017 Å larger than the one present in the starting geometry, and the interaction angle also decreased by a small value of 1.25°.

Another important point is the behavior of the attractive interactions formed between Met318 and the studied ligands. It can be observed that rebastinib forms two hydrogen bonds, both of which involve the polar backbone of the residue. The rebastinib:Met318(i) interaction is formed between the carbonyl group of the methionine structure and the N-H scaffold of the amide functional group of the ligand. Its length is reduced by 0.074 Å with optimization while also increasing the interacting angle by approximately 3°. A reduction in the interaction distance was not found in the case of the rebastinib:Met318(ii) interaction, which is formed between the amine function of the methionine residue and the nitrogen contained in the pyridine ring of rebastinib. This interaction shows an increase of 0.153 Å, while the angle also increased by 6.876°. The optimized structure interaction shows a geometry with an interacting angle closer to an ideal linear interaction and a shorter length compared to the crystallized rebastinib in complex with the BCR-ABL enzyme, specifically in the case of rebastinib:Met318(ii) [46].

The NCI isosurface for the rebastinib:Met318(ii) shown in Figure 1 is more intense compared to rebastinib:Met318(i), which can be associated with its geometry and the involvement of nitrogen in an aromatic ring, possibly increasing the attractive forces [71,72]. A similar effect may be associated with the hydrogen bond formed between the Met318 residue and ponatinib (ponatinib:Met318), which involves the amine functional group of the methionine backbone as a hydrogen-bond donor and nitrogen contained in the imidazopyridazine scaffold of the ligand. This interaction is well established as one of the most important in this system, as it provides access to the hydrophobic pocket of the BCR-ABL enzyme and ponatinib is particularly designed to favor this hydrogen bond [34,35]. This interaction increased its length by 0.002 Å and its angle by 0.25°, resulting in an optimized geometry that is more favorable for the formation of an intense hydrogen bond [68] than what is observed in both interactions involving the Met318 residue and rebastinib.

As for the interactions formed with the Asp381 residue in the DFG switch, it is possible to observe that rebastinib forms a single hydrogen bond, while ponatinib forms two. Rebastinib interacts with the Asp381 residue through the C=O scaffold of its carbamide functional group. This hydrogen bond was labeled rebastinib:Asp381, and its length increased by 0.196 Å in the ONIOM-optimized geometry. An increase was also observed in the interaction angle, which was 7.68° larger in the ONIOM-optimized geometry. Alternatively, ponatinib formed a conventional hydrogen bond through the oxygen of its amide functional group, which interacts with the amine function of the backbone of the Asp381 residue, ponatinib:Asp381(i). This interaction is found to be relatively conserved in the ONIOM-optimized geometry, with a modest increase in its length measuring 0.017 Å and a decrease of the interaction angle of 1.25°.

The additional interaction formed between ponatinib and the Asp381 residue, ponatinib:Asp381(ii), is an unconventional hydrogen bond established by the resulting attractive force between a partially charged methylene group in the 4-methylpiperazine ring and the carboxylate functional group in the side chain of the Asp381 residue. This interaction was made possible by the binding of the methylene group to a heavily polar nitrogen atom. As shown by Figure 1, this interaction is significantly intense and may be associated with effective inhibition of the BCR-ABL enzyme through the DFG switch [35,46,56].

Rebastinib shows an electron density distribution (Table 3) that is slightly more concentrated over the length of the interactions with Glu286 due to the smaller Laplacian value when compared to that of the ponatinib:Glu286 interaction, with a difference of 0.217 e·Å^−5^, as can be observed in Table 3. This difference indicates a stronger interaction that approaches a covalent interaction. This effect was shown to be modestly influential on the energy density (H^c^). The interaction formed with Glu286 by rebastinib and ponatinib has negative values, and following the concept of H^c^ introduced by Cremer and Kraka, negative values correspond to covalent interaction, while positive values show the electrostatic character of the weak interaction.

Alternatively, according to the NBO analysis, the smaller Laplacian value in the rebastinib:Glu286 interactions and consequent superior superposition of electronic density distribution along the interaction length could be associated with the higher value of the stabilization energy that the rebastinib:Glu286(ii) interaction shows when compared to ponatinib. However, the electron density observed in this interaction is smaller in rebastinib, albeit by a small margin.

Furthermore, according to Table 3, both NBO and BCP were found in the interactions formed with Met318 by ponatinib and rebastinib. As indicated by Figure 1 and Table 2, ponatinib forms a stronger interaction with the aforementioned residue, with a measured value for kinetic energy significantly larger than what was observed in both interactions between rebastinib and the Met318 residue, as well as almost double the electron density. The stabilization energy for the NBO in this interaction was measured as 19.96 kcal·mol^−1^, which is also significantly higher than that shown by rebastinib, which only reached 1.84 kcal·mol^−1^ in rebastinib:Met318(i) and 4.35 kcal·mol^−1^ in rebastinib:Met318(ii). All of the mentioned orbitals are of the LP-BD* type. This high energy difference can be explained by the intricate interplay involved in the electronic network of the interacting region, such as the fact that the LP-BD* hydrogen bond is enhanced by a larger orbital overlap and smaller orbital gap [71] when compared to its counterparts. It is associated with a more favorable orbital geometry under the lens of the Pauli repulsion [72,73]. It is also worth noting that the hydrogen-bond acceptor in this interaction is nitrogen contained within an intensely aromatic scaffold, an imidazopyridazine ring in conjugation with a triple bond, which is in resonance with another aromatic ring. This high aromaticity is a significant reinforcement of the hydrogen bond through $[eqn]$ resonance assistance [71,74,75].

As for the interactions with the Asp381 residue, the BCP also indicates stronger interactions in the complex formed by ponatinib. Both interactions established by this ligand have higher kinetic energy than that of the one formed by rebastinib, with an associated higher electronic density. The Laplacian values indicate an electronic density more distributed over the interaction length in the rebastinib:Asp381 interaction when compared to ponatinib:Asp381(i) but not when compared to ponatinib:Asp381(ii). Additionally, two NBOs of importance were identified in the complex formed by ponatinib, while only one was identified in the in the complex formed with rebastinib. The NBOs identified in the former are of the LP-BD* type in both the ponatinib:Asp381(i) and ponatinib:Asp381(ii) interactions, with the first showing a stabilization energy of 7.11 kcal·mol^−1^, while the other has an associated energy of 3.95 kcal·mol^−1^. All of these NBOs showed a higher energy values when compared to the one in the rebastinib:Asp381 interaction; they are individually intense and are additionally subject to the cooperative effects that reinforce interactions in a network, resulting in a general attractive force that is greater than the sum of the contributions of each attractive interaction [71,76,77,78,79].

Considering all of the presented data, it is possible to observe that the binding mode of ponatinib as obtained by the AutoDock Vina algorithm is very close to that which results from further optimization via ONIOM calculations, and it is also the one associated with the formation of stronger attractive interactions, such as the hydrogen bonds formed with the Glu286, Met318, and Asp381 residues. For the Met318 residue specifically, the stabilization energy that could be identified in the interaction with ponatinib indicates a significantly intense and stable attractive force, which has also been described in the literature as one of the driving forces for the effectiveness of ponatinib through its influence on the affinity and consequent access to the hydrophobic pocket of the BCR-ABL enzyme [34,35,36,80]. Although rebastinib is a ligand with affinity for the same binding site as rebastinib, its capacity to interact with the Glu282 and Arg386 dyad, is located far from the Met318 site [46], may result in a weakening effect regarding the entrance into the hydrophobic pocket and may result in poorer performance regarding the competitive inhibition of the BCR-ABL enzyme.

The additional attractive interaction can also be seen in Figure 2, which shows the RDG plots of the interactions formed between the ligands and the target enzyme. In these plots, sign $[eqn]$ represents the distribution of electronic density close to the formed interaction, with sign $[eqn]$ indicating high electronic densities and, therefore, attractive interactions [81,82,83,84].

According to Figure 2, the attractive interactions formed in the BCR-ABL + Ponatinib complex are more numerous and tend to be stronger. In contrast, the steric clashes and repulsive interactions indicated by the red spikes in the sign $[eqn]$ region are more similar to the profile exhibited by the isolated BCR-ABL enzyme, which may indicate that the interactions formed by ponatinib are less disruptive for the system as a whole, possibly due to its rigid and linear nature [35].

3. Methods

3.1. Docking Protocol

Aiming to obtain the starting geometries for the QM/MM calculations, the T315I isoform of the human BCR-ABL enzyme in the DFG-out conformation [46,59,85] was chosen as the target protein. Its crystal structure was obtained from the Protein Data Bank under code 3QRJ [46,86]. The structure was prepared for calculations with the removal of waters, the addition of hydrogens, the removal of doubled residues, and the repair of missing atoms through the AutoDock tool suite [87].

For the docking calculations, all rotatable bonds of the ligands were free, while the target protein was set to a rigid conformation. Considering that all ligands presented many rotatable bonds and elevated flexibility, the AutoDock Vina algorithm was used due to its adequacy for dealing with ligands with high degrees of flexibility [88], with the support of mathematical machine learning [88] techniques to compose its scoring function.

The search box was defined through the observation of the 3QRJ crystal structure, which is in complex with rebastinib (DCC-2036) [46]. The coordinates of this crystallized ligand were used to calculate a centroid, the origin point for the docking search, whose coordinates were set to $[eqn]$ , $[eqn]$ , and $[eqn]$ while measuring 22.25 × 23.00 × 25.25 Å in the x, y, and z axes, respectively. The exhaustiveness level of the search was set to 32 in an energy range of 3.0 kcal·mol^−1^ while searching for 1000 binding modes. These settings were validated through a redocking calculation using the crystallized DCC-2036 ligand, which provided a docking model with an RMSD value of 0.1438 Å, which is within the acceptable deviation range [55,89,90,91]. Visual analysis of the resulting docking conformations was performed using the Discovery Studio software package and the VMD suite [92,93]. The binding conformations obtained through the docking studies were used as the starting geometries for the ONIOM calculations.

3.2. ONIOM Calculations

The complexes comprising the BCR-ABL enzyme and the docked models of ponatinib and imatinib were then subjected to geometric optimizations using the two-layer ONIOM model [94] in order to find the equilibrium geometry. In each system, the Lys285, Glu286, Ala287, Ile315, Glu316, Phe317, Met318, Asp381, Phe382, and Gly383 residues were included in the high layer, along with the docked ligands, resulting in three different systems as follows:

• 1.The T315I isoform of the BCR-ABL enzyme in the absence of any ligands;
• 2.BCR-ABL(T315I) complexed with rebastinib;
• 3.BCR-ABL(T315I) complexed with ponatinib.

The Glu286, Ile315, Met318, and Asp381 residues are the key targets of modeling studies, as they are the most important in the inhibition of the BCR-ABL enzyme [34,56,85]. The electronic properties of the high layer were calculated according to the Density Functional Theory (DFT), using the B3LYP functional, along with the 6-31G(d) basis set [95,96]. The use of dispersion-corrected functionals is important for the strongly attractive, non-covalent interactions [97,98,99,100]. Nevertheless, the B3LYP hybrid functional is widely used in drug design studies due to its reasonable computational cost and reliable structural predictions [101,102,103,104,105]. The B3LYP results show reasonable agreement with respect to both the energy and geometries of hydrogen-bonded structures [105]. The hybrid approach of B3LYP and AMBER applying electronic embedding approximations has been used in the literature to obtain non-covalent results for the enzyme-inhibitor complex [106,107,108,109,110,111,112]. The optimized ONIOM and docking geometries of rebastinib and ponatinib are available in the Supplemental Materials.

The remaining residues and crystallographic waters were included in the low layer, whose properties were calculated through the AMBER force fields [113]. AMBER force-field studies are primarily based on this choice of basis set and functional: B3LYP and 6-31G(d) [114,115]. Despite the modest basis set used, the accuracy of the geometry for this basis set is reliable due to the large size and complexity of the studied system. Nevertheless, diffuse functions are required to yield appropriate energies. In this way, the possible size for the high-level ONIOM region is limited with large basis sets [116]. In contrast, using a large model adequate to improve the description of $[eqn]$ interactions, B3LYP with a small basis set probably performs comparatively well for optimization geometry with less computational effort. However, it was found that the ONIOM approach using B3LYP/6-31G(d) improves the relative energies without achieving chemistry accuracy (less than 1 kcal·mol^−1^) [116,117,118,119,120]. For biological systems, the receptor binding energies were correctly assessed using ONIOM with B3LYP/6-31G(d), including BSSE [70,121]. Furthermore, based on B3LYP/6-31G(d) basis set under the ONIOM approach, the smallest median absolute deviation (MAD) was found for the proton and electron affinities of 1.29 kcal·mol^−1^ and 0.11 eV, respectively [122]. The inclusion of water molecules characterized an explicit solvation model for the calculations, which were performed using Gaussian16 software [123]. The calculated energy of the systems was obtained through Equation (1).

$$[eqn]$$

3.3. NCI, BCP, and NBO Analysis

Along with the calculated energies, the distances and angles of the hydrogen bonds formed with Glu286 and Asp381 were evaluated. The attractive charge interaction with the Asp381 side chain, when present, was also tracked. A qualitative analysis of Non-Covalent Interactions (NCIs) between the residues and ligands in the high-layer ONIOM-optimized structure was also performed using the Multiwfn 3.8 software package [124]. This analysis was complemented by a quantitative study of relevant Bond-Critical Points (BCPs) according to the Quantum Theory of Atoms In Molecules (QTAIM) [125] method using the AIMAll (Version 16.01.09) software package [126], which allowed for the evaluation of each interaction, mainly regarding the hydrogen bonds formed in the high layers of the complexes.

A Bond-Critical Point (BCP) can be characterized as a point in the electronic density space ( $[eqn]$ ) where all the first derivatives of the density gradient ( $[eqn]$ ) are equal to zero [127,128]. The study of these points offers quantitative data for highly rigorous theoretical calculations, such as the strengths of intermolecular interactions [127].

Typically, the electron density at a bond-critical point ( $[eqn]$ ) exceeds 1.35 e·Å^−3^ in covalent bonds, whereas it is less than 0.67 e·Å^−3^ in closed-shell interactions such as intermolecular interactions. These include a range of forces, from van der Waals interactions to hydrogen bonds and ionic interactions [127,128].

Additional insights can be gained by examining the Laplacian of the electron density at a given point ( $[eqn]$ ). Negative values of the Laplacian indicate a concentration of electron density along the bond length, while positive values signify electron density concentration near the atomic nuclei. Consequently, in covalent bonds, $[eqn]$ , whereas in intermolecular interactions, $[eqn]$ . For example, the typical Laplacian value for a hydrogen bond is approximately $[eqn]$ e·Å^−5^ [127,128]. Additionally, the Laplacian can be used to determine energy densities at the BCP. The Cremer–Kraka energy density (H^c^) is useful to indicate a non-covalent closed-shell interaction (electrostatic/dispersion type) for positive values, while H^c^ < 0 indicates a covalent interaction [129,130,131,132,133].

Finally, measurements of the donor–acceptor interactions were obtained through second-order perturbative analysis of Natural Bond Orbitals (NBOs). Natural Orbitals (NOs) are so-called due to their property of spontaneously emerging from the wave function as ideal descriptors. In the context of this study, the most important types of NBO are Lone Pair (LP), Bonding (BD) and Antibonding (BD*) NBOs. The analysis of NBOs can describe the results of orbital interactions, which are profoundly related to the stabilization energies (E(2)) of molecules and complexes through intermolecular interactions (Figure 3).

4. Conclusions

The present study was performed to provide a deeper understanding of the energetics involved in the binding modes of ponatinib and rebastinib regarding the mutant BCR-ABL enzyme according to refined QM/MM calculations. The analysis of docking conformations using AutoDock Vina through the lens of ONIOM provided important insights into the behavior of rebastinib and ponatinib when bound to the T315I isoform of the BCR-ABL enzyme. The aspects may be related to the difference between an effective and an ineffective inhibition model. Ponatinib was shown to form additional stronger and synergic hydrogen bonds, as well as to conserve them throughout the optimization of the docking geometry without significant deviations. QTAIM results and NBO analysis revealed the high indices involved in pertinent intermolecular interactions. The electronic properties—specifically, the LUMO−HOMO gap—of the BCR-ABL enzyme bound to ponatinib also showed indications of significant stabilization in the context of chemical reactivity, as suggested by the higher value calculated for the aforementioned gap. This decrease results in the inhibition of catalytic activity and enzyme activation due to more intense energetic barriers. Although rebastinib is described as a capable binder regarding the ATP-binding site of the enzyme, it was shown to do so to a lesser extent when compared to ponatinib, which can indicate that the former would be an inferior competitive inhibitor when compared to the latter. This behavior was expected, since ponatinib was specifically designed as a competitive inhibitor to the BCR-ABL enzyme, while rebastinib was originally aimed at the Tie2 family of enzymes and is described to act on the BCR-ABL enzyme not only through competitive inhibition.