Study of the Differential Activity of Thrombin Inhibitors Using Docking, QSAR, Molecular Dynamics, and MM-GBSA

Non-peptidic thrombin inhibitors (TIs; 177 compounds) with diverse groups at motifs P1 (such as oxyguanidine, amidinohydrazone, amidine, amidinopiperidine), P2 (such as cyanofluorophenylacetamide, 2-(2-chloro-6-fluorophenyl)acetamide), and P3 (such as phenylethyl, arylsulfonate groups) were studied using molecular modeling to analyze their interactions with S1, S2, and S3 subsites of the thrombin binding site. Firstly, a protocol combining docking and three dimensional quantitative structure–activity relationship was performed. We described the orientations and preferred active conformations of the studied inhibitors, and derived a predictive CoMSIA model including steric, donor hydrogen bond, and acceptor hydrogen bond fields. Secondly, the dynamic behaviors of some selected TIs (compounds 26, 133, 147, 149, 162, and 177 in this manuscript) that contain different molecular features and different activities were analyzed by creating the solvated models and using molecular dynamics (MD) simulations. We used the conformational structures derived from MD to accomplish binding free energetic calculations using MM-GBSA. With this analysis, we theorized about the effect of van der Waals contacts, electrostatic interactions and solvation in the potency of TIs. In general, the contents reported in this article help to understand the physical and chemical characteristics of thrombin-inhibitor complexes.


Introduction
Thromboembolic diseases are among the principal causes of mortality in the world. Vein thrombosis can progress to pulmonary embolism. These disorders, identified with the term venous thromboembolism (VTE), affect several million people around the world [1]. VTE is the third leading cause of cardiovascular-related death, after myocardial infarction and stroke [2].
The central role of the serine protease thrombin in thrombosis and haemostasis makes it an attractive target for antithrombotic therapy [3]. Thrombin catalyzes the conversion of soluble fibrinogen to insoluble fibrin in the clotting cascade, and also acts on other substrates such as factor V, factor VIII, factor XI, and factor XIII. It is widely believed that an oral thrombin inhibitor (TI) could provide a new standard of care in anticoagulation therapy.
The discovery of small molecule TIs is an important goal for anti-thrombotic therapy [4]. In the last years, potent and selective inhibitors have been reported, such as pyridones, acetamides, oxyguanidines, aminopiperidines, amidines, and amidinohydrazones [5][6][7][8][9][10][11][12][13]. These sets have shown that subtle structural differences in compounds (due to the presence of similar scaffolds or the same scaffolds with different substituents) can lead to big differences in their thrombin inhibitory activities.
The knowledge of the relevant structural features that positively influence the activity of TIs is important for the design of potent compounds. Molecular modeling has demonstrated to be a powerful support to investigate bioactive compounds and their structure-activity relationship (SAR) with the main purpose of identifying the molecular features that contribute to a high bioactivity. These methods have been applied for studying TIs. Several quantitative structureactivity relationship (QSAR) models were reported using approaches such as classic QSAR [14], CoMFA/CoMSIA [15,16], topological descriptors [17], and artificial neural networks [18]. Other reports used docking and molecular dynamics (MD) simulations to study structural features of several TIs identified wit a high activity [19][20][21]. In general, these reports do not include an analysis of the interactions between inhibitors and different subsites of the thrombin binding site. An exception is the work of Nilsson et al. [22]. These authors designed a set of compounds to bind to the S 2 and S 3 subsites with no interactions in S 1 , and developed a classic QSAR model to study the chemical features that are important for the prediction of their binding constants. Other exception is the work of Bhunia et al. [23], which used three dimensional (3D) QSAR and MD simulations to profile structural determinants for the selectivity of representative diverse classes of thrombin-selective inhibitors.
In the current work, we applied some of the popular molecular modeling methods to study the interactions of TIs with S 1 , S 2 , and S 3 subsites of the thrombin binding site. We studied the orientations and SAR for 177 TIs by using a protocol that includes docking and the 3D-QSAR method CoMSIA. Additionally, we analyzed the dynamical behavior for some selected compounds (26,133,147,149,162, and 177) by using MD and free energy calculations. By means of a comparison of the selected systems, we gained further insight into the role played by different TI molecular constituents in the binding affinities due to interactions with different subsites in the thrombin active site.

