Thermochemical and Quantum Descriptor Calculations for Gaining Insight into Ricin Toxin A (RTA) Inhibitors.

In this work, we performed a study to assess the interactions between the ricin toxin A (RTA) subunit of ricin and some of its inhibitors using modern semiempirical quantum chemistry and ONIOM quantum mechanics/molecular mechanics (QM/MM) methods. Two approaches were followed (calculation of binding enthalpies, ΔH bind, and reactivity quantum chemical descriptors) and compared with the respective half-maximal inhibitory concentration (IC50) experimental data, to gain insight into RTA inhibitors and verify which quantum chemical method would better describe RTA-ligand interactions. The geometries for all RTA-ligand complexes were obtained after running classical molecular dynamics simulations in aqueous media. We found that single-point energy calculations of ΔH bind with the PM6-DH+, PM6-D3H4, and PM7 semiempirical methods and ONIOM QM/MM presented a good correlation with the IC50 data. We also observed, however, that the correlation decreased significantly when we calculated ΔH bind after full-atom geometry optimization with all semiempirical methods. Based on the results from reactivity descriptors calculations for the cases studied, we noted that both types of interactions, molecular overlap and electrostatic interactions, play significant roles in the overall affinity of these ligands for the RTA binding pocket.

Introduction

Ricin is a cytotoxic protein produced in the seed of the castor bean plant (Ricinus communis), and it belongs to the ribosome-inactivating protein family (type-2) and is one of the most potent biological toxins known.[1] This protein consists of two subunits joined together by a disulfide bond, namely, ricin toxin A (RTA) and ricin toxin B (RTB).[2−4] RTA (267 residues) is an N-glycosidase that inactivates eukaryotic ribosomes via depurination of a specific adenine located in the sarcin-ricin loop (SRL) motif of the 28S rRNA subunit.[5] RTB (262 residues) is a lectin that mediates the uptake of holoricin into cells via recognition of galactose and N-6-acetylgalactosamine.[6] Once internalized, holoricin undergoes vesicular retrograde transport until reaching the endoplasmic reticulum lumen;[7] then, an isomerase reduces the disulfide bond between the subunits, and RTA is translocated to the cytosol and subsequently efficiently attacks ribosomes.[8,9]

There is a global concern by world authorities regarding ricin toxicity due to its potential use as a chemical weapon, mainly by terrorist and activist groups. In addition, the absence of countermeasures to ricin poisoning further contributes to this concern. In this way, theoretical methods, such as molecular docking, have guided research groups to find antagonist scaffolds for the RTA active site. However, this search has not been a trivial task because the active site of RTA and its surroundings are largely polar, imposing polarity constraints that in turn are not accounted for by the simplest theoretical methods.

To date, theoretical and computational chemistries have played a key role in studying biological and/or biochemical systems, providing proper direction toward drug discovery.[10] Knowledge of the intermolecular interactions of protein–ligand systems is an important feature for the development of new drugs.[11] In this way, some approaches, e.g., molecular docking, have been used to predict the bioactive pose of ligands in the active sites of biological targets with therapeutic interest;[12] however, it is quite ineffective to predict the relative binding affinity and rank ligands according to experimental data.[10,13,14]

Currently, the estimation of the binding free energy (ΔGbind) and binding enthalpy (ΔHbind) of a ligand in a protein–ligand complex is a major challenge in computational chemistry.[10] To tackle this problem, several groups have been using computational chemistry methods to predict better energy profiles.[10,11,15−33] Among these approaches are those entirely based on classical force fields, such as MM-G(P)BSA,[34,35] LIE,[36] SMD,[37,38] and FEP,[39] as well as their versions involving hybrid potentials (quantum mechanics/molecular mechanics (QM/MM)) or those totally calculated by quantum mechanics (QM), such as QM-MM-G(P)BSA,[40,41] QM-MM-LIE,[23,42,43] QM-FEP,[25,44] and QM-SMD.[45,46] All of these approaches require protein–ligand flexibility, either at their end-points or at the free-energy surface, resulting in high computational costs. Therefore, the use of such approaches is limited to a small set of ligands and small- and medium-sized complexes.

An alternative that requires lower computational cost is using theoretical methods for a single structure, in which the flexibility of the receptor (R), ligand (L), and receptor–ligand (R–L) complex can be partially accounted for by carrying out geometry optimization for R, L, and R–L. In this scenario, if the interest is to increase the accuracy of the receptor–ligand-binding energy predictions, the use of QM methods is mandatory. There are two approaches for calculating the binding energy of a ligand in a biological target using QM methods. The first approach considers a selection of the R–L complex atoms (QM cluster), and the second considers all atoms of the R–L complex in the QM calculation.

