Exploring QSARs for inhibitory effect of a set of heterocyclic thrombin inhibitors by multilinear regression refined by artificial neural network and molecular docking simulations

Abstract

Several non-peptide heterocyclic compounds reported as potent thrombin inhibitors in vitro were chosen to carry out a QSAR study upon them using MLR and ANN analysis. In order to identify the best QSAR models, the input for ANN consisted of those subsets of descriptors used in the MLR models. The best QSAR models contained the SIC₀ descriptor as the main topological descriptor. To identify the physical and chemical properties involved in the ligand–thrombin complexes, an automated ligand-flexible docking procedure was used. The docking results suggest that the thrombin inhibition by these heterocyclic compounds is driven by π–π, hydrogen bonds and salt bridge interactions. The best Gibbs free energy of ligand binding was found at the thrombin sites S1 and D. We have shown that it is possible to build MLR models with geometries taken from two different sources (semi-empirical and MD geometries) and obtain similar results.

Keywords::

• Artificial neuron network
• docking
• quantitative structure–activity relationship
• heterocyclic thrombin inhibitors

Introduction

Thrombin is a multifunctional serine protease that plays a central role in coordinating the processes of blood coagulation¹, inflammation and wound repair following tissue injury^2–7. Recently, the search for a potent and orally administrated thrombin inhibitor for antithrombotic therapy has been one of the main challenges in medicinal chemistry since several problems were found in the early inhibitor types: for example, lack of selectivity, side effects or low oral bioavailability⁸. These inhibitors possessed reactive functional groups such as aldehydes, activated ketones or boronic acids⁸.

The non-peptide thrombin inhibitors with a heterocyclic core working as molecules with oral bioavailability gained special attention during the late 1990s decade^9–13. Recently, Ries et al.^14,15 have designed and developed orally active coagulation inhibitors with remarkably low molecular weight containing heterocyclic core structures as potent thrombin inhibitors in vitro. This compound design was based on the benzimidazole lead structure (Figure 1), by synthesis of alternative heterocyclic analogues and testing the activity. Through QSAR analysis, a wide variety of descriptors including constitutional, topological and quantum chemical have been reported^16–20, which have aided to develop knowledge-based drug-design methodologies.

Figure 1.  Benzimidazole lead structure.

Three recent papers proposed new methodologies that can be used to correlate the binding affinity data of several thrombin inhibitors with their molecular structure. Nicolotti et al.²¹ reported an approach to rationally interface structure- and ligand-based drug design through the automated generation of molecular docking alignments for 3D-QSAR using a genetic algorithm in a large series of 3-amidinophenylalanine inhibitors of thrombin. Fukunishi et al.²² proposed the hypothesis that “a model of active compound can be provided when integrating information of compounds high-ranked by docking simulation of a random compound library” where the inclusion of the true active compound in the high-ranked compounds is not necessary. Finally, Bonachéra and Horvath²³ reported the use of topological fuzzy pharmacophore triplets to establish and validate linear and non-linear QSAR models of a set of thrombin inhibitor models; in spite of the apparent success of the above mentioned methods, Nicolotti’s approach might be burdensome in terms of computational time, in the Fukunishi’s approach the QSAR model still provides higher performances than their pseudo structure–activity relationship (PSAR) models and in Bonachéra and Horvath’s approach their model revealed that some selected triplets make sense; however, their equations were unable to predict the activity of structurally different ligands.

Multiple linear regression (MLR) is still advantageous in building models with improved fit in the training set when the number of descriptors is small²⁴; further improvement can be achieved using artificial neural networks (ANN) a non-parametric non-linear modelling technique. The flexibility of ANN models for discovering more complex relationships has found a wide application in QSAR studies²⁵. Our working group has published some QSAR models^26–28 and also docking studies²⁹. Given that docking simulations yield important structural details about the ligand recognition on protein target²⁹, one can use the docking methodology³⁰ to corroborate the QSAR analysis results and to identify the recognition details of the ligand–target complex such as in pharmacophore groups³¹. Furthermore, computational docking is an efficient tool that has been applied in studies for structure-based drug design^32,33.