Data set
The primary structures and activities of 177 TIs were taken from the literature [5][6][7][8][9][10][11][12][13]. Inhibitory activities were collected and transformed into log(10 3 /Ki) values. Ki values are in nM and represent the enzyme inhibition constants. TI structures are in Fig 1 and their biological activities used in this study are summarized in Table 1. The chemical structures were sketched using the molecular editor of Maestro 9.0 software suite [24].

Docking
The use of docking has extended in the last decades because successful applications to research in medicinal chemistry. Typically, docking programs are able to closely reproduce the binding orientations of the inhibitors in complexes for which crystallographic data are available [25][26][27]. Docking was performed using Glide [28], which is part of the Maestro 9.0 software suite [24]. Protein coordinates were extracted from the crystal structure of the thrombin-inhibitor complex with the code 1T4U in Protein Data Bank [6]; the protein structure was refined and completed using the Protein Wizard Preparation module also available in Maestro software. A grid box of 30Å × 30Å × 30Å was centered on the center of mass of the inhibitor in this crystal structure. The module LigPrep 2.4 [29] was used to assign ionization states, stereochemistries, and ring conformations of the ligands. Docking parameters were used as in previous works [30], Glide standard (SP) and extra-precision (XP) modes were used. The better docking poses for each ligand were examined according to their relative total energy scores. Among docking poses, the more energetically favorable conformation was selected by considering the total energy value.

CoMSIA calculations
Compound structures and inhibitory activities as log(10 3 /Ki) values were considered for CoM-SIA calculations. The relative alignment of the compounds to set field calculations was inside the binding site. For this, the 3D conformations previously obtained by docking simulations were used. QSAR modeling was performed using the Sybyl 7.3 software of Tripos [31]. The data set was divided into two sub data sets (143 and 34 compounds were in training and test sets respectively) for external validation process. The molecules contained in the training set were placed in a rectangular grid extended beyond 4 Å in each direction from the more external coordinates of each molecule. The interaction energies between a hypothetical atom (a sp 3 hybridized carbon atom with +1 charge) and all compounds were computed at the defined points, using a volume-dependent lattice with 2.0 Å grid spacing. Then, partial least squares (PLS) method was applied using standard Sybyl parameters, and considering an optimal number of components determined by optimization of leave-one-out (LOO) cross-validation Q 2 value (using SAMPLS [32] sampling method). Similarity was expressed in terms of steric Compounds 151-163

Molecular dynamics and MM-GBSA calculations
Complexes of the compounds 26, 133, 147, 149, 162, and 177 inside thrombin active site were studied by MD simulations. All simulations were performed using the NAMD software package [33]. The initial coordinates for the MD calculations were taken from the docking experiments. To mimic the aqueous environment, an equilibrated water box with sides of 100 Å, centered on the mass center of each inhibitor, was used to solvate the thrombin-ligand system. Protein and inhibitors were described using the optimized potential for liquid simulations CHARMM force field [34] and the CGenFF force field [35], respectively. In turn, the water molecules belonging to the solvent box were described using the flexible TIP3P potential [36,37]. Energy minimization was performed on the models using conjugate gradient method (20.000 steps) to reduce any close contacts. Afterwards, the system was further equilibrated during 2.0 ns. The production run consisted of an MD simulation of 5.0 ns. In all cases, we applied a constraint to the backbone amino acids, we used the NPT ensemble at 300 K, we set a time step of 1 fs to solve the equations of motion, and we set a switched cutoff distance of 9.0 Å.
The free energy calculations were accomplished using the MM-GBSA method [38]. This method combines molecular mechanics energy and implicit solvation models at a reasonable computational cost, leading to outstanding results in several biological systems in recent years [39][40][41]. We applied this method to 500 snapshots extracted from the 5.0 ns production MD trajectories (explicit TIP3P water molecules and ions were removed for this).
Protein-ligand binding free energy using MM-GBSA was calculated as the difference between the energy of the bound complex and the energy of the unbound protein and inhibitor compound. The method allows for free energy decomposition into contributions originating from different types of physico-chemical interactions. Specifically, the energy is calculated for the protein-ligand complex, the ligand, and the protein, and their energies were computed using the CHARMM force field with the generalized Born implicit solvent model, in order to calculate the averaged binding free energy (ΔG) according to the following equation: where ΔE MM includes ΔE internal (bond, angle, and dihedral energies), ΔE elect (electrostatic), and ΔE vdw (van der Waals) energies; ΔG solv is the electrostatic solvation energy (polar and nonpolar contributions). The polar contribution is calculated using the Generalized Born (GB) model, while the non-polar energy is estimated using the solvent accessible surface area (SASA). The conformational entropy change -TΔS can be computed by normal-mode analysis on a set of conformational snapshots taken from MD simulations, but many authors have been reported that the lack of the evaluation of the entropy is not critical for calculating the MM-GBSA (or MM-PBSA) free energies for similar systems [38,42,43]. In this work MM-GBSA calculations were also achieved in NAMD software [33]; the entropy term -TΔS was not calculated to reduce computational time.