In the QM cluster strategy, a representative set of residues (ranging between 20 and 200 atoms) is selected and cut off from the active site. This set is studied separately by applying QM methods either in vacuum or in a continuous solvent model.[47,48] The advantage of this strategy is that few atoms are considered in the calculation, which allows us to assess a large set of ligands for the same active site. However, the disadvantage of this approach is that removing important residues during QM region selection can lead to inaccurate results.[49,50] In the strategy that considers all R–L atoms for QM calculations, it is possible to both hasten the computation of thermochemical properties and enable exploration of large ligand databases using linear scaling techniques coupled with graphics processing unit (GPU)-type accelerators.[51,52]

In the literature, some studies have performed molecular modeling and simulation involving ricin. Olson[53] applied molecular dynamics (MD) methods to analyze the structural and energetic aspects of three polynucleotide ligands (rRNA substrate analogues) that bound to the RTA active site. The results of the interaction energies showed that the overall binding is dominated by nonspecific interactions that occur in the region with a high degree of protein basicity through specific arginine contacts. The simulations of the three R–L complexes, as well as their comparison with experimental data, allowed a better understanding of the interaction of RTA with rRNA.

In another report, Olson and Cuff[54] expanded the previous study[53] by analyzing the free-energy determinants for the formation of RTA complexes with the rRNA substrate and several small ligands. The authors found that the absolute free energies of formation obtained for the RTA–RNA complex, as well as for several protein mutants, presented good agreement with the experimental data. In addition, it was observed that the terms of free energies presented unfavorable electrostatic contributions that were balanced by the favorable nonspecific hydrophobic effect, with free energies similar to those of protein–protein complexes. The individual components (by amino acid residue) of the binding free energy of the RTA–RNA complex revealed highly relevant electrostatic interactions arising from the charge–charge complementarity of the interfacial arginines with the RNA phosphate backbone. In addition, it has been observed that the hydrophobic complementarity of the domain is exerted by the base interactions of the GAGA loop structure.

Yan and co-workers[55] conducted studies on the interactions of small rings with the RTA active site to better understand how ricin recognizes adenine rings. The geometries and interaction energies were calculated using MM methods for some complexes between the RTA active site and tautomeric modifications of adenine, formycin, guanine, and pterin. The results indicated that the interaction energies between the pterin ring and RTA are stronger than those of formycin with RTA. It has also been found that formycin binds more strongly to RTA than adenine. This information presented good agreement with the experimental data. In addition, the results of experimental and molecular modeling work suggest that the binding site of ricin is quite rigid and can recognize only a small range of adenine-like rings.

In a recent study, Chaves and co-workers[12] carried out redocking of six known RTA inhibitors and then performed steered molecular dynamics (SMD) simulations to assess the relative binding affinity of these ligands. The molecular docking approach was able to predict the bioactive pose of the ligands; however, the score function was unable to rank these inhibitors according to experimental data. Steered MD simulation was used to decouple these ligands from the binding pocket, and the force profiles were estimated and presented a strong correlation with the experimental data (R2 > 0.9).

As a matter of fact, there are few molecular modeling or simulation studies about ricin, and in our searches, we found no study applying quantum chemical methods for whole RTA–ligand complexes. Thus, in this work, we performed calculations of ΔHbind and quantum descriptors between RTA and some of its inhibitors using semiempirical quantum chemical and QM/MM ONIOM methods to compare the results with experimental data (half-maximal inhibitory concentration (IC50)) and obtain local interaction information between the RTA residues and RTA inhibitors.

Methods

We retrieved the following RTA–ligand complex structures from the Protein Data Bank (PDB): 4HUP (RTA–19M),[56]4HUO (RTA–RS8),[56]4ESI (RTA–0RB),[57]4MX1 (RTA–1MX),[58]3PX8 (RTA–JP2),[59] and 3PX9 (RTA–JP3).[59]Figure presents the chemical structures of the six ligands of RTA studied in this work.

Figure 1

Structures of the six ligands of RTA studied in this work. The three-letter code corresponds to the PDB ID codes for these ligands.

All RTA–ligand complexes were relaxed and then equilibrated using molecular dynamics simulations. We set the MD simulations according to the settings used in ref (12). The last frame of each MD simulation was used to calculate binding enthalpies using semiempirical quantum mechanical (SQM) methods considering two scenarios: (i) a direct single-point energy calculation with RM1,[60] PM6,[61] PM6-DH+,[62] PM6-D3H4, and PM7[63] and (ii) after full-atom geometry optimization with PM6-DH+, PM6-D3H4, and PM7.