Generally speaking, the energy in docking studies is calculated as a function of five terms that include:

where the five ΔG terms on the right-hand side are coefficients empirically determined using linear regression analysis from a set of protein–ligand complexes with known binding constants³⁰. The summations are performed over all pairs of ligand atoms, i, and protein atoms, j. The in vacuo contributions include three interaction energy terms: a Lennard-Jones 12-6 dispersion–repulsion term, a directional 12-10 hydrogen-bonding term and a screened Coulombic electrostatic potential. The four terms correspond to the number of sp3 bonds in the ligand whereas the last term corresponds to the solvatation parameters, which are very important due to ligand recognition on proteins, which is carried out in the water environment³⁰.

In the present study, we have carried out a QSAR analysis over 50 thrombin inhibitor molecules (Table 1)^14,15,34 in order to derive a quantitative relationship between the thrombin inhibitory effects (IC₅₀, µM) and their chemical structural properties by MLR, ANN and docking studies. For these compounds, we have calculated a wide variety of descriptors such as quantum chemical, topological, Galvez topological charge index, 2D autocorrelation, RDF, 3D-MoRSE, WHIM and GETAWAY. In MLR analysis, the Leave-One-Out (LOO) and Leave-Group-Out (LGO) tests were used to select the most relevant QSAR models. A chance correlation study was performed on the optimal model via randomization test. The input for the ANN consisted of those subsets of descriptors used in different MLR models. Finally, an automated ligand-flexible docking (Autodock 4.0.1) was used to identify the chemical properties involved in the ligand–thrombin complexes.

Table 1.  Molecules considered in this study.

Code	Substituent	X	IC50 µM
MOL01		O	0.006
MOL02		O	0.018
MOL03		NH	0.037
MOL04		O	0.016
MOL05		O	0.6
MOL06		O	0.003
MOL07		O	0.006
MOL08		O	0.01
MOL09		O	0.092
MOL10		O	3.8
MOL11		O	1.3
MOL12		O	0.5
MOL13		O	0.066
MOL14		NH	0.067
MOL15		S	0.196
MOL16			1.5
MOL17			0.01
MOL18			0.027
MOL19			3.3
MOL20			0.016
MOL21			0.005
MOL22			0.008
MOL23			0.008
MOL24			0.012
MOL25			0.068
MOL26			0.013
MOL27			0.19
MOL28			0.015

Code	L	R1	R2	Ar	IC50 µM
MOL29	CH2	H	H	Phenyl	19
MOL30	CH2	CH2CH3	H	Phenyl	2.4
MOL31	CH2	CH2CH2CH3	H	Phenyl	4
MOL32	CH2	CH3	CH3	Phenyl	1.3
MOL33	CH2	CH3	H	3-Pyridyl	1.6
MOL34	CH2	CH3	H	2,5-Dimethoxyphenyl	1.8
MOL35	CH2	CH3	H	3,5-Bis-trifluoruromethyl-phenyl	>100
MOL36	CH2	CH3	H	1-Naphtyl	0.6
MOL37	CH2	CH3	H	2-Naphtyl	24
MOL38	CH2	CH3	H	5-Isoquinolyl	0.51
MOL39	CH2	CH3	H	8-Quinolyl	0.26
MOL40	CH2	CH3	CH2COOH	8-Quinolyl	0.12
MOL41	CH2	CH3	(CH2)3COOH	8-Quinolyl	0.13
MOL42	CH2-CH2	CH3	CH2COOH	8-Quinolyl	0.032
MOL43	CH2-O	CH3	CH2COOH	8-Quinolyl	0.058
MOL44	CH2-NH	CH3	CH2COOH	8-Quinolyl	0.011

Code	A	n	X	IC50 µM
MOL45	CH2	1	CH	0.054
MOL46	O	1	CH	0.33
MOL47	NH	1	CH	0.010
MOL48	NH	2	CH	0.0054
MOL49	NH	3	CH	0.010
MOL50	NH	2	N	0.0093

Methodology

Activity data and descriptor generation

