Investigation on the binding mode of benzothiophene analogues as potent factor IXa (FIXa) inhibitors in thrombosis by CoMFA, docking and molecular dynamic studies

Abstract

Recently, benzothiophenes attract much attention of interest due to its possible inhibitory activity targeting FIXa, a blood coagulation factor that is essential for the amplification or consolidation phase of blood coagulation. To explore this inhibitory mechanism, three-dimensional quantitative structure–activity relationship (3D-QSAR), molecular docking and molecular dynamics (MD) studies on a series of 84 benzothiophene analogues, for the first time, were performed. As a result, a highly predictive CoMFA model was developed with the q² = 0.52, r² = 0.97 and r²_pred = 0.81, respectively. The CoMFA contour maps, the docking analysis, as well as the MD simulation results are all in a good agreement, proving the reliability and robustness of the model. These models and the information, we hoped, would be helpful in screening and development of novel drugs against thrombosis prior to synthesis.

Keywords::

• FIXa
• benzothiophene
• docking
• CoMFA
• MD

Introduction

Thrombotic diseases, such as myocardial infarction, stroke, deep vein thrombosis and pulmonary embolism¹, are leading causes of mortality. Currently, these diseases are most often treated with heparin or warfarin, both of which carry significant risks of bleeding complications¹. Therefore, there is a great medical need for superior anti-thrombotic agents with a wider therapeutic index. In the past few years, researchers have found that FXa plays a central role in the regulation of normal homeostasis and abnormal intravascular thrombus development. Thus, inhibition of FXa furnishes effective approaches for controlling blood coagulation process and for the treatment of thrombotic diseases.

As part of the effort for developing fast and low-cost tools for facilitating the discovery of new anticoagulants, several computational methods have been explored for predicting the FXa inhibitors. Molecular docking method has been used to find potential inhibitors by identifying those compounds that can be docked into the inhibitor binding site of the 3D structure of FXa, which achieves accuracies of 80% for inhibitors and 85% for non-inhibitors based on the test of 112 ligands². Specific structural and physicochemical properties of the known inhibitors have been used to derive 3D quantitative structure–activity relationships with a reported inhibitor prediction accuracy of 84–88% based on the test of 279 inhibitors³. Although many selective and potent FXa inhibitors have been advanced to the clinical development, the potential risk of bleeding and a narrow therapeutic window remain for these newer anticoagulants, and their safety profiles are yet to be established in clinical trials⁴.

FIXa is an important blood coagulation factor that is essential for the amplification or consolidation phase of blood coagulation that results in thrombus formation via FX activation to FXa⁵. A recent study shows that patients who are heterozygous gene carriers of FIXa have decreased coagulability and are protected from ischemic heart diseases⁶. As FIXa is activated both by the stimulation of the intrinsic system and by low level of TF, it is believed that selective inhibition of FIXa would be effective in anti-thrombosis in low tissue factor sites but not in high tissue factor environments, which will allow for intravascular anticoagulation with maintenance of extravascular hemostasis⁷. Consequently, at the initiation stage of the cascade, the upstream inhibition of FIXa would represent a wonderful method for the development of anticoagulants with good selectivity and safety. Recently, several classes of compounds, such as 5-amidinoindoles¹, pyrazole analogues^8,9, 5-amidinobenzo[b]thiophenes¹⁰, benzothiophenes^11,12 and so on, have been reported as FIXa inhibitors. Among them, the benzothiophenes, designed on the basis of FIXa X-ray crystallography, illustrate high potential and selectivity for FIXa. However, to our best knowledge, there is still no report of in silico modelling on these benzothiophene-based FIXa inhibitors up to date. Therefore, it should be beneficial to explore the quantitative structure–activity relationship of structurally diverse benzothiophenes as potent FIXa inhibitors by computational approaches.

Three-dimensional quantitative structure–activity relationship (3D-QSAR) methods, especially, the popular comparative field analysis (CoMFA)¹³ provides a direct way to explore and visualize the structure–activity relationship of the molecules being studied. But in all these studies, the molecule alignment and active conformation determination, as is well-known, are very important that they greatly affect and sometimes even determine the success of a model. In most cases, the bound protein–ligand complex is not available; therefore, a computation method has to be deployed to determine the conformations and alignment of the set of ligand molecules based on the ligand-based methods, such as the atom-based alignment¹⁴, database and field fit alignments^15,16, pharmacophore-based superimposition¹⁷ and so on. However, once the crystallography presents, docking is an attractive way to align the molecules. It has been shown that the order of preference for the alignment selection for 3D-QSAR model development may be DCBA (docked conformer-based alignment) over GMCBA (global minima energy conformer-based alignment)^18,19. Many successful applications of docking alignment with CoMFA have been reported^20–22. In the present work, since the structural information of FIXa in the PDB database (the protein data bank) is available (PDB ID: 3LC5), the determination of the ‘active’ conformation of each molecule and the molecular alignment are carried out in the present work using the flexible docking program, Surflex-dock²³. Additionally, the stability of conformation for docked complex structure was tested by molecular dynamics (MD) simulation. Thus, a combination of receptor docking, 3D-QSAR and MD simulation approaches would be desired for the investigation of the ligand libraries as FIXa inhibitors.