For all semiempirical calculations, we used the linear scaling algorithm MOZYME,[64] available in the MOPAC2016 package.[65] For the SCF convergence criteria, we used standard settings and a cutoff radius of 10 Å for the MOZYME algorithm. For full-atom geometry optimizations, we used the limited-memory BFGS algorithm considering a norm of gradient of 5.0 kcal·mol–1·Å–1 as stop criteria for complexes and noncomplexed proteins and 0.01 kcal·mol–1·Å–1 for the free ligands, where no restrictions of movement were considered for the atoms. For all calculations (single-point and geometry optimization calculations), we considered the conductor-like screening model (COSMO) implicit solvent field with a relative permittivity of 78.4 and an effective solvent molecule radius of 1.3 Å. The binding enthalpies for the ligand–RTA complexes were calculated according to eq . For all noncomplexed proteins and R–L complexes, we assumed a total charge of +2 e.The RTA has 4200 atoms and the ligands 19M, RS8, 0RB, 1MX, JP2, and JP3 have 67, 47, 30, 43, 20, and 31 atoms, respectively. The summation between the number of atoms for RTA and each ligand gives the number of atoms for the respective RTA–ligand complex.

For the sake of comparison and to test the performance of semiempirical methods, we also carried out single-point calculations considering hybrid QM/MM ONIOM method[66−68] using the same geometries from DM for the six RTA–ligand complexes. We used Gaussian 09 program[69] with the hybrid GGA B3LYP functional[70,71] and the hybrid GGA ωB97X-D functional,[72] which includes DFT-D2[73] dispersion correction. In both cases, we used the 6-31+G(d) basis set[74−76] for the QM region in both RTA protein and RTA–ligand complexes. UFF universal force field[77] was used for the remainder of the protein, which is denoted as the MM part. As requested, hydrogen atoms were used as link atoms connecting these two parts. We used all default criteria for ONIOM QM/MM calculations in Gaussian 09 program.

The QM part of the RTA protein includes the following residues: Glu-177, Arg-180 (important for enzyme catalysis), Tyr-80, Val-81, Gly-121, and Tyr-123 (important for attachment and recognition of the ligand at the active site).[78] In addition, we include residues Asn-122, Ser-176, Asn-209, Gly-212, and Arg-213 because they are about 3.0 Å apart from any of the six ligands, thereby forming hydrogen-bond interactions, producing a model with 189 atoms, 716 electrons, and +1 charge. For the RTA–ligand complexes, we also included the respective ligands in the QM part. Therefore, models were generated with (256 atoms, 1008 electrons and +1 charge), (236 atoms, 930 electrons and +1 charge), (219 atoms, 864 electrons and +1 charge), (232 atoms, 908 electrons and +1 charge), (209 atoms, 822 electrons and +1 charge), and (220 atoms, 864 electrons and +1 charge) for RTA–19M, RTA–RS8, RTA–0RB, RTA–1MX, RTA–JP2, and RTA–JP3 complexes, respectively. The binding energies for the ligand–RTA complexes were calculated according to eq .In Figure , we present the geometry for the RTA–19M complex with emphasis on the residues of the active site and the ligand (in yellow).

Figure 2

(a) Geometry of RTA–19M complex used in the ONIOM QM/MM calculation. (b) Zoomed view of the active site showing the 11 residues and 19M ligand (in yellow).

In Figure , we present QM models used for all complexes. We considered all atoms for residues Tyr-80, Val-81, Gly-121, Asn-122, Tyr-123, Asn-209, Gly-212, and Arg-213 because there are interactions between atoms of protein backbone and the ligands for such residues. For residues Ser-176, Glu-177, and Arg-180, we consider only side chains since intermolecular interactions with ligands occur only in this region of amino acids.

Figure 3

QM models for ONIOM QM/MM calculations for complexes: (a) RTA–19M, (b) RTA–RS8, (c) RTA–0RB, (d) RTA–1MX, (e) RTA–JP2, and (f) RTA–JP3.

Reactivity Descriptor Calculations

The most successful quantum descriptors come from conceptual density functional theory (CDFT), which through the mathematical development of density functional theory has given valid quantitative definitions for well-known chemical concepts,[79] such as hardness and softness. From this theory, the main interactions between two molecules are summarized in two types of processes: mutual polarization followed by molecular orbital overlap interactions and Coulombic forces.[80] The propensity of these effects in molecules can be traced locally using the Fukui function for the first type of interaction and local hardness for the second.[81]

To locally represent the molecular orbital overlap interactions, we used the Fukui function.[81] This descriptor was further divided into two functions, namely, the left Fukui function (f–), specific for electrophilic attack susceptibility (which for the frozen orbital approximation is equal to the highest-energy occupied molecular orbital (HOMO) density[80]), and the right Fukui function (f+), specific for nucleophilic attack susceptibility (which is equal to the lowest-energy unoccupied molecular orbital (LUMO) density).