There are at least 50 compounds that have shown thrombin inhibitory effect on factor XA (IC₅₀) in vitro, the reported data were taken from refs. (14,15,34). These compounds are characterized for containing several functional groups like quinoline, quinoxalinone, sulphonamides, carboxamide, aromatic rings, both hydrogen bond donating/accepting groups and an amine group.

All molecular structures were built as neutral species with the HYPERCHEM software³⁵. The molecular geometries were optimized under the PM3 semi-empirical approach³⁶ coupled to the Eigenvector following method³⁷ included in MOPAC 06 package³⁸. A conformational analysis was carried out to localize the global minimum of each molecule. Frequency calculations were performed to confirm the nature of the stationary point. Vibrational frequencies were unscaled.

In this study, a total of 1488 molecular descriptors were calculated, 1481 molecular descriptors belonging to 18 different types of theoretical descriptors were calculated for each molecule using the DRAGON software³⁹, and the rest seven quantum chemical descriptors were calculated using PM3 semi-empirical approach.

Data reduction, MLR and ANN modelling

The calculated descriptors were analysed for the existence of constant or near constant variables in about 85% of their elements⁴⁰, and those found were eliminated. After this process, 1195 descriptors were selected. In addition, for each one of the analysed data sets (vide infra) the correlation among the descriptors was investigated. Those descriptor pairs with high collinear relationship (r ≥ 0.90) were eliminated, leaving 366 descriptors to work with. This process was performed with a subroutine written in MATLAB (Version 5.3)^a.

The training models were built using the stepwise MLR subroutine of the SPSS for Windows 13.0 (2004)^b. The optimality of the linear models was determined by considering the t-statistical (models with all absolute values >4.0 were accepted) as well as the overall F-statistical test. The accepted linear models were subjected to a cross-validation (CV) analysis by LOO and LGO procedures to ensure that the models are not over-fitted or under-fitted. The LOO-CV and LGO-CV algorithms were written in MATLAB.

To demonstrate the absence of chance correlation on the best models obtained with the above procedure, we performed a Y-scrambling test^41,42 where the output value (the biological response) of the compounds was shuffled randomly and the scrambled data set was re-examined by MLR method against real (unscrambled) input descriptors to determine correlation and productivity of the resulting ‘model’. This procedure was repeated 30 times on many different scrambled data sets. R²_CC and LOO cross-validation (R²_CV-CC) of the resulting models were averaged and compared with the real model⁴³. The corresponding algorithm was written in MATLAB.

The accepted MLR models were further optimized using a feed-forward neural network (NN) with back propagation of error algorithm. This ANN subroutine was written in MATLAB (version 7.2.0)^c. After several tests, the best results were obtained with a (N × 2N × 1) NN with one hidden and one output layer, where N is the number of descriptors obtained from the best MLR models.

Molecular modelling procedure

The 50 known thrombin inhibitors previously mentioned were also studied by docking simulation employing AutoDock 4.0.1³⁰. The 3D structure of the ligands in their minimum-energy conformation was obtained by PM3 semi-empirical approach³⁶ coupled to the Eigenvector following method³⁷ included in MOPAC 06 package³⁸.

To identify the ligand–thrombin recognition-binding sites, all the possible rotatable bonds, the torsional degrees of freedom, the atomic partial charges (Gasteiger) and merge-non-polar hydrogens of the ligands were assigned using AutoDockTools 1.5.4, which is included in the AutoDock program³⁰.

The 3D structural model for human thrombin was obtained from the Protein Data Bank (PDB code: 1A4W). Afterwards, hydrogens of human thrombin at pH ~7 were added. The Kollman charges for all atoms of the enzyme were assigned with AutoDockTools 1.5.4. The ligands were docked under the blind docking procedure (grid box of 124 × 124 × 124 Å with the grid points separated by 0.375 Å) on the enzyme by means of AutoDock software under the hybrid Lamarckian Genetic Algorithm, with an initial population of 100 randomly placed individuals and a maximum number of energy evaluations of 1 × 10⁷. The resulting docked orientations within a root-mean square deviation of 0.5 Å were clustered together. The lowest energy cluster returned by AutoDock for each compound was used for further analysis. All other parameters were maintained at their default settings. All ligand–enzyme complexes were visualized using the Visual Molecular Dynamics v.1.8.6 (VMD) program⁴⁴.