Validation is another crucial aspect of QSAR modelling. In this work, to establish a statistically robust model capable of making accurate and reliable predictions of biological activity of compounds, the developed model is also validated carefully according to Tropsha and so on^24–27.

Thus, the aim of the present study is to quantitatively investigate the structure–activity relationships of benzothiophene-based FIXa inhibitors with regard to the future development of potential new drug leads against coagulation by using the obtained statistically significant CoMFA analysis based on the docked conformer-based alignment. Furthermore, our work presented the first 3D-QSAR study for a set of 84 benzothiophene compounds, which might provide a platform for the screening of novel inhibitors and enables the interpretation of their binding modes to FIXa.

Materials and methods

Dataset

A diverse dataset of 84 benzothiophene analogues with the experimental values for the K_i of the FIXa inhibitors were taken from literatures^11,12 that are published by the same research group. The activity of all the molecules was measured by the same assay on a recombinant form of human factor IX, which was prepared to overcome potential problems with the use of ethylene glycol in a factor IXa inhibition assay¹¹. Here, the converted molar pK_i (−logK_i) values, ranging from 4.57 to 8.70 moles, were used as the dependent variables in the QSAR regression analysis to improve the normal distribution of the experimental data points. A span of 4 log units of the pK_i values and a large dataset made it highly appropriate for a QSAR analysis. As for the division rule for training and test sets, previous studies have provided some valuable strategies such as the random selection²⁸, activity–range algorithms, K-means cluster-based selection, self-organizing maps^29,30 and sphere–exclusion algorithm^27,31–33. Of the investigations they have made, the separation rule was all emphasized to the key point that how to make the training set to represent the entire dataset to the most. Hence here, 16 compounds (a ratio of about 4:1 between the training and the test sets) were selected as the test set by considering the criterion that the test set should represent the structural diversity and a range of FIXa inhibitory activity similar to that of the training set. In this work, the mean of the biological activity of the training and test sets is 6.88 and 6.81, respectively. Table 1 depicts the structures and activity of all the compounds in the dataset.

Table 1.  The structures, experimental and predicted pK_i activity of benzothiophene analogues.

No.	Substituent		Experimental	Predicted	Residual	Ref.^a
Training set
1	4-I		6.02	5.76	−0.26	11
3	4-OBn		5.82	5.91	0.09	11
4			6.1	6.28	0.18	11
5	4-Ph		5.8	5.75	−0.05	11
6	5-Ph		5.1	5.30	0.20	11
7			5.23	5.30	0.07	11
10	6-OCH2COOEt		5.42	5.24	−0.18	11
11			4.62	4.54	−0.08	11
14	4′-F		6.15	6.04	−0.11	11
15	2′-Cl		5.14	5.21	0.07	11
16	3′-Cl		6.05	6.04	−0.01	11
17	4′-Cl		6.35	6.09	−0.26	11
18	2′-MeO		4.57	4.68	0.11	11
19	3′-MeO		6.35	6.47	0.12	11
21	4′-Me		6.05	6.22	0.17	11
22	4′-CO2Me		4.95	5.09	0.14	11
25	3′,4′-Cl,F		7.51	7.14	−0.37	11
26	PhCH2		4.67	4.50	−0.18	11
27	3-Py		5.98	6.32	0.34	11
28	2-Thiophen		6.46	6.43	−0.03	11
	Ar	R