A convenient representation of the Fukui functions is where the values are assigned for each k atom center in the molecule. The f– defined in eq is the sum of the squared ν atomic orbital (AO) coefficients that comprise the HOMO and belong to the kth atom, plus the sum of the product between the coefficients of indices ν, μ, and the corresponding overlap matrix element Sμν.[82] An equivalent definition for the f+ is presented in eq , using the LUMO instead of the HOMO.In this study, the overall effects tallied in these two Fukui functions were used in combination via the average Fukui function, defined in eq .These quantum chemical descriptors have been identified as useful for determining the toxicity and biological activity of several potential ligands.[83,84] For the enzymes, these descriptors have been used to determine reactivity sites, as score functions in molecular docking, and to search protein native structures[85] and find key structures in reaction path simulations.[86,87]

The electronic structure of protein systems presents several relevant molecular orbitals for chemical interactions with the electronic energy near the HOMO and LUMO.[84,86] Thus, in this study, we computed the Fukui functions not only by using the HOMO and LUMO orbitals but also by considering all of the molecular orbitals from a range of 3 eV from the HOMO and LUMO. As Fukushima and co-workers showed in their work, this value can vary from 1 to 5 eV depending on the system under study.[88]

We used the local hardness (H(k)) as a quantum descriptor to compute the Coulombic force interactions. This particular definition is based on the electron–electron contribution of the molecular electrostatic potential, as shown in eq ;[89] the calculation of local hardness at the kth atomic center is defined as a sum of the left Fukui function for each lth atomic center divided by their Euclidean distance R. These theoretical quantities were calculated using the PRIMoRDiA software.[90]

Results and Discussion

In Figure , we present the ΔHbind results of all of the RTA–ligand complexes studied in this work, considering both strategies are carried out: (i) direct single-point energy calculation and (ii) calculation of ΔHbind after full-atom geometry optimization of R, L, and R–L.

Figure 4

Binding enthalpies, ΔHbind, for the complexes (RTA–0RB, RTA–1MX, RTA–19M, RTA–JP2, RTA–JP3, and RTA–RS8). In (a), the results are depicted from single-point calculations, and in (b), the results are after full-atom geometry optimization of R, L, and R–L.

Looking at the results for ΔHbind calculated using the single-point strategy (Figure a), we can see that PM6-DH+, PM6-D3H4, and PM7 predicted negative ΔHbind values for all studied RTA–ligand complexes. The range of ΔHbind values is also consistent with the expected results from ligand–enzyme thermochemical assays, approximately a dozen of kcal·mol–1.[91]

PM6 and RM1 showed low performance, presenting positive values for ΔHbind for all RTA–ligand complexes. By comparing the performance among the PM6, PM6-DH+, and PM6-D3H4 methods, we verified that there are differences in binding enthalpies when corrections for dispersion and hydrogen bonding are considered. These results are consistent with recent studies suggesting that the quality of semiempirical methods presents significant improvements when such corrections are considered.[16,18,92] These studies also indicate that SQM methods incorporating dispersion and hydrogen bonding yield results similar to those of D-type corrections for density functional theory (DFT) methods.

Figure b shows the results of ΔHbind after full-atom geometry optimization for all complexes using the PM6-DH+, PM6-D3H4, and PM7 methods. We verified that the methods followed the same tendency of the single-point calculations, i.e., the calculated binding enthalpies were negative for all RTA–ligand complexes. The only exception was the RTA–JP3 ligand, which presented a positive ΔHbind value with the PM6-D3H4 method. This result suggests that the PM6-DH+, PM6-D3H4, and PM7 methods are able to describe, at least qualitatively, the binding enthalpies of the studied systems. We observed that the full-atom geometry optimization showed a tendency to decrease the ΔHbind value of the complexes. From the 18 ΔHbind calculations performed (optimization of six complexes with three distinct semiempirical methods: PM6-DH+, PM6-D3H4, and PM7), 12 presented a decreased ΔHbind value compared to those obtained by single-point calculations. Exceptions occurred for the RTA–0RB (ΔHbind = −60.50 kcal·mol–1), RTA–1MX (ΔHbind = −34.41 kcal·mol–1) and RTA–JP2 (ΔHbind = −33.23 kcal·mol–1) complexes with the PM7 method, RTA–JP3 (ΔHbind = 18.35 kcal·mol–1) and RTA–RS8 (ΔHbind = −65.33 kcal·mol–1) complexes with the PM6-D3H4 method and the RTA–19M (ΔHbind = −20.97 kcal·mol–1) complex with the PM6-DH+ method.

In Table , we present the IC50[56−59] and ΔHbind values for each RTA–ligand complex evaluated in this work.

Table 1