Data analysis

The selected 50 thrombin inhibitors were arranged in four different data sets, first to test their independence on the biological data source and second to test the prediction strength of MLR analysis versus ANN analysis upon a general quinoline, quinoxalinone derivative data set.

Having in mind that the Ockham’s razor⁴⁵ suggests that the larger the ratio of number of molecules to the number of properties, the better expected predictive power of the model. In the following procedure, we applied the commonly accepted lower limit of 5 for this value⁴⁶.

The first data set considered all molecules tested (50 molecules) as the training sample. This data set was called M50. To construct the next three data sets, five molecules were eliminated each time at random from the M50 set and allowed them to form part of the external prediction set. The remaining molecules were used to construct three different sets by taking at random 36, 32 and 27 compounds, respectively, as the training sample; the remaining 9, 13 and 18 molecules in each set were part of the prediction sample. These sets were called M36, M32 and M27, respectively.

Despite that many authors consider a high q² (for instance, q² > 0.5) as an indicator or even as the ultimate proof that the model is highly predictive⁴⁷, to define the relative importance of the QSAR models, we have also performed several other tests beyond the q². These include the Akaike’s information criterion (AIC)⁴⁸, the Kubinyi fitness function (FIT)⁴⁹ and the calculation of the correlation coefficient for regression through the origin (RTO, R²₀) were performed for the training and predicted activity sets⁵⁰.

Finally, we have also taken into consideration the suggestion made by Guha and Jurs⁴⁰ who have shown that compounds that are very similar in the training set, generally exhibit smaller residuals and standard errors of prediction. These authors suggested that the cutoff value must be selected so that the outliers are approximately 20–30% of the whole data set.

Results

Pool of descriptors

The first step was to select the best descriptors for data set M50. After application of the stepwise MLR subroutine from the SPSS software, 21 descriptors were chosen from 10 different groups (see Table SI1 in the supporting information). The best model obtained was built with 17 variables from the reduced descriptor pool; however, the relation molecule number:descriptor is 2.94, which is lower than the ideal value of 5⁴⁶, in this sense the maximum number of variables should be 10 (see Equation 1). In the resulting equation, the main descriptors correspond to one topological descriptor and several chemical descriptors weighted by atomic van der Waals volumes.

or

1

R = 0.9526, SD = 0.3275, P < 0.0001, N = 50, SEE = 0.3816, q²_LOO = 0.6502, q²_LGO = 0.6650, FIT = 2.5725, AIC = 0.1420, F = 38.218.

The main descriptor in this equation is SIC₀, a topological descriptor developed by Basak et al.⁵¹; this descriptor takes into account all atoms in the constitutional formula (included hydrogen atoms) and considers the information content provided by various classes of atoms based on their topological neighbourhood. The SIC_r descriptors (r is the order level and could be 0 to 6) are not intercorrelated with other topological indices. Particularly, at the zero-order level the atom set is partitioned solely on the basis of its chemical nature. The SIC₀ descriptor has a regression coefficient of R = 0.693 with respect to the experimental biological activity of the 50 compounds.

The correlation matrix of this model showed no colinearity among the chosen descriptors. The largest colinearity value is found between SIC₀ and ATS6v descriptors (R = −0.715). The ATS6v descriptor is a 2D autocorrelation descriptor; the computation of this descriptor involves the summations of different autocorrelation functions corresponding to different fragment lengths and leads to different autocorrelation vectors corresponding to the lengths of the structural fragments⁵². As a result, this descriptor addresses the topology of the structure or parts thereof in association with a selected physicochemical property.

For a cutoff = 1 used in all experiments performed in this study, only MOL19 was present as an outlier (see Figure SI1 in the supporting information). No apparent geometrical characteristic is responsible for its behaviour. The structure of MOL19 is very similar to that of MOL18, which includes an ester group that is more voluminous than the acid group of MOL19. It is important to mention that the IC₅₀ of MOL18 and MOL19 are 0.027 and 3.3 µM, respectively; this means that in terms of inhibition of thrombin, the ester group is more important than the acid group.

