AutoGPA-Based 3 D-QSAR Modeling and Molecular Docking Study on Factor Xa Inhibitors as Anticoagulant Agents

The three-dimensional-quantitative structure activity relationship (3D-QSAR) studies were performed on a series of direct factor Xa (FXa) inhibitors using AutoGPA-based modeling method in this paper. A training set of 38 molecules and a test set containing 10 molecules were used to build the 3D-QSAR model and validate the derived model, respectively. The developed model with correlation coefficients (r) of 0.8564 and cross-validated correlation coefficients (q) of 0.6721 were validated by an external test set of 10 molecules with predicted correlation coefficient (rpred) of 0.6077. Docking study of FXa inhibitors and FXa active site was performed to check the induced pharmacophore query and comparative molecular field analysis (CoMFA) contour maps using MOE2012.10. It was proved to be coincidence with the interaction information between ligand and FXa active site and was rendered to provide a useful tool to improve FXa inhibitors.


Introduction
Thrombosis can occur in the arterial and the venous or arise from the heart and lead to most cardiovascular events.The pathological mechanisms and risk factors that may lead to arterial thrombosis and venous thrombosis ware reported and approaches for antithrombotic therapy were described.[1] Vitamin K antagonists (VKAs) and heparin are most widely used anticoagulants to date.VKAs are the only available orally active anticoagulants for long-term use which function by reducing the synthesis of the vitamin K-dependent coagulation factors while heparin is used for surgery because the ability to bind to the protein antithrombin and increase the capacity of it to inhibit factor Xa and thrombin.[2] However, both of them shared same disadvantages such as narrow therapeutic window, slow onset of action, uncontrollable therapeutic level and food and drug interaction.[3,4] Therefore drugs that can be given in fixed doses without coagulation monitoring are a mainstay of cardiovascular therapy.[5] As alternatives to VKAs, direct thrombin inhibitors (DTI) and direct factor Xa inhibitors ware proved to be effective at reducing arterial and venous thrombosis in patients with cardiovascular disease because it has become evident that inhibition of factor IIa (FIIa, or thrombin) and factor Xa (FXa) is expected to reduce or prevent clotting.[6] Thrombin is a serine protease, which converts fibrinogen to fibrin and activates factor V, factor VIII and factor IX. [7] FXa belongs to a family of trypsin-like serine proteases, and plays a central role in the production of thrombin and the activation cascade of the blood coagulation system, linking the extrinsic and intrinsic activation pathways.[8] Since the pivotal position of FXa in the coagulation cascade, direct FXa inhibitors are supposed to be efficacious and safe compared with thrombin inhibitors.[9] Hence, the FXa inhibitors gradually become the largest family of new anticoagulants for long-term use.FXa inhibitors has become direct focus by most pharmaceutical companies resulting in the exploration of a number of small molecules, orally active FXa inhibitors, such as rivaroxaban and apixaban.[10,11].
Although FXa inhibitors have been synthesized up in different structures to the present.[12][13][14][15] it is time consuming to synthesis and assay every compound to find some initial hits.In the last decade, Computer-aided drug design (CADD) method including dock and pharmacophore search has been applied on search and design new FXa inhibitors.Recently, 3D-QSAR has been applied on the FXa inhibitors discovery in order to reveal the binding information between the 3D structure of the molecules and biological activity values.Currently, comparative molecular field analysis (CoMFA) [16] and comparative molecular similarity indices analysis (CoMSIA) [17] are the most successful and widely used 3D-QSAR techniques in ligand-based drug design (LBDD).Steric and electrostatic requirements of a biologically active molecule are revealed in both CoMFA and CoMSIA, Additionally, CoMSIA also gives hydrogen bonding and hydrophobic fields.Besides, 3D pharmacophore query [18] was reported as an alternative modeling method using compounds containing different molecular scaffolds, and compounds with different structures and common chemical features can be found by ligand pharmacophore map while CoMFA and CoMSIA methods were restricted to use compounds based on common molecular scafflods.Several models [19,20] derived by CoMFA and pharmacophore query method were reported to understand the relationship between structures and FXa protein, but information was still in urgent need.