Binding Enthalpies (kcal·mol–1) for the Six RTA–Ligand Complexes Evaluated in This Study

		ΔHbind via single-point calculations					ΔHbind via geometry optimizations
ligand	IC50 (μM)	PM6	PM6-DH+	PM6-D3H4	PM7	RM1	PM7	PM6-DH+	PM6-D3H4
19M	15 (1)	9.58 (3)	–54.37 (1)	–66.47 (2)	–90.76 (1)	41.27 (5)	–156.11 (1)	–20.97 (6)	–86.38 (1)
RS8	20 (2)	6.10 (2)	–47.31 (2)	–68.23 (1)	–75.59 (2)	41.12 (4)	–93.33 (2)	–64.95 (3)	–65.33 (4)
0RB	70 (3)	5.55 (1)	–42.88 (3)	–65.53 (3)	–70.5 (3)	32.74 (3)	–60.5 (3)	–78.21 (2)	–80.12 (2)
1MX	209 (4)	26.28 (6)	–18.06 (5)	–33.01 (5)	–53.65 (4)	45.24 (6)	–34.41 (5)	–92.76 (1)	–59.54 (5)
JP2	230 (5)	9.71 (4)	–28.5 (4)	–39.53 (4)	–52.39 (5)	–29.51 (1)	–33.23 (6)	–56.1 (4)	–79.94 (3)
JP3	380 (6)	20.54 (5)	–8.12 (6)	–5.16 (6)	–24.89 (6)	24.78 (2)	–54.16 (4)	–45.87 (5)	18.35 (6)

When comparing the ΔHbind results obtained by single-point calculations using the PM6, PM6-DH+, PM7, and RM1 methods with IC50 experimental data for the ligands (19M, RS8, 0RB, 1MX, JP2, and JP3),[12] we found that the PM6 method showed correlation of 0.688 and RM1 method showed no correlation with the IC50. On the other hand, the PM6-DH+ and PM7 methods were able to identify the 19M ligand, which has the lowest IC50 value (15 μM), as the best ligand.

The PM6-D3H4 method switched the rank order of the 19M and RS8 ligands (ΔHbind = −66.47 and −68.23 kcal·mol–1, respectively), but this may have occurred because the IC50 values of these two ligands are very close (15 and 20 μM) so that the PM6-D3H4 method was not sensitive enough to rank the best ligand for small IC50 variations. The PM6-DH+ method also showed inversions of binding enthalpy values between the 1MX (ΔHbind = −18.06 kcal·mol–1) and JP2 (ΔHbind = −28.50 kcal·mol–1) methods. The PM7 method was able to correctly rank all of the ligands without inversions, being more sensitive to the small variations of IC50.

When we performed the statistical treatment between the IC50 and ΔHbind data obtained by single-point calculations, we found R values of 0.962, 0.975, and 0.984 for the PM6-DH+, PM7, and PM6-D3H4 methods, respectively. Thus, the PM6-D3H4 method presented the best correlation with the IC50 data. So, it is worth mentioning that the results for PM6-DH+ and PM6-D3H4 are quite satisfactory. In Figure , we present the correlation graphs between the IC50 and ΔHbind data obtained by single-point calculations when the PM6-DH+, PM6-D3H4, and PM7 methods were used.

Figure 5

Correlation graphs between IC50 and ΔHbind data obtained by single-point calculations for the (a) PM6-DH+, (b) PM6-D3H4, and (c) PM7 methods.

When analyzing the binding enthalpy data with full-atom geometry optimization (Table ), we verified that the PM6-DH+ method, although presenting good results with single-point calculations, did not show the same performance. This method was not able to identify 19M as the best ligand or present a good correlation with the IC50 data.

The PM7 method could rank the best ligands (from 1 to 3), but there were inversions in the ranking of the other ligands. Besides, there is a very marked distortion between the variation in the binding enthalpy results and IC50 values. The PM6-D3H4 method was able to identify the best ligand (19M) but showed inversions in the ranking between the RS8 and 0RB ligands and 1MX and JP2 ligands. However, compared to the PM7 method, the PM6-D3H4 method presented lower distortions between the ΔHbind and IC50 data.

When the statistical treatment was performed between the IC50 and ΔHbind data, we found R values of −0.085, 0.672, and 0.789 for the PM6-DH+, PM7, and PM6-D3H4 methods, respectively. Therefore, the PM6-D3H4 method also presented the best correlation with the IC50 when we performed a full-atom geometry optimization strategy. Moreover, we observed that in the full-atom geometry optimization, the correlation between the IC50 data and binding enthalpies decreased considerably.