Finding the best thrombin inhibitor model

To determine the optimum number of descriptors arising from data set M50, several training:prediction ratios were studied (see ‘Data analysis’ section). The compounds included in these ratios were taken randomly. For each data set, three different subsets were built; in total, nine subsets were studied by MLR methodology. The best inhibitor of thrombin MLR models were selected according to the best statistical parameters. Finally, the ANN analysis was applied to the best MLR models.

MLR models

The maximum number of descriptors for data set M36, M32 and M27 will be 7, 6 and 5, respectively; on the other hand, the best models of each data set were built using 4 (M36-MLR), 3 (M32-MLR) and 3 (M27-MLR) descriptors (see Table SI2 in the supporting information). These three models present the same structure on the descriptor type, that is, the SIC₀ descriptor is the main variable as in M50 MLR model (see ‘Pool of descriptors’ section). The RBF constitutional descriptor¹⁹ is taken to form part of these three models. This descriptor corresponds to the rotatable bond fraction in the molecule; the negative term associated with the descriptor indicates that the fractional decrease in the rotatable bond in the molecule is conductive for the thrombin inhibitory activities of the derivatives; for example, molecule MOL06 is the most potent compound whose X-ray structure with thrombin revealed an unexpected binding mode of this lead novel coagulation inhibitor structure¹⁴, where the ester group interacts with the S4 subsite. Two of the cyclopentyl methylene groups point into the bulk solvent indicating that they do not contribute to the enzyme affinity. The rigidity of the cyclopentyl group is an important factor of this enzyme affinity.

Statistically, the M27-MLR model is the finest model that has a correlation coefficient of 0.8725, q²_LOO value of 0.6737, FIT parameter of 2.2784 and AIC parameter of 0.2528; therefore, the M36- and M32-MLR models show statistical values close to M27-MLR values (see Table SI2 in the supporting information).

In the three models, the cutoff was defined at ±1.0 logarithmic units from the ideal line; this value is less than the common cutoff value that is three times the standard deviation (in these cases, this value should be between 1.395 and 1.512 units). A small number of outliers were found in each MLR model. The M36-MLR model shows two outliers from the Training set (MOL40 and MOL04), whereas the M32-MLR model also shows two outliers, one from the training set (MOL37) and one from the prediction set (MOL31), and finally the M27-MLR model presents four outliers: one from the training set (MOL37) and three from the predicted set (MOL13, MOL14 and MOL31). In these three models, validation sets molecules are always inside of the cutoff range.

To obtain the best statistic representation and to decide on the best model, the ANN techniques were employed.

The Y-scrambling test showed that the average for R²_CC and R²_CV-CC are 0.1220, 0.1287, 0.1165 and 0.3375, 0.1836, 0.1198 for M36-MLR, M32-MLR and M27-MLR, respectively. The lower values in the predicted set in M36-MLR decrease the probability of chance correlation effect. From these results, we concluded that increasing the number of molecules in the training set or decreasing the number of molecules in the prediction set increases the probability of chance correlation.

ANN models

To seek improvement on the MLR results, the ANN methodology was applied to M36-MLR, M32-MLR and M27-MLR models. We used a feed-forward NN with a back propagation of error algorithm. After several trial and error tests, the best models were obtained for a N × 2N × 1 network, with one hidden and one output layer, where N is the number of descriptors of the best MLR models.

Table 2 shows the equations and statistical parameters of the three models obtained using the ANN methodology. The statistical parameters for the three models were improved when using the ANN methodology where the correlation coefficients were increased to values around 0.90 units. Moreover, the q² values were increased.

Table 2.  The three best models obtained by the ANN methodology.