Materials and methods
A dataset of 48 structures reported by Chuen Chan, Nigel S. Watson, Robert J. Young, Savvas Kleanthous et al. [21][22][23][24] The Ki values were converted to pKi(-lgKi).The Ki values which were used to build 3D-QSAR model cover wide range from 0.3nm to 285nm.In this paper, an AutoGPA-based 3D-QSAR model was developed, which combined pharmacophore, QSAR equation (electrostatic coefficient and steric coefficient) and allowable space, using known FXa inhibitors in MOE 2012.10.The AutoGPA module in MOE2012.10 was firstly reported in 2012 and proved to be more accurate and predictive than CoMFA.The information revealed by AutoGPA was believed to be more comprehensive than that by CoMFA method.[25] Since the biological activity of compounds synthesized by chemists was measured by different methods, it was advisable to build our 3D-QSAR model with biological data assayed by same method.The inhibitor activities of these compounds were measured by using Rhodamine 110, bis-(CBZ-glycylglycyl)-L-arginine amide as the fluorogenic substrate.The selected compounds were divided into training set (38 molecules) and test set (10 molecules), the former is to generate 3D-QSAR model and the latter to valid the derived model.Structures and corresponding pKi of training set and test set are shown in Table 1 and Table 2.In order to improve the quality of the model, it was necessary to consider that the training set should consist of highly active, moderately active and less active and the activities of test set should also in a wild range.The process of the model is shown in Figure 1.All of our studies including energy minimization, conformational searching and 3D-QSAR modeling were performed by MOE 2012.10 on windows 7 operating system.All structures of compounds were built and energy minimized under Merck molecular force field 94x (MMFF94x).Conformational searching for each compound was generated by conformation import method, which calculates low-energy conformations of molecules by decomposing each ligand into constituent overlapping fragment, performing stochastic conformational search on each fragment, and then assembling the fragments into unique conformers.Finally, a total of 5909 and 2454 conformations for training set and test set were generated respectively, and stored in respective database.The each generated conformation of training set was assigned to the typical pharmacophoric features such as hydrogen bond acceptor, hydrogen bond donor, positively and negatively ionized areas and aromatic center.Common pharmacophore queries were extracted from training set, in which some or most compounds were satisfied the queries, and most of conformers of compounds satisfied the four-point or five-point queries, resulting in a total of 33 pharmacophore queries, then all of these compounds were aligned to the generated pharmacophore queries above and the database contained the pharmacophore query information was stored for further AutoGPA-based 3D-QSAR model construction.Moreover, AutoGPA module improved 3D-QSAR model with allowance space around the aligned compounds makes it convenient to screen large molecule library.The AutoGPA-based 3D-QSAR model employed the algorithm of the CoMFA, as a consequence, all aligned molecules were placed in 3D grid box with a 2.0Å separation, following the calculation of electrostatic and steric interaction energies between the molecule and a series of probe atom (an sp3 carbon atom with charge +1.0) at each grip point which was to evaluated molecular fields.After all the grid potential fields were calculated, the partial-least-squares (PLS) analysis was introduced to study the relationship between grid potential fields and biological values to develop AutoGPA-based 3D-QSAR model.The optimal number of compounds up to ten is identified by leave-one-out (LOO) cross-validation and ranked by q2.The calculation method of q2 is shown below: q 2 = 1-MSE/YVAR (1) where: MSE means mean squared error, YVAR is the variance of experimental data.
The best q 2 value was saved as an AutoGPA model.In order to check the robustness, accuracy, reliability and predictability of the developed model, internal validation and external validation were employed.The LOO crossvalidation is used as internal validation while the prediction of test set as external validation.The validation was judged by r 2 and was calculated by a following equation: r 2 = (SD-PRESS)/SD (2) where: SD is the sum of squared deviations of biological activity from their mean, and predictive sum of squares (PRESS) is the sum of the squared differences between the actual and predicted biological activity.

Results
Ten different AutoGPA models were generated and the results were summarized in the Table 3.It was ranked by q 2 and q 2 >0.5 suggested that this model is reasonable and should have a good predictive ability.Thus model 1 was considered as the best model as it has highest q 2 (0.6721) value.Moreover, the correlation coefficient r 2 between observed and predicted activity was recommended to be another criteria for a given QSAR model.In this model, the r 2 of training set calculated by non cross-validation and LOO cross-validation was 0.8564 and 0.6726 respectively.External validation of test set was employed to test the predictability of this model, and the r 2 was 0.6077 which proved that the derived QSAR model is acceptable.Result of the most satisfying pharmacophore query and corresponding allowance space is shown in Figure 2. Four pharmacophores have been generated with three hydrophobiccentroids (green) and one projected acceptor (blue).The projected acceptor was located along implicit hydrogen direction and were used to annotate the location of possible hydrogen bond.Namely, an H-bond donor might located here and could form an H-bond with the ligand in a special direction.The three hydrophobic centroids indicated that the corresponding portions of ligand should decorated with hydrophobic ring or atom.The allowance space was introduced as a portion of 3D-QSAR model for the first time and helped to reject dropped conformers and evaluate the activity value in ligand screen procedure cooperating with pharmacophore query.An allowable space was suggested to be coincidence with sharp of active pocket and actually the allowance space derived by model was coincidence with the sharp of active pocket described in references.However, it was noted that no pharmacophore was generated in the p4 motif of the ligands.representing electron donating and electron withdrawing groups would be favorable respectively.A large region of red contours around the p4 motif of the ligand indicated that the presence of an electronegative substituent such as oxygen is very important for decreasing activity (pKi) while an electropositive substituent such as nitrogen increases activity.For example, the activity value of compound 1 with morpholino amide orienting to the red region is -0.7782 while compound 33 with benzylic amine having a nitrogen atom pointed to the red region shows an enhanced activity of -0.3010.
In the steric contour map as shown in Figure 3b, green contours showed well tolerance to bulky groups, whereas yellow contours implied regions of unfavorable steric effects.The two larger green contours around p1 motif of the ligand indicated that steric features have a positive effect on activity.For example, compound 13 with 6-indoly converted to compound 12 with 6chloronaphthyl, meanwhile, the activity value increased