29	4′-F-Ph	Me	4.95	5.59	0.64	11
30	4′-F-Ph	CONH2	5.7	5.58	−0.12	11
31	4′-F-Ph	CO2Me	7.4	7.13	−0.27	11
32	4′-F-Ph	CO2H	7.68	7.78	0.10	11
33	4′-Cl-Ph	CO2Me	7.7	7.70	0.00	11
34	4′-Cl-Ph	CH2CO2Me	7.55	7.30	−0.25	11
35	4′-Cl-Ph	(CH2)2CO2Me	7.44	7.59	0.15	11
36		CO2Et	8.52	8.37	−0.15	11
37	—		6.57	6.78	0.21	12
38	—		5.29	5.49	0.20	12
39	—		6.30	6.21	−0.09	12
40			5.08	4.98	−0.11	12
41			6	5.98	−0.02	12
42	Ph		6.77	6.83	0.06	12
44	2′-MeO-Ph		6.85	6.87	0.02	12
45	3′-MeO-Ph		6.89	6.86	−0.03	12
46	4′-MeO-Ph		6.85	6.60	−0.25	12
47			6.38	6.37	−0.01	12
48	4′-(Me2NCO)-Ph		6.43	6.64	0.21	12
49			6.54	6.45	−0.09	12
50	NH2		6.8	6.63	−0.17	12
51	NHCOMe		6.47	6.33	−0.14	12
52	NHBn		7	7.04	0.04	12
54	NHCONHPh		7.54	7.23	−0.31	12
55	OH		6.57	6.69	0.12	12
56	OCONHMe		7.01	6.99	−0.02	12
59	OCONHBn		7.21	7.25	0.04	12
60	2′-F		8.3	8.18	−0.12	12
61	3′-F		7.6	8.08	0.48	12
62	4′-F		8.1	8.34	0.24	12
63	2′-Me		8.1	8.03	−0.07	12
64	4′-MeO		8.3	8.41	0.11	12
66			8	8.18	0.18	12
67	2′-CH2NH2		8.4	8.67	0.27	12
69	2′-CH2NHBu-i		8.7	8.44	−0.26	12
70	2′-CH2NMe2		8.1	8.13	0.03	12
71			8.52	8.41	−0.11	12
73	2′,6′-F,F		8.15	8.23	0.08	12
74			7.92	7.61	−0.31	12
75			8	7.67	−0.33	12
76			8.05	7.96	−0.09	12
77			8.3	8.23	−0.07	12
78			8.1	8.09	−0.01	12
79			7.96	7.92	−0.04	12
80			8.3	8.48	0.17	12
81			8.52	8.60	0.08	12
83	2′-CH2NH2		8.36	8.45	0.09	12
84	2′-CH2NMe2		8.09	8.07	−0.02	12
Test set
2	4-OCH2COOEt		4.92	5.27	0.35	11
8	5-OBn		5.72	5.89	0.17	11
9	6-OBn		5.42	5.26	−0.16	11
12	2′-F		4.64	5.17	0.53	11
13	3′-F		6.17	5.61	−0.56	11
20	4′-MeO		6.6	6.06	−0.54	11
23	3′,4′-F,F		6.35	5.73	−0.62	11
24	3′,4′-F,Cl		7.12	6.18	−0.94	11
43	Bn		6.41	6.27	−0.14	12
53	NHCO2Ph		7.3	7.09	−0.22	12
57	OCONHEt		7.55	6.98	−0.57	12
58	OCONHPh		8	7.30	−0.70	12
65	4′-O(CH2)2NMe2		8.4	7.56	−0.84	12
68	2′-CH2NHMe		8.52	7.98	−0.54	12
72	2′,4-F,F		8.15	7.82	−0.33	12
82	H		7.81	7.91	0.10	12

Reproduced with permission from refs. [11] and [12].

^aFrom the corresponding reference.

In the modelling process, taking the crystal structure of IZX (i.e. compound 37 in this work) as the starting conformation, the rest of the molecules were created using the sketch molecule module of SYBYL 6.9 package (Tripos Associates, St. Louis, MO). A constrained minimization followed by full minimization was carried out on these molecules to prevent the conformations from moving to false regions. Partial atomic charges were calculated by the Gasteiger–Hückel method³⁴, and energy minimizations were performed by using the Tripos force field³⁵ and the Powell conjugate gradient algorithm with a convergence criterion of 0.05 kcal/mol Å.

Docking simulation and structural alignment

Docking simulations of benzothiophene analogues into the FIXa-binding pocket were performed using the Surflex-dock module of SYBYL package in this study. Surflex-dock utilizes a so-called “whole” molecule alignment algorithm based on morphological similarity between the ligand and target³⁶. The docking method aligns the ligand to a “protomol” or idealized ligand in the active site of the target. The detailed algorithm for Surflex-dock is presented in the literature²³. For our studies, the X-ray crystal structure of FIXa with high resolution (2.62 Å) was retrieved from RCSB Protein Data Bank (PDB ID: 3LC5). Prior to docking, all ligands were extracted from the crystal structure with water molecules considered, and hydrogen atoms were added to the protein in standard geometry using the biopolymer modulators. In this study, Ligand-based Mode was adopted to generate the protomol in Surflex-dock program, and two parameters that significantly affect the size and extend of the protomol generated are the threshold and the bloat values. In the present work, the threshold and bloat values were set to 0.5 and 0, respectively. During the present molecular docking process, the protein was considered to be rigid, and the ligand molecules were flexible. Other parameters were established by default values in the soft. In the current work, 20 conformations were obtained though Surflex-dock for each ligand and all conformations were extracted from the optimized inhibitor–FIXa complex. The conformations with the highest total scores for each ligand of the training set were aligned automatically together inside the binding pocket of FIXa and used directly for CoMFA to explore the specific contributions of the electrostatic and steric effects on the molecular bioactivities.