Model	Equation
M36-ANN	log IC50^cal=−0.1918 (±0.0946) + 0.8258 (±0.0651) log IC50^Exp R = 0.9085, SD = 0.4105, P < 0.0001, N = 36, q^2LOO = 0.8578
M32-ANN	log IC50^cal=−0.0852 (±0.1088) + 0.8892 (±0.0772) log IC50^Exp R = 0.9031, SD = 0.4726, P < 0.0001, N = 32, q^2LOO = 0.9100
M27-ANN	log IC50^cal=−0.2599 (±0.1190) + 0.7983 (±0.0757) log IC50^Exp R = 0.9036, SD = 0.4475, P < 0.0001, N = 27, q^2LOO = 0.7892

Each equation is represented as function of log IC₅₀^exp.

In terms of outliers, the behaviour is very similar. The M36-ANN and M32-ANN models continue with two outliers, whereas the M27-ANN model shows a decrease in this number. The outliers in ANN models are almost the same found in the MLR models; in M36-ANN model, only the MOL39 molecule is a new outlier in place of MOL04 molecule.

To establish which one of these models is the best model, RTO calculations were performed upon the training and predicted sets (Table 3), the best model will be the one that shows the correlation coefficients close to 1 (R₀ and R₀′ for training and predicted set) and slopes close to 1 too. Three models show correlation coefficients and slopes close to 1 and therefore it is not possible to determine with these parameters which model is the best. It is clear that the molecule randomization does not have influence in the different models; at this point, the selection of the best model should be by the number of descriptors and the number of molecules in the model. After performing the MLR, ANN, RTO and descriptor number analysis, we concluded that the best model is M27-ANN (see Figure SI2 in the supporting information).

Table 3.  The three best models obtained by the RTO methodology.

Model	Equation
M36-RTO	log IC50^cal= 0.9910 log IC50^Exp
	R0 = 0.9532, P-param = 0.0000, N = 36, F = 348.2199
	log IC50^cal = 1.1221 log IC50^Exp
	R0′ = 0.9536, P-param = 0.0000, N = 9, F = 80.2186
M32-RTO	log IC50^cal = 0.9579 log IC50^cal
	R0 = 0.9428, P-param = 0.0000, N = 32, F = 248.0019
	log IC50^cal = 1.2486 log IC50^cal
	R0′ = 0.9765, P-param = 0.0000, N = 13, F = 246.3107
M27-RTO	log IC50^cal = 0.9899 log IC50^cal
	R0 = 0.9503, P-param = 0.0000, N = 27, F = 242.3091
	log IC50^cal = 1.0608 log IC50^cal
	R0′ = 0.9311, P-param = 0.0000, N = 18, F = 110.6911

Each equation is represented as function of log IC₅₀^exp, first equation corresponds to training set and second equation corresponds to predicted set.

Docking results

This protein is characterized by having sites named D-pocket, P-pocket and the most important binding site the Ser195 amino acid, this is a site S or S1. In this work, we identified that the vast majority of the ligands tested reached the binding site S1 by making several hydrogen bonds, electrostatic interactions and π−π interaction as has been reported by Frédérick et al.⁵³ Most ligands can reach the three sites at the same time being more common site S1 and the D pocket (see Table SI3 in the supporting information). This is because these two sites are very close to each other and the ligands have a volume that cover the entire cavity and they also share the same pharmacophore group. In addition, some ligands do not reach the reported sites as lowest energy complexes (see Table SI3 in the supporting information).

The binding energy in the enzyme–ligand complex is calculated from the addition of the next five energies (see Table SI4 in the supporting information): the ΔG_Intermol (final intermolecular energy); ΔG_vdW (van der Waals energy + hydrogen bond energy + desolvation energy); ΔG_Elect (electrostatic energy); ΔG_internal (final total internal energy); ΔG_Torsional (torsional free energy)³⁰. Due to the lack of several biological properties, it is very common not to find relationships between the theoretical and experimental data. However, after studying each energy contribution, a relationship between four descriptors taken from the M50-MLR and ΔG_internal were found as suitable, this model is called M50-MLR-docked (Figure 2, Equation 2). It is important to mention that the main descriptor in the M50-MLR model is the same descriptor as in the M50-MLR-docked model, that is, the QSAR models from semi-empirical and docking calculations are similar. The equation obtained for the M50-MLR-Docked model is:

Figure 2.  Training set of the M50-MLR-docked data set. Calculated equation by multiple linear regression (solid line), ideal QSAR line (dotted line) and residual line limits (dashed line, cutoff = |1.0|). Outliers are also depicted.

or

2

R = 0.883, SD = 0.469, P < 0.0001, N = 50, SEE = 0.5544, q²_LOO = 0.7175, q²_LGO = 0.6580, FIT = 2.0762, AIC = 0.3006, F = 31.100.

In this model, three molecules are considered as outliers (Figure 2), again MOL19 is an outlier as in M50-MLR model.

Discussion

Brandstetter et al.⁵⁴ and Banner and Hadváry⁵⁵ studied the interaction between the bovine–thrombin complex with the benzamidine-based inhibitor NAPAP, an indicator of what molecule type could present some biological activity. This analysis was performed by comparison of several molecular interactions, for instance, the amidine group forms a bidentate salt bridge with the carboxylate of Asp189, hydrophobic interactions are formed by the piperidine ring and the naphthyl moiety with the proximal and the distal pocket of the thrombin active site. Additionally, the bridging glycine moiety of NAPAP, in analogy to peptidic serine protease substrates, forms the canonical hydrogen-bonding patterns with residue Trp215 and Gly216 at the rim of the S1 pocket. On the other hand, MOL16 and human α-thrombin form a complex³⁴ with the same salt bridge between Asp189 and MOL16 as found in the NAPAP complex, where the terminal phenyl ring is bound by a hydrophobic interaction to the D-pocket and the N-methyl group fits nicely into the P-pocket. The hydrogen bond between the inhibitor and the enzyme is not present; this last interaction should be very important in the inhibitor activity because the NAPAP and MOL16 show IC₅₀ = 0.2 µM and IC₅₀ = 1.5 µM, respectively.

In the case of NAPAP–MOL06, the first hydrophobic interaction is formed by the quinoline ring and its methyl substitution and the second hydrophobic interaction is formed by ethyl moiety or cyclopentyl moiety with the proximal and the distal pocket of the thrombin active site. These types of interactions could be very important in setting up the biological activity. On the other hand, NAPAP–MOL35 shows that p-CF₃ (and o-CF₃) substitution in the phenyl moiety is important to decrease the biological activity, that is, the higher polar substitution in the hydrophobic interaction, the less biological activity (see Figure SI3 in supporting information).

In M27-MLR model, the leader compound (MOL06) presents the smaller value of the SIC₀ descriptor, whereas the MOL35 compound presents the worse biological activity that shows a higher value of the SIC₀ descriptor, that is, the SIC₀ descriptor is the main descriptor in the biological activity correlation. The linear regression between the biological activity and the SIC₀ descriptor, gives a correlation coefficient of 0.748.

The leader compound is located in the prediction set and the worse biological activity compound is located in the training set. In the linear regression, residual values are |0.6172| and |0.4648|, respectively, whereas in the M27-MLR model, these values change to |0.7559| and |0.2662|, in both cases, the residuals are inside the outlier region. The best residual of leader compound is found in the M27-ANN model with |0.4406|.

Docking

Several types of interaction drive the ligand–protein recognition in the thrombin system. We discussed what happens when MOL01 is the ligand (see Figure 3A) for which the binding energy is ΔG = −7.6 kcal/mol. In this figure, several ligand–protein interactions shown explain the best affinity; for example, there are several π−π interactions with Tyr228, Trp60D, Tyr60A, and the R of the compounds (that is an aromatic ring) also makes π−π interactions with Trp215. This confirms that the π−π interactions are very important for stabilizing these molecular complexes as has been reported by others⁵⁶. In addition, a cluster of hydrogen bonds is made with the backbone of Gly216, Gly219 and Phe227. Finally, as have been reported elsewhere, electrostatic interaction of salt bridge is carried out between the amine of the compound with carboxylic acid of the Asp189. This means that aromatic substituents with few hindrance effects could be necessary to increase the affinity, chemical characteristics that MOL01 has. The same is observed for ligands from MOL01 to MOL15.