Docking results
First, we applied docking methodology in order to reproduce the crystal structures of thrombin-ligand complexes for compounds 13, 24, 128, 151, and 165 (accession codes in PDB: 1T4U, 1T4V, 3C27, 2R2M, and 3LDX respectively). This initial test was used to assess the quality of the docking method to reproduce known structures. In Fig 2 is shown that the docked structures fitted in an acceptable way with available inhibitor X-ray crystal structures; all the inhibitors were adequately oriented. The values of the root mean square deviation (RMSD) for the docked structures with respect to the co-crystal inhibitor structures considering all heavy atoms were < 2.0 Å in almost all the cases analyzed. Considering that 2.0 Å is the threshold value that differentiates between correct and incorrect docking solutions [44], we can state that Glide found the correct binding mode of the ligands in four of the five cases analyzed. Despite compound 13 had the RMSD value above 2.0 Å, it had the expected orientation in the thrombin binding site. The high RMSD value for compound 13 was due to a small displacement of P 2 group, and a different orientation of the P 3 group of the docked conformation with respect to the conformation in the crystallographic structure. The remaining compounds were docked following the same docking protocol. The analysis of the docking poses obtained for all 177 ligands shows that all compounds adopt the same binding mode (Fig 3; mol2 files are in S1 File). This similar binding mode was expected since  The S 1 subsite is a large cavity bounded in one part by the residues contained in the sequence VSWGEGC (Val213-Ser214-Trp215-Gly216-Glu217-Gly218-Cys219). The backbone atoms of these residues and the side chain atoms of Val213 form a semicircle controlling the shape of the S 1 subsite (Fig 3B). The other part of the S 1 subsite is bounded by the residues contained in the sequence DACE (Asp189-Ala190-Cys191-Glu192), where Asp189 contributes with a negative charge, Cys191 forms a disulfide bond with Cys219 from the VSWGEGC motif, and the residues Ala190, Cys191, and Glu192 contribute with their backbone atoms ( Fig  3C). The hydroxyl of the residue Ser195 and the CH 2 of Gly226 also face the subsite S 1 . In general, all compounds kept the electrostatic interactions between a HB donor group in P 1 and the side chain carboxylate of Asp189. The backbone carbonyle groups of Gly218 and Ala190 are close to the site where interactions with Asp189 are formed; therefore they contribute with negative density to attract the P 1 group (according to crystal structures previously reported, Gly218 can form HB with the P 1 group). The S 1 subsite can bound linear or cyclic groups that can adopt several conformations. The residue Ser195 (serine from the catalytic triad) has the side chain hydroxyl facing the entrance of S 1 ; polar groups in this zone can interact with this residue (for instance, NH group in compounds 128, 151, and 165).
The S 2 subsite is delimited on one side by the residues Tyr61 and Trp64 from the YPPW insertion loop, and the His57 from the catalytic triad. On the other side S 2 is delimited by the residues Trp215 and Gly216 from the motif VSWGEGC (Fig 3D). The S 2 subsite was occupied by central aromatic group in all the docked compounds. The residues His57, Tyr61, and Trp64 form a small hole that is occupied by short substituents such as CH 3 , CN, or halogens. When bigger substituents (for instance OCH 3 in the low active compound 91) are present, the ligand central aromatic ring changes its orientation to place the bulky substituent out of the hole. The residue Trp215 exposes Cα and Cβ (CH and CH 2 groups) to the S 2 subsite; therefore it contributes to establish hydrophobic interactions in the portal between S 2 and S 1 . Finally, the residue Gly216 exposes backbone CO and NH groups, and several TIs form HBs with them (for instance, the highly active compounds 164-177 that contain the 3-aminopyrazin-2(1H)-one scaffold form HB with both Gly258 backbone groups).
The S 3 subsite is an ample pocket exposed to the solvent bounded in one part by the residues contained in the sequence WRENL (Trp95-Arg96-Glu97-Asn98-Leu99). The backbone atoms of these residues and the side chain atoms of Leu99 form a wall where three CO groups (from Trp95, Arg96, and Glu97) are facing the S 3 subsite (Fig 3E and 3F). Very close to this wall, the side chain hydroxyl group of Tyr61 faces to the S 3 subsite. The other part of the S 3 subsite is delimited by the hydrophobic chemical groups of the side chains of the residues Ile174, Trp215, and Glu217. Interestingly, S 3 subsite has a hydrophilic wall (backbone CO groups of Trp95, Arg96, and Glu97, and side chain hydroxyl group of Tyr61) and a hydrophobic wall (side chains of the residues Leu99, Ile174, Trp215, and Glu217). Docking results show that the S 3 subsite tolerates the presence of different aromatic groups. Very hydrophobic P 3 groups were located close to the hydrophobic wall and groups with some polarity were located close to the hydrophilic wall. Very large P 3 groups were also oriented towards the hydrophobic wall because a bigger available space in this area. For compounds that have a P 1 moiety with white loop, Ser195 is represented with yellow sticks. (D) Subsite S 2 , residues His57, Tyr61, and Trp64 are represented with violet spheres, and residues Trp215 and Gly216 from the VSWGEGC motif are represented with white spheres. (E, F) Subsite S 3 , residues of the WRENL motif are represented with brown spheres, residues Trp215 and Glu217 from the VSWGEGC motif are represented with white spheres, and Ile174 is represented with yellow spheres. In f several compounds that have a P 1 moiety with two branches are represented: compounds with two hydrophobic branches are represented with green sticks, and compounds with one of the branches containing a polar group are represented with pink sticks. doi:10.1371/journal.pone.0142774.g003 two branches (compounds 24-77) we found that branches are distributed according to polarity. If the two branches are hydrophobic, both groups are located close to the hydrophobic face (compounds in green in Fig 3E), but if one of the branches contains a polar group, it is located close to the hydrophilic face and the other branch is located close to the hydrophobic face (compounds in pink colour in Fig 3E).