CoMFA analysis

A training set of 68 compounds (Table 1) was selected from the existing database, representing the diversity of structures and activity of the dataset. After alignment, the molecules were inserted as rows of a QSAR table along with their respective pK_i values. CoMFA steric and electrostatic fields were calculated as described below and entered as columns in the QSAR table. Standard steric and electrostatic CoMFA field energies of each compound were calculated using an sp³ probe atom with a charge of +1.0 |e| at all intersections in regularly spaced (2.0 Å) grids surrounding each inhibitor. Lennard-Jones and Coulomb-type potentials within the Tripos force field, and dielectric constant 1/r, were used in the calculation. The grid box dimensions were determined by the automatically created features in the CoMFA module within the SYBYL program. The same grid box was used in all calculations. An energy cutoff of 30 kcal/mol for both the steric and electrostatic contributions was set as threshold, and the electrostatic terms were dropped within regions of steric maximum, that is, 30 kcal/mol. Sixteen additional inhibitors (Table 1) were selected as a predictive set to test the robustness of the resulting model.

Partial least square analysis and statistical validation

The CoMFA descriptor served as independent variables and pK_i values as dependent variables in partial least square (PLS) regression analysis for deducing the 3D-QSAR models. PLS is an extension of multiple regression analysis in which the original variables are replaced by a small set of their linear combinations^37,38. These latent variables generated are used for multivariate regression, maximizing the commonality of explanatory and response variable blocks. PLS has been useful in cases in which the number of descriptors is greater than the number of samples as is the case with CoMFA.

The predictive value of the models was evaluated first by leave-one-out (LOO) cross-validation process^39,40. The cross-validated coefficient, q² (also called r²_cv), was calculated using Equation (1):

1

Where y_i, ŷ_i and are the observed, predicted and mean values of the target property (pK_i), respectively, for the training set. is the predictive residual sum of squares (PRESS). The optimal number of components obtained from the cross-validated PLS analysis were used to derive the final QSAR model using the compounds in the training set without cross-validation. Then, a non-cross-validation analysis was carried out, and the Pearson coefficient (r²) and standard error of estimates (SEE) were calculated. Finally, the CoMFA results were graphically represented by field contour maps, where the coefficients were generated using the field type ‘Stdev*Coeff’.

The ref. 24 has reported that the low value of q² for the training set can, although, serve as an indicator of a low predictive ability of a model, the opposite is not necessarily true. That is, the high q² is necessary, but not sufficient, for a model with the high predictive power. Thus, the external validation is the only way to establish a reliable QSAR model. Therefore, in the present work, for the external validation based on the validation set, the following criteria were used:

2 3 4 5

In Equation (2), the predictive r², r²_pred, is defined as follows:

6

where SD is the sum of the squared deviations between the experimental activity of the compounds in the test set and the mean activity of the compounds in the training set, and ‘PRESS’ is the sum of the squared deviations between predicted and experimental activity for every compound in the test set. In Equation (3), r²_test is the conventional correlation coefficient between experimental values and model predictions in the test set. The r²_o is a quantity characterizing linear regression with the Y-intercept set to zero (i.e. described by Y = kX, where Y and X are the actual and predictive activity, respectively). The k in Equation (5) is the slope of regression lines (predicted versus observed activities) through the origin. The definitions of the aforementioned statistical indices are presented in detail in references^24–27.

MD simulations

To identify a functionally validated complex from protein docking and the most potent molecule 69, we performed 4 ns MD simulations to investigate the conformational changes in the complex induced by the ligand 69. The software AMBER 10⁴¹ was used for the MD simulations. Partial charges for the inhibitors were calculated with AM1⁴² in MOPAC7 and fitted into the AM1-BCC type of charge⁴³. The force field parameters for these molecules were assigned by the Antechamber tool⁴⁴ in AMBER10. Hydrogen atoms were added to the protein with xleap module from AMBER. With the minimum distance from the surface of the complex to the faces of the box set to 12 Å, the system was then put in to a rectangular box of TIP3P water molecules⁴⁵, and this solvated system contained ~52,000 atoms, 16,000 water molecules. The final box dimensions were 640,000 ų. MD simulation was carried out by the SANDER module and a continuous system was simulated by the use of periodic boundary condition. The SHAKE algorithm⁴⁶ was applied to fix all bond lengths involving hydrogen bonds, permitting a 2-fs time step. The particle-mesh Ewald (PME) method⁴⁷ was preformed to handle Coulombic interactions, and a 8 Å cutoff was applied on all vander Waals interactions. After 1000 steps minimization and equilibration for 500 ps, the final 4 ns MD simulations were carried out at 300 K for 4 ns, using Berendsen temperature coupling⁴⁸ and constant pressure (1 atm) with isotropic molecule-based scaling.