Sulimov et al.[93] recalculated the energies of approximately 8000 best-ranked structures in the molecular docking with the MMFF94 force field through single-point calculations and ligand optimization using the PM7 method, considering an implicit COSMO solvent model. The authors observed that the energies obtained by single-point calculations greatly improved the ranking of structures close to the native state. However, when optimizing the ligands with the PM7 method in vacuum and recalculating their energies with the COSMO model, the authors observed a worsening of the results for several complexes compared to those with PM7 (single-point energy calculation). In accordance with those results, we observed more comprehensively that the geometry optimizations reduced the performance of all considered semiempirical methods.

Although the correlation for full-atom geometry optimizations is lower than that obtained through single-point energy calculations, we emphasize that if the molecular dynamics simulation is not carried out properly to find an equilibrated structure, then this optimization is necessary to obtain some correlation with experimental data. We have indicated PM7 and PM6-D3H4 as the best semiempirical methods for obtaining binding enthalpies and ligand rankings. We are not stating that single-point QM calculations on the end frame from MD simulation are better for ligand-binding energies than ones that use geometry optimization at the QM level. However, this was true for our study. Effects due to low sampling, entropy lacking, or cancelation of errors could be behind these results, but the investigation of these issues is outside the scope of this study. Anyhow, an excellent discussion about these issues can be found in the review by Ryde and Soderhjelm.[10]

In Table , we present the root-mean-square deviation (RMSD) data between the MD structures and those optimized with the PM6-DH+, PM6-D3H4, and PM7 methods.

Table 2

RMSD Results (in Å) between the End Frame of the MD Simulations and the Optimization of These Same End Frames by Semiempirical Quantum Methodsa

complex	PM6-DH+ (Å)	PM6-D3H4 (Å)	PM7 (Å)	crystallographic (Å)
RTA–19M	1.554	1.296	2.108	2.176
RTA–RS8	1.747	1.533	2.024	2.266
RTA–0RB	1.567	1.401	1.950	1.839
RTA–1MX	1.832	1.468	1.919	2.676
RTA–JP2	1.581	1.626	2.036	1.917
RTA–JP3	1.745	1.511	2.057	2.103

In the last column, we present RMSD results (in Å) between the end frame of the MD simulations and each crystallographic structure.

When analyzing the data in Table , we verified that the PM6-D3H4 method presented the lowest RMSD values, except for the RTA–JP2 complex, where the PM6-DH+ method presented a lower RMSD value. The PM7 method presented the highest RMSD values for all of the analyzed complexes. This observation suggests that the PM7 method is able to relax the initial structure more than the PM6-DH+ and PM6-D3H4 methods. Figure shows the superposition between the initial structure and the structures optimized via the PM7 method.

Figure 6

Superposition of the RTA–ligand optimized complexes: (a) RTA–0RB, (b) RTA–1MX, (c) RTA–19M, (d) RTA–JP2, (e) RTA–JP3, and (f) RTA–RS8. The optimized structures are represented in green, and initial structures are represented in yellow.

From the analysis of the RTA structures in the complexes, we verified that there were no significant changes after optimization, with conservation of the secondary structures. All R–L complexes presented RMSDs close to 2.0 Å, which is consistent with the resolution of the experimental data.[94]

When comparing the ligand poses via full-atom geometry optimization with structures from molecular dynamics, we observed that the 0RB ligand (Figure a) presented rotation of the triazole group, being approximately perpendicular to the preferred position of this group in the protein. The ligand 1MX (Figure b) presented small pterin group torsion and benzene ring variation. Due to its high conformational flexibility, the 19M ligand (Figure c) underwent a large modification after geometry optimization, mainly regarding torsions in the pterin moiety, benzene rings, and carboxylic acid group. The JP2 ligand (Figure d) exhibited rotation of its carboxylic acid group and small deviations in other groups of the structure. The JP3 ligand (Figure e) showed modification of its structure, with displacements in the positions of atoms along the same plane. The RS8 ligand (Figure f) exhibited a small variation in comparison to that observed in the last frame obtained from the molecular dynamics simulations.

In Table , we present the RMSD data between the ligands of complexes from the MD simulations and the ligands of the complexes optimized with the PM6-DH+, PM6-D3H4, and PM7 methods.

Table 3

RMSD Results (in Å) between the Ligands of the Complexes Optimized by Semiempirical Methods and the Ligands of the Complexes of the MD Simulationsa

complex	PM6-DH+ (Å)	PM6-D3H4 (Å)	PM7 (Å)	crystallographic (Å)
19M	0.675	0.657	1.530	6.005
RS8	0.808	1.291	0.598	5.819
0RB	0.761	0.755	1.062	2.856
1MX	0.932	0.477	0.812	4.318
JP2	0.889	0.500	0.909	1.957
JP3	0.581	0.510	0.612	2.597

In the last column, we present RMSD results (in Å) between the ligands of the MD simulations and the ligands from crystallographic structures.