CoMSIA results
CoMSIA models were carried out to provide information about the structural features affecting the thrombin inhibitory activity of the compounds under study. The use of conformations aligned using the binding site (obtained by docking) allows identifying the relevant pharmacophoric features required to best match in the binding site [45]. In this sense, the best model reported below accounts for the desired features that characterize the most potent TIs.
Firstly, the CoMSIA models were developed by including one field, and then, these fields were combined and the statistical quality of hybrid models was analyzed by considering Q 2 values [46]. According to this analysis, we detected that compounds 5, 18, 33, 90, and 140 presented large residuals in all models; in this sense, they were excluded as outliers. Outliers are those compounds that have unexpected biological activities and are unable to fit in a QSAR model [47].
The list of CoMSIA models using different fields, without considering outlier compounds, is presented in Table 2. All the combinations were tested but only the more relevant ones are reported in Table 2 (models including one field, models including two fields with Q 2 > 0.2, NC is the number of components from PLS analysis; R 2 is the square of the correlation coefficient; S is the standard deviation of the regression; F is the Fischer ratio; Q 2 and S cv are the correlation coefficient and standard deviation of the leave-one-out (LOO) cross-validation, respectively. The best model is indicated in boldface.
doi:10.1371/journal.pone.0142774.t002 models including three or more fields with Q 2 > 0.3. The internal predictability of the models was the criterion that was used to select the best QSAR model. We observed that CoMSIA models using one field were statistically unacceptable (Q 2 < 0.5). The analysis of the hybrid models yielded the model CoMSIA-SDA with the best Q 2 value of 0.543. This model combines steric, HB donor, and HB acceptor fields. We included other fields and we observed that it does not produce an improvement in the internal validation of the model CoMSIA-SDA, since models including more fields had lower Q 2 values ( Table 2). The model CoMSIA-SDA was derived by using eight components and showed contributions of the steric field of 18.6%, HB donor field of 35.0%, and the HB acceptor field of 46.4%. In addition, it explains 87.1% of the variance, has a low standard deviation (s = 0.322), and a high Fischer ratio (F = 109.20). The predictions of log(10 3 /Ki) values for the 138 TIs from the training set using model CoM-SIA-SDA are shown in Table 1. The correlation between the calculated and experimental values of log(10 3 /Ki) (from training and LOO cross-validation) is shown in Fig 4. Plots of the LOO cross-validation predictions reveal that the proposed model is able to discriminate between the most active and the less active compounds. Considering that the CoMSIA approach was applied to a big dataset, we consider that the value of Q 2 = 0.543 reflects that there is no redundancy in the training set since each member of the training set is important for the model [48]. We also used model CoMSIA-SDA to predict the TI inhibitory activities of the test set compounds. The values are given in Table 1 and correlation between the calculated and experimental values are also represented in Fig 4. This analysis reveals that the proposed model also predicted adequately all the compounds in the test set.
The contour plots of the CoMSIA steric, HB donor, and HB acceptor fields are presented in Fig 5 for the best model CoMSIA-SDA. To aid in visualization, compound 134 is displayed in the maps with a superposition of the contour plots on active-site residues. The representation of contour plots and the thrombin active site in the same 3D space helps in delineating the relationship between the relevant chemical properties of the more active compounds and the chemical groups of the residues in the active site.
The colored isopleths in the map represent the 3D space where the structural properties changes are related to the changes in thrombin inhibitory potency. Green and yellow isopleths in Fig 5A indicate regions where bulky groups increase and decrease the inhibitory activity, respectively. A large region of yellow contour at S 1 subsite suggests that bulky groups are not desired in this zone. This feature indicates that linear groups should be prioritized rather than branched groups in S 1 . Other yellow isopleths at S 3 subsite near the residues Ile174 and Tyr61 and green isopleths near the residues Glu217 and Leu99 indicate that bulky groups extended to the more hydrophilic part of S 3 decrease the inhibitory activity, but groups extended to the hydrophobic part of S 3 increase the inhibitory activity. The predictions of chemical features of substituents in S 3 according to the analysis of steric fields included in CoMSIA-SDA model are difficult since it is difficult to foresee their occupation and interactions in this wide cavity.
Cyan and purple isopleths in Fig 5B are in regions where HB donor groups favor and disfavor the activity, respectively. On the other hand, in Fig 5C magenta isopleths indicate regions where HB acceptors enhance the activity, and red isopleths indicate regions where HB acceptors decrease the activity. Interestingly, all the cyan and purple isopleths are located at S 1 subsite. The cyan and purple isopleths near Asp189, Gly218, Gly226, and Tyr228 represent 3D spaces that the model identified that are affected by the presence of HB donor groups in the deeper region of S 1 . The purple isopleth near Ser214, the red isopleth near Gly216, and the magenta isopleth near Cys191 and Glu192 indicate that the entrance to S 1 is affected by the presence of HB donor and acceptor groups. Two red isopleths in S 3 indicate that HB acceptor groups are not necessary in this zone. The first one is located near Asn98 and the other is exposed to solvent between Glu217 and Trp64. The magenta isopleth in S 2 near His57 indicates that HB acceptor groups such as CN are desired in this zone. Compounds 26,133,147,149,162, and 177 were selected as a sample for deriving free energy calculations using MM-GBSA method. The study includes a previous conformational sampling using MD simulations to consider averaged properties and to get more realistic conditions. Stability of the MD trajectories using the RMSD of the positions for all the protein atoms as a  Docking, QSAR, MD, and MM-GBSA of Thrombin Inhibitors function of simulation time was evaluated; RMSD was almost constant after the equilibration process (2.0 ns) for all the systems (Fig 6).

Free energy calculations results
The above mentioned compounds were selected from data set having in mind the differences in their activities against thrombin and their differences in the substitution and structural characteristics of the groups located at S 1 , S 2 and S 3 . In this sense, MM-GBSA paradigm allowed finding the role of physico-chemical terms in the binding potency when different groups located at S 1 , S 2 , and S 3 are changed. Guanidinooxy groups as substituent of N-ethylacetamide or n-propanol, and 2-amino-6-methylpyridin-5-yl groups were considered in S 1 . Methyl, Cl, F, CN, and oxo groups were considered in the aromatic ring in S 2 . 2,2-Difluoro-2-arylethyl, (1-arylcyclopropyl)methyl, arylethyl, cyclopentyl groups were considered as chemical features in S 3 .
MD allows analyzing the effect of the water media, the movement of the ligand inside the binding site, the movement of the residues in the binding site, and the stability of HBs. In general, from the analysis of all the simulations, it is possible to identify different water characteristics in the pockets S 1 , S 2 , and S 3 . There are many water molecules in S 3 which is the more exposed site to the water media, most of these molecules move freely and the other ones have occasional interactions with polar groups forming HBs. On the other hand, there are no water molecules in S 2 during MD simulations. Finally, there are some water molecules at S 1 trapped in HB networks formed between the ligands and the enzyme; these water molecules have a very limited movement. Interestingly, it is easy to differentiate S 1 , S 2 , and S 3 pockets according to the distinct water dynamical behavior, and this is highly related to the TI characteristics in each pocket disclosed by docking analysis. The dynamic of the studied ligands using MD simulations revealed that groups at S 1 form stable HB interactions, groups at S 2 are anchored in a hydrophobic cage, and groups at S 3 have a free movement.
The MM-GBSA free energy calculations were employed to get quantitative estimates for the binding free energies of the inhibitors inside the thrombin active site. We found a high correlation (R 2 = 0.751, Fig 7) between calculated ΔG values and the experimental ones (the experimental values were expressed as ΔΔG values with respect to the most active selected compound 133). However, only six compounds were considered; due to reduced data, we cannot conclude that MM-GBSA was able to explain the trend of the modeled TIs, but we can analyze the free energy components such as van der Waals (VDW), electrostatic, and solvation contributions to give detailed molecular information about the selected systems. The purpose is to evaluate the role of the physico-chemical features when different substituents are presented. The results of the predicted ΔG values and the physico-chemical components ΔE vdw , ΔE elect , and ΔG solv for the complexes were summarized in Table 3. To get a better view on which energy terms have more impact on the inhibitory potency, these individual energy components were compared. From Table 3, it can be seen that ΔE vdw has the major favorable contribution to the total free energy, but there is no big difference among this value for different complexes. In this sense, the VDW energy term is not primarily responsible for differentiating the binding affinity of the selected TIs. However, some of the most negative ΔE vdw values were obtained for the most Docking, QSAR, MD, and MM-GBSA of Thrombin Inhibitors active compounds 133 and 149. The solvation term also has the same effect in most of the complexes with a favorable contribution to the global free energy. Strangely, this term is positive for the most active compound 133. On the contrary, the electrostatic term was only favorable for compound 133, and the worst contribution of this term was for the less active compound 26.
The analysis of the MM-GBSA terms indicates that not any physico-chemical property has a preponderant role in the potency of TIs. It is expected that this behavior could be maintained if we consider the remaining compounds considering the complexity of the thrombin binding site in S 1 , S 2 , and S 3 pockets. We previously observed that these pockets have very different characteristics; therefore, each of them could be modulated considering different physicochemical features. However, the binding affinity can be the result of the combination of the particular effect of these features in each pocket.

Conclusions
In this work, we report a theoretical study on drug design area for thrombin inhibitors (TIs). We selected a data set of 177 compounds with diverse groups at motifs P 1 , P 2 , and P 3 , and explored the 3D positioning of them inside the thrombin active site pockets S 1 , S 2 , and S 3 by docking experiments. Our approach reproduced the previously reported position of compounds 13, 24, 128, 151, and 165 inside thrombin's active site. The remaining compounds were oriented in a similar manner. The interactions established for all compounds within thrombin binding site were carefully described.
Additionally, predictive QSAR models were built by using CoMSIA method. We found that CoMSIA-SDA including 138 compounds in the training set was the best among the explored models. This model included steric, HB donor, and HB acceptor fields. It was validated by LOO cross-validation and then it was used in the prediction of an external set which contains 34 compounds that were not included during the training process. The CoMSIA model could be used to predict novel candidates; moreover, the interpretation of the CoMSIA fields makes it possible to draw conclusions concerning the most appropriate features for novel analogues.
Finally, six compounds were selected and were subjected to MD and free energy calculations using MM-GBSA method. The dynamical behavior of the systems were revealed by MD; it was interesting to note that water molecules have completely different dynamical behavior in pockets S 1 , S 2 , and S 3 , which is related to the desired characteristics of motifs P 1 , P 2 , and P 3 of TIs. The free energy calculations using MM-GBSA allowed studying the role of physico-chemical parameters such as van der Waals interactions, electrostatic interactions, and solvation. It was found that not any component was preponderant in the study of the differential potency of inhibitors.