Results and discussion

Docking reproduction of the crystal complex structure

To validate the ability of the docking protocol to reproduce a FIXa-binding structure of the crystal complex, the docking simulation was conducted against the crystal complex. First, the crystal complex of compound IZX/FIXa (PDB ID: 3LC5) was used to validate the docking protocol. The IZX (i.e. compound 37 in the present work) was docked back into the binding pocket of this enzyme using the protocol described in the Materials and methods section. As shown in Figure 1, the alignment of the redocked IZX and the crystal complex matches each other very well, showing that they take almost the same binding position in the active sites. The root-mean-square deviation (RMSD) between the conformations of the co-crystallized and the re-docked IZX is 0.43 Å, suggesting that Surflex-dock, in this work, is able to generate the experimentally observed binding mode for FIXa inhibitors and that the parameters set for the Surflex-dock simulation are appropriate for reproduction of the X-ray structure. Moreover, the moderate correlation coefficient between the total scores predicted by Surflex-dock and the activity is 0.4 (data not shown) indicating further the reliability of the simulation parameters. Thus, the rest series of benzothiophenes could be docked into the binding pocket of the FIXa (PDB ID: 3LC5) using the current docking protocols.

Figure 1.  Conformational comparison of IZX (coloured by atom type) from docking and the original ligand (magenta-coloured) in the crystal structure (PDB ID: 3LC5).

Molecular docking

The crystal structure of IZX with human FIXa shows some characteristics of the binding pocket¹² divided into four parts as shown in Figure 2. The pocket sandwiches the inhibitor between Gln192, Cys191, Ser190, Asp189 on the left and Trp215, Gly216, Phe174 on the right. There were two hydrogen bonds involved: Asp189 and Ser190 with the bound ligand IZX in the crystal structure 3LC5 at the P1 region, where residues were within 3 Å around the inhibitor. To be specific, the amidine group on the benzothiophene of the inhibitor forms two hydrogen bonds with the O atoms of Asp189 and Ser190, with bond lengths of 2.70 and 2.63 Å, respectively. A hydrophobic cavity consisting of side chains of Cys191, Gln192 on the left and Trp215, Gly216 and Phe174 on the right accommodates both the benzene rings in the regions P3 and P4, respectively, and the benzothiophene ring in the region P2 interacts well with these hydrophobic residues within 4 Å around the ligand IZX (Figure 2).

Figure 2.  The binding models of IZX (i.e. compound 37 in the present work) with FIXa (PDB ID: 3LC5). The inhibitor is represented as stick model and carbon atoms are coloured cyan by element; P1, P2, P3 and P4 refer to the corresponding positions in the binding pocket. All hydrogen atoms are removed for clarity.

The docking results could illustrate us the interaction modes between inhibitors and the serine protease enzyme (FIXa). In the present work, the most potent compound 69 is taken as the representative to elucidate this point in detail. As illustrated in Figure 3, compound 69 locates at the active site of the enzyme and binds primarily through the hydrophobic and polar interactions. And the majority of hydrogen bond (H-bond) and hydrophobic modes are the same to those observed in the crystal structure of IZX complexed with FIXa (PDB ID: 3LC5), which validates the reliability of the docking model.

Figure 3.  The binding models of the most potent ligand 69 (K_i = 0.002 μM) with FIXa (PDB ID: 3LC5). The inhibitor is represented as stick model and carbon atoms are coloured cyan by element; the key residues are denoted as stick and coloured green by element. The yellow dash lines denote the hydrogen bonds. P1, P2, P3 and P4 refer to the corresponding positions in the binding pocket. Only the hydrogen bonds within 3 Å are shown and all hydrogen atoms are removed for clarity.