When analyzing the data of Table , we verified that the poses of the ligands presented several variations after the geometry optimization for all R–L complexes. For the PM6-DH+ method, the 1MX ligand presented the highest variation (RMSD = 0.932 Å), and the JP3 ligand presented the lowest variation for this method (RMSD = 0.581 Å). For the PM6-D3H4 method, it was observed that the RS8 and 1MX ligands presented the highest and lowest variations, respectively, with RMSDs equal to 1.291 and 0.477 Å. For the PM7 method, the 19M ligand showed the highest variation (RMSD = 1.530 Å), and the lowest variation for this method was for the RS8 ligand (RMSD = 0.598 Å). We observed that the PM7 method showed higher RMSD values for four of the six ligands, and the two exceptions were the RS8 and 1MX ligands.

In Table , we present the IC50[56−59] and ΔEbind values for each RTA–ligand complex calculated via ONIOM QM/MM with B3LYP and ωB97X-D functionals.

Table 4

Binding Energies, ΔEbind, (in hartrees) for the Six RTA–Ligand Complexes Calculated with ONIOM QM/MM Single-Point Calculations

complex	IC50 (μM)	ΔEbind/B3LYP (au)	ΔEbind/ωB97X-D (au)
19M	15 (1)	–0.136 (1)	–0.230 (1)
RS8	20 (2)	–0.119 (2)	–0.202 (2)
0RB	70 (3)	–0.106 (3)	–0.178 (3)
1MX	209 (4)	–0.062 (6)	–0.130 (4)
JP2	230 (5)	–0.070 (4)	–0.129 (5)
JP3	380 (6)	–0.057 (5)	–0.088 (6)

When analyzing the data in Table , we found that the single-point ONIOM QM/MM calculations with B3LYP functional were able to identify the three best ligands (19M, RS8, and 0RB). The only exception was the inversion between the ligands 1MX and JP2. The correlation between the IC50 and ΔEbind data presented an R-value of 0.929. Although the ONIOM QM/MM results with B3LYP functional showed good correlation with IC50, PM6-DH+, PM7, and PM6-D3H4 showed better results.

For the sake of testing, we also carried out single-point ONIOM QM/MM calculations with the B3LYP functional considering a smaller QM region for each complex. In this test, we considered in the QM part only the six most important residues of RTA plus the ligand. When we did this the correlation decreases dramatically to 0.693.

When we changed the B3LYP functional to the ωB97X-D functional, which includes DFT-D2[73] dispersion correction, we observe that the R-value increased from 0.929 to 0.972, showing a slight improvement in the correlation with IC50. In addition, calculated results obtained with this functional, all ligands were ranked correctly according to ΔEbind and IC50 values.

Considering these results and the low computational cost of semiempirical methods, we can state that methods PM7, PM6-DH+, and PM6-D3H4 showed better performance than the ONIOM QM/MM calculation with B3LYP functional for RTA complexes studied in this paper. When using the ωB97X-D functional, the QM/MM method presented R-value and performance similar to the PM7 method but with a higher computational cost. The advantage of ONIOM QM/MM is that one can always improve the results by increasing the QM region and using larger basis sets, and/or taking solvent effects and dispersion corrections into account, but this ends up being computationally expensive.

As a final experiment with binding affinities for these RTA complexes, we carried out new calculations using the dataset of binding affinities produced in the study carried out by Chaves and his collaborators.[12]