Discussion
Docking study was carried out to better understand the developed model.2y7x.pdb was downloaded via www.rcsb.organd loaded into MOE2012.10.Structure preparation was performed to correct errors in the PDB file, followed by protonate 3D to assign ionization states and position hydrogens in a macromolecular structure.Water molecules have been retained because of its possibility of binding the ligand and imparting the binding affinity.The ligand of 2y7x.pdb was separated, stored in database file and went through conformation import procedure to generate conformations.Redocking was performed with default options to verify whether the position of ligand in pdb file is right and the docked conformation in the active pocket of FXa protein and ligand interactions were shown in Figure 4a and Figure 4b respectively.The active site of FXa protein is comprised of S1 pocket and S4 pocket.S1 pocket is deep, narrow and mainly formed by Tyr228 on the bottom, Asp189 and Gly226 on the side.S4 pocket is aromatic box which lined by aromatic side chains of Tyr99, Phe174 and Trp215.Present research has proved that when a p4 motif embodies a cation atom, commonly nitrogen atom under protonation state, can form an arene-cation interaction to improve binding affinity, [26,27] thus the p4 motif interacted with any aromatic side chain of Tyr99, Phe174 and Trp215 can improve binding affinity resulted in the increase of Ki value.Therefore, it made the p4 motif to be various, which explained the reason why there is no pharmacophore generated in the p4 motif.Furthermore, two water molecules existed in S4 pocket were bind to the nitrogen atom with charge +1 of p4 motif to form two strong H-bonds, where corresponded to the large region of red contours in the electrostatic contour map of the model.As we all see, three blue sections were placed around the 3-aminopyrrolidin-2-one means ligand exposure, so a hydrophobic scaffold was needed to enhance the binding affinity.Therefore, a hydrophobic center was emerged in the pharmacophore query.Moreover, an H-bond between oxygen of the ligand and Cys220 was observed which further corroborated the evidence that a projected acceptor in active site exists and interacts with ligand through an H-bond indicated by pharmacophore query.In addition, steric contours of the model was also proved referring to the docked conformations of compound 24 and 23 by the above docking method in Figure 4c and 4d.Chlorine was well pointed towards the center of the Tyr228 aromatic ring in compound 24 owing to the replacement of 6indoly(compound 23) by 6-chloronaphthyl, which is the so-called chloro binding mode.[28] In conclusion, despite that different types of 3D-QSAR models of FXa inhibitors have been reported by using SYBYL and Discovery Studio, a new type of 3D-QSAR model has developed by using AutoGPA module in MOE2012.10, which was based on CoMFA method and combined pharmacophore and allowance space in our study.The model showed good prediction and reliability with q2 (0.6721) and r2 (0.8564).In addition, docking study was carried out to check the reliability of the electrostatic and steric contour map.The results of docking and the developed 3D-QSAR model were complementarity.Thus this model provided a reliable tool for screening lead compound and guiding the design of novel and potent FXa inhibitors.

DOI: 10
.1051/ C Owned by the authors, published by EDP Sciences, 201

Figure 1 .
Figure 1.The process of AutoGPA model.

Figure 2 .
Figure 2. Pharmacophore query and allowance space of the developed model.Green ball means hydrophobic centroid (Hyd) indicating that a hydrophobic ring or atom may exist.Projected acceptor (Acc2) is represented by blue ball and stands for an acceptor forming an H-bond with ligand.

Figure
Figure3ashows the electrostatic contour map in the developed model, as indicated by red and blue contours

Figure 3 .
Figure 3. Electrostatic and steric contours map from AutoGPAbased model, compound 12 is shown as reference (A) Electrostatic contour map, red contour suggests electropositive charge and red contour electronegative charge.(B) Steric contour map, green represents bulky groups are favorable to increase activity and yellow indicates decrease activity.

Figure 4 .
Figure 4. Docking study of active pocket in 2y7x.pdb.(A) Redocked conformation of ligand separated from 2y7x.pdb showing a well fitness for active pocket, colored by blue.(B) Ligand interactions graph drawn with MOE2012.10.(C) P1 motif of compound 24 colored by purple from side aspect.Chlorine was pointed to the center of the Tyr228 aromatic ring.(D) P1 motif of compound 23 colored by blue from side aspect cannot well fitted into S1 pocket with 6-indoly.

Table 1
Structures and pKi of Compounds of Training Set and Predicted pKi of the Model No .

Table 2
Structures and pKi of Compounds of Test Set and the Predicted pKi of the Model

Table 3
Ten AutoGPA models based on the training set.