The hydrophobic cavity for the interaction of compound 69 with FIXa, as seen from the figure, consists of the side chains of His147, Arg143, Cys191, Cys220 and Gln192 on the left, and Ile213, Ser195, His57 on the bottom, as well as Tyr99, Ile176, Thr175, Trp215, Lys98, Gly226, Gly216 and Phe174 on the right, which interact with the aromatic rings and aliphatic chain in the regions P2, P3 and P4, respectively (Figure 3). Compared with the interactions between IZX and FIXa (Figure 2), the most potent compound 69 forms two additional hydrogen bonds with the oxygen atoms of Asp189 and Gly216 with bond length of 2.9 and 2.8 Å, respectively. Additionally, the aliphatic tail linked to the benzene ring in the region P4 produces the hydrophobic interactions with the residues of Ile176, Lys98, Trp215, Tyr99, Phe174 and Thr175. It could be observed that the amidine group on the benzothiophene template (P1) is essential for FIXa inhibition, which might play a central role in fixing the orientation of the ligand in the pocket⁹ forming polar interactions with Asp189 and Ser190. In addition, by analysing this series of compounds docked in the serine protease enzyme (PDB ID: 3LC5), we observed that favoured substituents in P3 region are actually all aromatic groups due to the strong π–π interactions among the aromatic group and the aromatic rings of His147. In fact, the most potent compound, 69, docked in the FIXa shows a distance of 3.2 Å. This might be a possible reason that compound 39 with a phenyl group in this position exhibits over 10-fold FIXa activity than the corresponding methyl analogue 38.

Through the investigation of the series of benzothiophenes, an interesting common characteristics was observed for most of the inhibitors that they all have a low molecular weight with the rather flat, hydrophobic heterocyclic rings and distal amidine group, which makes their binding pocket recognitions comparable even if these inhibitors present different chemical structures. Figure 4 depicts the binding conformations of all the 84 benzothiophenes in the binding pocket of the FIXa, which are derived from the docking simulations using the Surflex-dock module of SYBYL software.

Figure 4.  The binding conformations of the benzothiophene analogues displayed inside the active site of the FIXa (PDB ID: 3LC5). The FIXa surface is rendered with electrostatic potential. P1, P2, P3 and P4 refer to the corresponding positions.

CoMFA results

Sixty-eight FIXa inhibitors were selected as the training set for constructing the CoMFA model, the remaining 16 ones served as the test set for the validation of the model. PLS analysis was carried out for the 68 binding conformations and the results are listed in Table 2, where the CoMFA model exhibits a cross-validated q² of 0.52 with six components. The non-cross-validated PLS analysis reveals a conventional r² value of 0.97, F = 381.44, and an estimated standard error of 0.20. The steric field descriptors explain 39.8% of the variance, whereas the electrostatic descriptors explain 60.2%. The fitted activity for the 68 inhibitors versus their experimental values with their residuals is listed in Table 1. The scatter plot between the fitted and the experimental activity for the CoMFA model of FIXa inhibition derived from the training set is depicted in Figure 5. All the results demonstrate that the CoMFA model is fairly fitted.

Table 2.  Summary of CoMFA analysis.^a

Parameters	CoMFA
q^2	0.52
SEP	0.86
PCs	6
SEE	0.20
r^2	0.97
F	381.44
Fraction
S	0.398
E	0.602
r^2pred	0.81
r^2test	0.90

^aq², cross-validation correlation derived from SAMPLES; SEP, standard error of prediction; PCs, principal components; SEE, non-cross-validated standard error of estimate; r², regression coefficient in the training set; F, Fisher’s F-value; S, steric field; E, electrostatic field; r²_pred, predictive r² in the test set; r²_test, conventional correlation coefficient in the test set.

Figure 5.  Fitted predictions versus actual pK_i for the CoMFA model of FIXa inhibition derived from the training set.

CoMFA contour

The CoMFA result is usually represented as 3D coefficient contour. It shows regions where variations of the steric and electrostatic nature in the structural features of the different molecules contained in the training set cause the increase or decrease in the activity. Figure 6 shows steric and electrostatic field contour plots of the model. In addition, the 3D-QSAR model was mapped back to the binding site of the FIXa, to get a better understanding of the vital interactions between the benzothiophene derivatives and the protease. To aid in visualization, the most potent ligand 69 is displayed in the map as a reference structure. Green-coloured regions indicate areas where increased steric bulk is associated with enhanced FIXa inhibitory activity and yellow regions suggest areas where increased steric bulk is unfavourable for the activity. Blue polyhedrons show areas where positive charges are favoured to FIXa inhibitory activity. On the contrary, red polyhedrons show areas where negative charge would benefit the activity.

Figure 6.  CoMFA StDev*Coeff contour plots combined with the most potent ligand 69. (A) Steric contour map. Green contours indicate regions where bulky groups increase the activity (favoured level 80%); yellow contours indicate regions where bulky groups decrease the activity (disfavoured level 20%). (B) Electrostatic contour map. Red contours indicate regions where negative charges increase the activity (favoured level 20%); blue contours indicate regions where positive charges increase the activity (disfavoured level 80%). The inhibitor is shown in magenta stick.