Figure 3.  Ligand interactions of the compounds with best predicted affinity on thrombin (pdb code: 1A4W) classified according to their chemical cores, (A) sample belonging to MOL01–MOL15 subset, (B) sample belonging to MOL17–MOL28 subset, (C) sample belonging to MOL29–MOL44 subset and (D) sample belonging to MOL45–MOL50.

Other set of compounds that share the same core was the set of compounds from MOL17 to MOL28, being MOL23 the one with the best affinity (ΔG_binding = −7.89 kcal/mol). As it occurs with MOL01, there are several π−π interactions with aromatic residues such as Tyr60A and Trp60D and hydrogen bonds with backbone of Gly216, Gly219 and Phe227. Also a π–cation interaction with the amine group and the aromatic ring of the Tyr228 could take place. Additionally, hydrophobic interactions are found between the cyclohexane group of the compound and the aromatic ring of Trp215, since the carbonyl group is close to the cyclohexane it can make a hydrogen bond with the NH group of the aromatic ring of Trp215 and with the amine group of the Gly216 backbone. These ligand–protein interaction descriptions are shown in Figure 3B. With these results, it is possible to infer that the chemical core is necessary to fit in the thrombin protein, that the R group that normally is bound to Trp215 is required to be hydrophobic and not necessarily an aromatic ring; however, it is important that there must not be high hindrance effects.

Figure 3C depicts the best compound belonging to the set from MOL29 to MOL44; this is MOL41 with ΔG_binding = −8.14 kcal/mol, which shares the same chemical core while changing a six member cycle by a five member cycle (Table 1). In this compound, a binding mode occurred for which we did not find a previous report. As one can see, the aromatic system could make several interactions with Trp215, Trp60D and Tyr60A. Some hydrogen bonds are also observed between the carbonyl group from S and Gly216 backbone, and the amine group of five cycle member with the hydroxyl group of Tyr60A.

Finally, Figure 3D shows the ligand that had the best affinity on the thrombin-binding site; this is MOL46 (ΔG_binding = -8.18 kcal/mol), which belongs to the group of compounds from MOL45 to MOL50. These compounds have the same chemical core present in the above described compounds (MOL29 to MOL44); however, the set of ligands from MOL45 to MOL50 had more hindrance effects due to the aromatic ring. MOL46 makes several π−π interactions with Trp60D, Tyr60A and Trp215 and several hydrogen bonds with the hydroxyl group of Ser 195, Lys60F and the hydroxyl group of Tyr60A.

Conclusion

In this study, we have explored a set of 50 heterocyclic molecules under QSAR and docking studies. Several models were studied with MLR and ANN methodology where it is found that the chemical nature of some atoms (SIC₀ descriptor) in the molecules was the most important descriptor. QSAR shows that the main interaction between the heterocyclic and the target protein is by hydrogen bond interactions, whereas docking studies show electrostatic and π−π interactions are important. These results suggest that it is possible that other hydrophobic substituents in the phenyl moiety could be very important in decreasing the biological activity.

The observed difference in the QSAR models obtained with MLR might be the effect of the strategy used in assigning the type and given number of molecules present in the training/predicted sets. Such effect has been reported by the Hemmateenejad group⁵⁷.

Finally, we have shown that it is possible to build MLR models with geometries taken from two different sources, such as semi-empirical or MD methodologies, and obtain similar results. An option to improve the prediction is to separate the different contributions of the binding energy and to take them all as independent descriptors.

Acknowledgements

RGJ acknowledges partial support from CONACyT grant 46061-R. JCB acknowledges partial support from CONACyT grant 62488, ICyTDF and Fundación Miguel Aleman AC, EDI/SIP-SIBE/COFAA-IPN.

Declaration of interest

The authors report no conflicts of interest.

Notes

a MATLAB Version 5.3 1.29215a (R11.1), The MathWorks Inc., Natick, MA, 1999.

b SPSS Version 13 for Windows, Chicago, IL, 2004.

c MATLAB Version 7.2.0.232 (R2006a). The MathWorks Inc., Natick, MA, 2006.