Our inspiration was a study carried out by Gupta and collaborators, which achieved a significant improvement for the binding affinity predictions when they considered molecular descriptors in different molecular docking approaches for the inhibitors of microsomal prostaglandin E synthase-1.[95] In their study, molecular descriptors such as topological polar surface area (TPSA), partition coefficient (log P), volume (Vol), and the number of rotatable bonds (Nrtb) were used as docking scores. With this, the authors took into account desolvation penalties and conformational free-energy changes when a ligand binds to a protein. In their study, the re-score routine was considered as an empirical summation of normalized docking affinities. Similarly, we also incorporated a re-score routine to the binding affinity predictions provided by Vina software[96] for the same set of ligands (19M, RS8, 0RB, 1MX, JP2, and JP3). These binding affinity predictions were retrieved from another validation study carried out by our group; therefore, details about docking protocol can be reviewed in ref (12). In summary, the ligand molecular descriptors such as TPSA, Vol, and Nrtb were calculated using the web service Molinspiration (https://www.molinspiration.com). We also normalized the docking scores and the molecular descriptors to values between 0 and 1. Then, the following equation was used to perform re-score calculationsWe present in Figure the Pearson correlations between pIC50 vs Vina score and pIC50 vs re-score, respectively. In short, our results showed that the consideration of such molecular descriptors improves the relative binding affinity predictions provided by Vina software, with an increase from 0.500 to 0.814 in the Pearson correlation.

Figure 7

Pearson correlation between pIC50 vs Vina score and pIC50 vs re-score that is defined in eq .

Reactivity Descriptors

We calculated the reactivity descriptors for the protein–ligand complexes and used to perform a graphical analysis to theoretically characterize the most important interactions that occur between the inhibitors and the active site residues. As shown in Figure , the local hardness has higher values in the binding pocket for the RTA–19M complex, with its highest value at Arg-180 and Trp-211. For the RTA complex with the JP3 ligand, the electrostatic interactions are not placed within the ligand region. Significant Fukui function values were computed in the residues of the binding pocket for the three most active ligand complexes and with small values for the JP2 ligand. Figure depicts the Fukui function for the second and third most active cases. For the RTA–RS8 complex, the Fukui function localizes at one of the ligand rings of the pterin group and in Tyr-80. In the RTA–0RB complex, all atoms in the pterin group present high values of the given descriptor, as does Arg-180.

Figure 8

Calculated local hardness for RTA complexes with the 19M and JP3 ligands (orange). (A) Labels for the nearest residues (green) from 19M; (B) local hardness for the RTA–19M complex; (C) labels for the nearest residues (green) from JP3; and (D) local hardness for the RTA–JP3 complex.

Figure 9

Calculated Fukui function for RTA complexes with the RS8 and 0RB ligands (orange). (A) Labels for the nearest residues (green) from RS8; (B) Fukui function for the RTA–RS8 complex; (C) labels for the nearest residues (green) from 0RB; and (D) local hardness for the RTA–0RB complex.

In Figure , the Fukui function is presented for 1MX and JP2, showing the distribution of the descriptor at the pterin group and none at the closest residues. The molecular structure of these ligands is based on the pterin group,[58] the same group present in adenine that is hydrolyzed in the ribosome by the action of RTA catalysis.[5] In the complexes considered by calculations, the reactivity in this group is always indicated by the Fukui function.

Figure 10

Calculated Fukui function for RTA complexes with the 1MX and JP2 ligands (orange). (A) Labels for the nearest residues (green) from 1MX; (B) Fukui function for RTA–1MX complex; (C) labels for the nearest residues (green) from JP2; and (D) local hardness for the RTA–JP2 complex.

The models built from thermochemical properties correlations could rank the best inhibitors, but only reactivity descriptors can pose, which are the molecular structure features that new ligands must have to be more efficient.

Conclusions

The main objective of this work was to present a study where two computational approaches were applied (calculation of ΔHbind and reactivity quantum chemical descriptors) to theoretically gain insights into the interactions between the RTA subunit of ricin and some of its inhibitors.

In our study, we observed that single-point energy calculations of ΔHbind with the PM6-DH+, PM6-D3H4, and PM7 methods presented an excellent correlation with the IC50 data. In addition, although the PM7 method presented a slightly lower correlation than the PM6-D3H4 method, only the PM7 method was able to correctly rank all of the ligands analyzed. This result suggests that the PM7 method is more sensitive to small IC50 variations.

When comparing the IC50 data with the calculated ΔHbind values obtained after full-atom optimization with the PM6-DH+, PM6-D3H4, and PM7 methods, we find that the correlation decreases significantly. Although the correlation was greatly reduced with the optimization of the structures, the PM6-D3H4 and PM7 methods were able to correctly rank the best ligand. This result suggests that if the molecular dynamics is not carried out properly to obtain an equilibrated structure, it is necessary to carry out a full-atom geometry optimization to obtain some correlation with the experimental data.

We also observed that in geometry optimization, the structure adopts a minimum local conformation that contributes less to the preferential position of the ligand at the active site. On the other hand, the representative structure from molecular dynamics represents the preferred position of the ligand in the active site. Thus, we conclude that for the dataset studied, it is preferable to use equilibrated MD structures to perform single-point energy calculations to obtain ΔHbind.

ONIOM QM/MM single-point calculations with the B3LYP functional showed good correlation with IC50; however, the methods PM6-DH+, PM6-D3H4, and PM7 showed better performance. When using the ωB97X-D functional that includes DFT-D2[73] dispersion correction, the QM/MM method presented R-value and performance similar to the PM7 method. Besides, ONIOM QM/MM calculations for ΔHbind incurred a higher computational cost compared with semiempirical methods, about 250 times higher (B3LYP) and 1400 times higher (ωB97X-D).

In addition, we found for the cases studied, reactivity descriptors pointed out that both types of interactions, molecular overlap and electrostatic interactions, play an important role in the overall affinity of these ligands for the RTA binding pocket. From a relatively simple computational protocol, this approach opens the possibility to reveal additional information using structures treated with earlier simulations.