In Figure 6A, the common scaffold benzothiophene is sandwiched by several yellow-coloured contours in P2 region indicating that bulk groups are disfavoured in this position. This is in agreement with our docking result that the residues Gln192, Cys191 on the left and Trp215, Gly216 on the right, just as two walls, clamp the benzothiophene core and, in this way, produce strong hydrophobic interactions among them. This might be one of the reasons why the researchers tend to select this kind of rigid template as the core structure^1,8–10. A small yellow-coloured contour in the middle of the two walls (P2) suggests that the group with a big bulk in this position is disfavoured. This conclusion is well-illustrated by compounds 82–84 (Table 1), where the introduction of an additional fluoro group at the 6-position of the benzothiophene ring led to a marginal reduction of the FIXa activity, compared with the corresponding nonfluorinated analogues 58, 67 and 70, respectively. In addition, our docking results indicate that there is a change of ligand binding vector probably due to the repulsive force between the additional fluoro substituent and electron-rich−OH of Ser195, which also might be one of the reasons for the slight loss of activity. This, to some degree, is in agreement with the conclusion suggested by Wang et al.¹².

Three yellow-coloured contour maps can be observed at the P1 position suggesting that small groups in these regions benefit the biological activity, which observation is consistent with our docking result. In these regions, the hydrogen bonds can be formed between the residues of Asp189 and Ser190 with the amidine group on the benzothiophene template with very close distance, and any larger substitute group might lead to a collision with these residues in the pocket, therefore bulk groups here are forbidden. The yellow belt contours shown at the 4′-position of the benzene ring of Figure 6A at the P4 region represent that they are disfavoured regions for the steric interactions. For example, compared with compound 64 with the−OMe group, the piperazine-substituted analogue 66 exhibits 2-fold reduced potency probably due to the presence of the larger substituent in position 4′ of the benzene ring in the region P4. It can be seen that compound 60 exhibits relative high potency, whereas compound 72 just gives slightly low activity probably because of the existence of additional−F group at the 4′-position compared with compound 60. On the contrary, at the 2′-position of Figure 6A in P4 region, there is a big green contour, where, at the same position of FIXa, the residues Trp215, Tyr99, Lys98, Ile176 and Phe174 are located, which are characterized as a hydrophobic pocket. This also explains why compound 69 shows high activity, for the big−CH₂NHiBu substituent at this position interacting with these hydrophobic residues well (Table 1). Another green-coloured contour locating the flexible chain in the region P4 indicates that bulk groups may increase the potency. A close investigation on analogues 2–5 finds out that the 4-position-substituted groups well-locate at these positions. The compound 4 with a large substituent exhibits relative high activity among them. In addition, an increase of the length of this flexible chain could cause the distal substituent extending to the hydrophobic pocket, which forms hydrophobic interactions with residues Trp215, Tyr99, Lys98 and Phe174. This is in agreement with the experimental results. Compounds 56–59 with long-chain substituents exhibit more potential than compound 55 that possesses only a small−OH group at the same region. It is noteworthy that a green-coloured contour embedded in the bottom-left yellow maps suggests that the group with an appropriate bulk benefits the inhibitory activity of FIXa, which explains the experimental results that the appropriate CO₂Me substitute on R (compound 33) is better than other extended ester analogues 34 and 35 (Table 1). Exactly, the inhibitory activity of ligands 33 (substituted by−CO₂Me group), 34 (–CH₂CO₂Me) and 35 (–(CH₂)₂CO₂Me) is in an order of 33 > 34 > 35, which agrees well with the size order of the substitute on R (Table 1). In addition, an investigation of the inactive compound 26, with the large−PhCH₂ substituent positions the forbidden yellow region, which might cause its less potency.

In the electrostatic contour map (Figure 6B), areas where electronegative groups increasing the potency are represented in red, whereas areas where electronegative groups decreasing the potency are represented in blue. Using derivative 69 as reference, a close inspection of the electrostatic contour plot reveals that electropositive groups are more preferred surrounding N atoms. At the P1 region, three blue-coloured contours surrounding the amine groups indicate the favour for the electropositive substituents at the position. This observation is well-supported by our docking studies. As seen from Figure 3, the amidine group on the benzothiophene template in the reference molecule 69 forms several H-bonds with the oxygen atoms of residues Asp189 and Ser190, respectively. These H-bonds at the P1 position could be found in all these benzothiophene analogues, which may be one classical characteristics for almost all FIXa inhibitors (including not only the benzothiophene-based ones) with an amidine tail in the P1 position. This investigation is consistent with the previous reports^1,11,12 and the X-ray result (PDB ID: 3LC5), and recent studies have demonstrated that in the absence of the amidine, for any of the serine proteases, they are completely lack of potency due to the loss of key interactions with Asp189 in the S1 pocket⁸. In P4 region, two big blue-coloured regions exist suggesting the possibility of increasing the potency by an introduction of electropositive substituent at the region. Also, it can be observed from Figure 3 that a strong H-bond could be formed between the−NH substituent and the O atom of the residue Gly216. The red contours in the vicinity of−NH in region P1 and at the upper surface of the O atom linked between the benzethiophene and the phenyl ring in region P4 as well as flying near the phenyl ring indicate that electron-rich groups are beneficial to the activity. For example, at the P4 region, the docking modelling shows that possible H-bond interaction between the carbamate carbonyl and the NH of Glu219 could be formed, which is significant because it is a glycine residue at the same position in FXa^49,50. Thus, for improving the selectivity of benzothiophenes, more flexible linked groups in this position are needed (compounds 68, 69, etc.). It can be seen in Figure 6A, in the same position, a green-coloured region locates there, supporting the experimental fact that a long flexible group here would increase the potency.

Validation of the 3D-QSAR models

Validation implies a quantitative assessment of the robustness and predictive power of a QSAR model. The predictive power of a model can be defined as its ability to predict accurately the modelled property (e.g. biological activity) of new compounds. The details about the validation criteria for QSAR models have already been discussed in many literatures^24–27 and thus are directly employed here. In the present work, 16 randomly selected compounds (Table 1) were used as the test set to verify the constructed CoMFA models. The calculated results are listed in Table 1 and displayed in Figure 7. The predicted pK_i values derived from the QSAR model are in good agreement with the experimental data with a statistically tolerable error range (Table 2). The CoMFA model has also passed successfully Tropsha’s recommended tests for predictive ability, where

Figure 7.  Predicted versus actual pK_i plot of the CoMFA model derived from the test set.

From the above results, it can be observed that the CoMFA model would be reliably used in new inhibitor design for developing drug against coagulability.

MD simulations

In the present work, with the docked FIXa–compound 69 complex structure, 4 ns MD simulations were carried out, with purpose to dynamically explore the conformation changes that take place in aqueous solution between FIXa and the ligand. Figure 8A depicts the RMSDs of the trajectory with respect to the initial structure of the complex with a range from 0.56 to 1.27 Å. Clearly, a metastable conformation occurs after 2.5 ns of simulation, when the RMSD of the docked complex reaches and remains at about 0.88 Å throughout. The average structure of ensemble for the whole 4 ns simulation and the docked structure are superimposed as shown in Figure 8B, where the red and green ribbons represent the initial structure for the docked complex and the average structure, respectively. The ligand 69 is represented in stick for the initial complex and line for the final average complex, respectively. Clearly, no significant difference is observed between the average structure extracted from MD simulations and the docked model of the complex, validating the reliability of the docking model. Thus the docked model, we hoped, would be of help for virtual screening of novel potent FIXa inhibitors.

Figure 8.  MD simulation result. (A) Plot of the root-mean-square deviation (RMSD) of docked complex versus the MD simulation time in the MD-simulated structures. (B) View of superimposed backbone atoms of the average structure of the whole 4 ns of the MD simulation (green) and the initial structure (red) for compound 69–FIXa complex (PDB ID: 3LC5). Compound 69 is represented in cyan stick for the initial complex and cyan line for the final average complex. All hydrogen atoms are removed for clarity.

Conclusions

In this study, using the alignment scheme generated from the docking study, a highly predictive CoMFA model is developed based on 84 benzothiophene analogues as FIXa inhibitors. A statistically satisfactory model is obtained, with LOO cross-validation q² = 0.52, conventional r² = 0.97 and r²_pred = 0.81, respectively. The binding conformations of the benzothiophenes are further determined and predicted by molecular docking. The consistency between the CoMFA field distributions, docking and MD simulations proves the robustness of the QSAR model.

Our conclusions are that: (1) benzothiophenes form hydrogen bonds to the backbones of amino acids Asp189 and Ser190 as the classical hinge interactions; (2) some implications are drawn to improve the activity and selectivity of benzothiophene as FIXa receptor inhibitors, including that the small electropositive group substituent at the P1 position and the groups with an appropriate volume at the C-6 of the benzothiophene ring in region P2, the introduction of aromatic groups at the P3 region and bulk hydrophobic substituent at the C-2 and small electronegative groups at the C-4 position of the benzene ring at the region P4 might be preferable to produce higher activity. These results demonstrate the power of combining docking/QSAR approaches to explore the probable binding conformations of the compounds at the active sites of the protein target, and further provide useful information to understand the structural and chemical features of FIXa in designing and finding new potential inhibitors.

Acknowledgements

We are very grateful to Prof. Ling Yang for access to SYBYL software.

Declaration of interest

This work is financially supported by the National Natural Science Foundation of China (Grant No. 10801025).