Discovery of Potential Orthosteric and Allosteric Antagonists of P2Y1R from Chinese Herbs by Molecular Simulation Methods

P2Y1 receptor (P2Y1R), which belongs to G protein-coupled receptors (GPCRs), is an important target in ADP-induced platelet aggregation. The crystal structure of P2Y1R has been solved recently, which revealed orthosteric and allosteric ligand-binding sites with the details of ligand-protein binding modes. And it suggests that P2Y1R antagonists, which recognize two distinct sites, could potentially provide an efficacious and safe antithrombotic profile. In present paper, 2D similarity search, pharmacophore based screening, and molecular docking were used to explore the potential natural P2Y1R antagonists. 2D similarity search was used to classify orthosteric and allosteric antagonists of P2Y1R. Based on the result, pharmacophore models were constructed and validated by the test set. Optimal models were selected to discover potential P2Y1R antagonists of orthosteric and allosteric sites from Traditional Chinese Medicine (TCM). And the hits were filtered by Lipinski's rule. Then molecular docking was used to refine the results of pharmacophore based screening and analyze the binding mode of the hits and P2Y1R. Finally, two orthosteric and one allosteric potential compounds were obtained, which might be used in future P2Y1R antagonists design. This work provides a reliable guide for discovering natural P2Y1R antagonists acting on two distinct sites from TCM.


Introduction
Human P2Y receptors are a family of nucleotide activated G protein-coupled receptors (GPCRs) which comprises eight subtypes [1]. According to the difference of coupling protein, eight subtypes are subdivided into two groups, named Gqcoupled P2Y1R and Gi-coupled P2Y12R. Gq-coupled P2Y1R (P2Y1, P2Y2, P2Y4, P2Y6, and P2Y11) can activate phospholipase C-(PLC ) and promote intracellular calcium mobilization, while Gi-coupled P2Y12R (P2Y12, P2Y13, and P2Y14) can inhibit the adenylyl cyclase and activate certain ion channels [2]. Among the two groups, P2Y1R and P2Y12R play a significant role in ADP-induced platelet aggregation. The blockade of either receptor could effectively inhibit ADP-induced platelet aggregation and decreases thrombosis formation [3,4]. Nowadays, antithrombotic drugs mainly act on P2Y12R but have some potential side effects, such as a long bleeding time and dyspnea [5][6][7]. Comparing with P2Y12R inhibitors, P2Y1R antagonists may offer an advantage of reducing bleeding liabilities [7]. Therefore, P2Y1R is becoming an especially attractive antithrombotic drug target. While still in the development stage, the antithrombotic drugs of P2Y1R antagonists are necessary for treatment thrombosis.
Consistent with all of the GPCRs, P2Y1R possesses the characteristics of the seven-transmembrane domain (7TMD). And it has two distinct sites, including orthosteric and allosteric site [8]. Orthosteric site located within the 7TMD, while allosteric site is situated at the lipidic interface of the 7TMD. Previous study [9] showed that orthosteric antagonists can completely block the ADP-mediated platelet aggregation and effectively reduce the formation of arterial thrombosis with only a slight prolongation of bleeding time. And allosteric antagonists [10] will accelerate the dissociation of its endogenous agonist 2MeSADP and substantially reduce platelet aggregation. It can also significantly control 2 Evidence-Based Complementary and Alternative Medicine bleeding time with a limited effect on bleeding. Fortunately, the complex of the P2Y1R with orthosteric and allosteric antagonists was revealed, respectively, by Dr. Zhang et al. [8], which can help us understand the mechanisms of inhibiting ADP-induced platelet aggregation and provides a basis for discovering potential target compounds with less side effects and higher security. However, the crystal structure of P2Y1R with allosteric antagonists BPTU was discovered recently. Before then, some of compounds with similar biological activity as allosteric antagonists were regarded as orthosteric antagonists, which will have an adverse impact effect on the virtual screening based on ligand. Therefore, it is a primary problem to classify antagonists bound with orthosteric and allosteric site. In this paper, 2D similarity search was used to distinguish antagonists of two distinct sites. Ligand-based similarity search is one of the common computed techniques to find similar structure of molecules based on 2D fingerprint [11]. The searching principle states that molecules with similar structures are likely to have similar biological activity. In this process, template molecule with certain biological activity was chosen as a query to search for similar molecules and the top-ranked molecules are likely to exhibit the required bioactivity [12][13][14]. Pharmacophore based screening and molecular docking would also be used in combination to discover potential P2Y1 antagonists acting on orthosteric and allosteric sites. Ligand-based pharmacophore models [15] can be established by extracting common chemical features that are responsible for their specific bioactivity, and the optimal model will be used to screen compounds which are arranged in the same relative orientation with the same features. Structure-based molecular docking is mainly employed to refine the results of pharmacophore based screening and analyze the interactions between proteins and small molecules.
Some constituents in Chinese herbal medicine have proven effectiveness in inhibiting ADP-induced platelet aggregation. For example, ginsenoside Rg2 exhibited stronger anticoagulation effects [16]. Therefore, it is possible to discover potential P2Y1 antagonists from Chinese herbs. The purpose of this study was to screen potential P2Y1 antagonists from Traditional Chinese Medicine Database (TCMD, version 2009) by using a series of molecular simulation methods. 2D similarity search was utilized to identify allosteric antagonists and orthosteric antagonists; Genetic Algorithm with Linear Assignment of Hypermolecular Alignment of Datasets (GALAHAD) was used to separately construct pharmacophore models of orthosteric and allosteric antagonists. Validated by built-in parameters and test set, two optimal pharmacophore models were selected as templates to search potential P2Y1 antagonists from the TCMD. Then molecular docking was utilized to refine the results of pharmacophore and analyze the receptor-ligand interactions. The results of virtual screening provide a basis for the discovery of potential antithrombotic lead compound acting on P2Y1R.  [17] was used to measure the similarity by comparing the selected fingerprint property. Theoretically, the range of Tanimoto coefficient value is zero to one. The value is closer to one indicating the greater similarity between the two fingerprints of compound and reference compound. By entering "human P2Y1 receptor antagonists" as a search term in the Binding Database (http://www.bindingdb .org/), 357 antagonists were obtained. MRS2500 and BPTU, two initial compounds which, respectively, bound with orthosteric and allosteric site, were selected as queries to search similar structures from 357 antagonists, respectively. All the parameters were set to the default value automatically.

Materials and Methods
Then 357 antagonists were successfully divided into two parts by 2D similarity search. The 2D similarity search results of orthosteric and allosteric compounds were listed in Tables 1 and 2, respectively. There is no industry standard for selecting the threshold of Tanimoto coefficient. For example, 0.3 and 0.5 were set as threshold in some literatures to identify compounds [18,19]. From Tables 1 and 2, compounds with Tanimoto coefficient value higher than 0.7 can account for more than 90% of the total compounds. False-positive rate may be reduced by a relatively higher similarity threshold. Hence, 0.7 was selected as the threshold of Tanimoto coefficient to classify orthosteric and allosteric antagonists from all compounds. Then 55 orthosteric antagonists and 278 allosteric antagonists were obtained and there is no duplicated compound between the two parts. To insure the comparability of the data, 246 compounds with Ki value were selected from 278 allosteric antagonists. And those compounds will be employed to build train set and test set.

GALAHAD Pharmacophore Hypotheses Generation.
The compounds which obtained by 2D similarity search had similar structures. Therefore, GALAHAD was selected as the pharmacophore modeling method, which is a typical modeling method based on compounds with similar structure. Based on the results of 2D similarity search, 7 P2Y1R orthosteric antagonists were selected to build orthosteric pharmacophore hypotheses. The remaining orthosteric antagonists were defined as active compounds of test set. Similarly, 6 P2Y1R allosteric antagonists were selected to generate allosteric pharmacophore hypotheses, and the remaining allosteric compounds were defined as active compounds of test set. The chemical information of the orthosteric compounds and allosteric compounds in the training sets was shown in Figures 1 and 2, respectively. All the compounds were added into Tripos force field and then were minimized by Powell method with 1000 iterations. Then GALAHAD models were constructed based on the train set. Initially, the compounds in the training set aligned to each other in torsional space. And a genetic algorithm (GA) was operated to identify a set of compounds conformations with minimized strain energy and maximized pharmacophoric similarities. Then, the optimal set of conformations was aligned in Cartesian space as rigid-bodies. During the GA run, pharmacophore and steric quartets were employed to evaluate overlap. Twenty models were generated with different features, conformations, and overlay of the molecules [20,21].

Validation of the Pharmacophore Model.
Internal parameters were calculated to evaluate the generated pharmacophore models firstly. The models were selected by three criteria [22]: (1) number of "hits" ( Hits) should approximately equal the number of compounds which are used to generate the pharmacophore models; (2) the model needs to have lower Energy; (3) the models should have higher Specificity. The models which satisfied all the criteria would be selected and evaluated by the test set using UNITY module.
One test set including 48 active compounds and 144 inactive compounds was used to evaluate the orthosteric pharmacophore models. And the other test set was used to evaluate the allosteric models which contained 240 active and 720 inactive compounds. Both of the inactive compounds were selected from the Binding Database. Part of the structures of compounds in the testing sets and the active values are shown in Figures 3 and 4.
The external parameters are presented as follows [23,24]: HRA (the effectively hit ratio of active compounds), IEI (identified effective index), and CAI (comprehensive appraisal index). Considering all factors, the optimal pharmacophore model of orthosteric and allosteric antagonists was obtained, respectively.

Database
Search. The optimal pharmacophore model of orthosteric and allosteric antagonists was utilized as queries to search the potential P2Y1R orthosteric and allosteric antagonists from TCMD. The "flexible database search" was carried out to perform the virtual screening process. Then, the hit compounds were filtered by Lipinski's rule of five (≥4), including MWT ≤ 500, A Log ≤ 5, H-bond donors ≤ 5, and H-bond acceptors ≤ 10 [25]. Compounds which meet the requirements were remained. Then, two lists of compounds including orthosteric and allosteric antagonists, with druglike properties, were obtained. Finally, the two lists would be further analyzed in molecular docking study, respectively.  algorithm with the smallest RMSD was selected for further employing.

Molecular Docking
In addition, in order to further validate the rationality of pharmacophore model and active pocket, the initial compounds were used to match the optimal model and analyze the interactions with active pocket of P2Y1R. And then compounds hit by two optimal pharmacophore models were docking into the crystal structure. Finally, potential P2Y1 orthosteric and allosteric antagonists which got higher docking score and formed favorable interaction with amino acid residue were obtained, respectively.

Pharmacophore Model Studies
3.1.1. Pharmacophore Model Studies of Orthosteric Antagonist. Twenty models were produced by GALAHAD module based on a training set including seven active orthosteric compounds. Internal parameters such as Hits, Specificity, Energy, and Pareto Ranking were used to evaluate the models based on three criteria mentioned in Section 2.1.3. Firstly, ten models with values of " hits" more than 5 were displayed in Table 3, which indicated that at least six active compounds Evidence-Based Complementary and Alternative Medicine 7 Moreover, Pareto rank of these models was zero, which suggested that those models were not superior to each other. The four pharmacophore models were validated by the test set of orthosteric antagonists; the results of validation were shown in Table 4. The HRA values of all those models were 100%, which indicated that all the models have the best ability to identify active compounds from test set. What is more, Model O-01 achieved the highest IEI and CAI, which indicated that Model O-01 had the best ability to identify active compounds from the inactive compounds. Thus, Model O-01 was chosen as the optimal pharmacophore model of orthosteric antagonists to screen the TCMD. The optimal model was shown in Figure 5 which included thirteen features: three hydrogen bond donors (DA 1, DA 2, and DA 9), seven hydrogen bond acceptors (AA 3, AA 4, AA 5, AA 6, AA 10, AA 11, and AA 12), and three hydrophobic features (HY 7, HY 8, and HY 13). Among them, DA 2 and AA 3 were generated in the same position, while AA 6 and DA 9 were also generated in the same position.

Pharmacophore Model Studies of Allosteric Antagonist.
Based on a training set including six allosteric antagonists, twenty GALAHAD models were derived. Similar to validation method of P2Y1R orthosteric antagonist, internal parameters were utilized to evaluate the models firstly. Six models were selected based on three criteria mentioned in Section 2.1.3, and then the models were validated by the test set of allosteric antagonist constructed previously.
The results of validated by test database were shown in Table 5. Model A-17 has highest value of HRA, which indicated that it has the strong ability to identify active compounds. However it has extremely low value of IEI, which suggests that it has poor ability to distinguish active compounds from inactive ones. Therefore, it cannot be regarded as the optimal model. The rest of five models with approximate values of HRA and the ability to identify active compounds are similar. However, Model A-07 has the highest value of IEI and CAI, which suggested that Model A-07 has the best ability to identify active compounds from the inactive compounds. Thus, Model A-07 was regarded as the optimal pharmacophore model of allosteric antagonists for screening the TCMD.
The optimal model was shown in Figure 6, which consisted of ten features, including three hydrogen bond donors (DA 1, DA 6, and DA 7), three hydrogen bond acceptors (AA 2, AA 3, and AA 8), and four hydrophobic features (HY 4, HY 5, HY 9, and HY 10). In Model A-07, a green sphere covered a purple sphere because DA 1 and AA 2 were generated in the same position.

Database Search.
The Model O-01 of P2Y1R orthosteric antagonists and Model A-07 of P2Y1R allosteric antagonists were treated as the optimal pharmacophore models which were employed to screen P2Y1R antagonists from TCMD, respectively. Then, two lists were obtained including 3480 and 6229 compounds, respectively. During this process, the QFIT value was computed to assess the matching degree between compounds and pharmacophore models, and a higher value of QFIT indicated that corresponding compound was mapped well with the pharmacophore model. Then further filtering based on "Lipinski's rule of five" was implemented; 798 and 961 potential drug-like P2Y1R orthosteric antagonists and allosteric antagonists were obtained, respectively. Finally, the antagonists would be docked into corresponding active pockets of P2Y1R, respectively, by using molecular docking algorithm.

Orthosteric Antagonists Molecular
Docking. The orthosteric active pocket was created with a sphere radius of 9.50Å around the initial compound MRS2500 in 4XNW. Then, the MRS2500 was redocked into the crystal structure by using LibDock and CDOCKER, respectively, and the RMSD values of corresponding two docking algorithms were 3.81Å and 0.58Å. Therefore, CDOCKER with the value of RMSD less than 2.00Å was chosen for docking study of orthosteric antagonists based on the smaller RMSD.
In addition, the initial compound, MRS2500, was used to match the optimal pharmacophore model (Model O-01) and analyze the interactions with orthosteric active pocket of P2Y1R. The mapping results of MRS2500 with Model O-01 and the interactions results between MRS2500 and amino acid residue were shown in Figure 7. The initial compound MRS2500 was mapped well with all the features of pharmacophore model. The purine ring which mapped features of HY 7 and HY 13 also formed hydrophobic interactions with LEU44 and ARG287. And the phosphate group formed hydrogen bond interactions with TYR306, THR205, and ASP204, which were also mapped with features of DA 2, AA 3, and AA 12, respectively. Thus, the analysis results of pharmacophore modeling and molecular docking were almost consistent which suggested that the study is reliably.
Then, 798 drug-like compounds which were filtered by the optimal pharmacophore model were docked into  a Ha is the number of active compounds screened by a pharmacophore model. b Ht is the total number of compounds screened by a pharmacophore model. c HRA represents the ability to identify active compounds from the test set. d IEI represents the ability to identify active compounds from nonactive compounds. e CAI is the comprehensive appraisal index.  the orthosteric active pocket, resulting in a hit list of 285 compounds. The compounds docking score was expressed by -CDOCKER ENERGY, and a high -CDOCKER ENERGY corresponds to a favorable binding mode between compound and P2Y1R. And the compounds, which had higher docking score and QFIT, were selected as potential compounds for further research. Finally, among the potential compounds, (2R)-hydroxy-4-(9-adenyl)butyric acid and ganoderpurine were regarded as potential P2Y1R orthosteric antagonists. The two candidates have got higher -CDOCKER ENERGY score, QFIT value, and similar interaction modes with P2Y1R-MRS2500. In addition, Lentinus edodes, the source herbal medicine of (2R)-hydroxy-4-(9-adenyl)butyric acid, proved to be effective in inhibiting the aggregation of platelets [27]. More specifically, (2R)-hydroxy-4-(9-adenyl)butyric acid got a higher QFIT 81.24 and was mapped with eight features of the optimal pharmacophore model. The compound got a higher -CDOCKER ENERGY score 30.96 and formed hydrogen bond interactions with ASP204, THR205, ASN283, and TYR306 and also formed hydrophobic interactions with LEU44 and ARG287. What is more, when the purine ring of the compound was mapped with features of HY 7 and HY 13, the amino acids residues LEU44 and ARG287 also formed hydrogen bond interactions with the compounds. In terms of hydrophobic interactions, the nitrogen atom formed hydrogen bonds with TYR303 and was also mapped with features of AA 4 (Figure 8(a)).
Ganoderpurine well mapped seven features with the Model O-01 and the QFIT value was 90.19. Besides, ganoderpurine achieved a higher -CDOCKER ENERGY score 26.48. It also formed hydrogen bond interactions with TYR303 and formed hydrophobic interactions with TYR203. Furthermore, the atoms which formed hydrogen bonds and hydrophobic interactions with TYR303 and TYR203 were also mapped with AA 4 and HY 8, respectively (Figure 8(b)). Therefore, the result of pharmacophore modeling study agrees with the molecular docking study, which suggested that the above results are credible. Besides, two compounds, eckol and nodakenin which can inhibit the aggregations of platelets [28,29], were hit by present study. And the source herbal medicines of the other four compounds have been found pharmacological effects of inhibiting the aggregation of platelets [27,[30][31][32]. The four compounds and source herbal medicines were as follows: icariol A1 (Epimedium sagittatum), lentysine (Lentinus edodes), boeravinone E (Boerhavia diffusa), and 6-demethoxycapillarisin (Artemisia capillaris-derived). To some extent, this further proved the reliability of the results of present study.

Allosteric Antagonists Molecular Docking.
Similarly, the initial compound BPTU was employed to define the active pocket of allosteric site of P2Y1R. And the radius of the binding pocket was 10.16Å. The RMSD was calculated after the initial compound BPTU was redocked into the crystal structure using LibDock and CDOCKER, respectively, and the corresponding RMSD values were 0.74Å and 0.61Å. In that case, both docking algorithms could reproduce the P2Y1-BPTU binding model reasonably. However, the closer the RMSD is to zero, the more reliable the results of the corresponding docking procedure is. Thus, CDOCKER was selected for further docking study of allosteric antagonists.
To contrast result of pharmacophore model and docking, the initial compound of allosteric site, BPTU, was used to match the optimal pharmacophore model (Model A-07) and analyze the interactions with allosteric active pocket of P2Y1R. The mapping results of BPTU with Model A-07 and the interactions results between BPTU and allosteric binding site were shown in Figure 9. From Figure 9, it found that the initial compound BPTU was mapped well with all the features of the Model A-07. Hydrogen bonding and hydrophobic interactions were the major interactions between ligand and receptor. To be specific, BPTU formed hydrogen bond interactions with LEU102 and hydrophobic interactions with PHE62, PRO105, LEU126, and MET123, which was consistent with literature [8]. What is more, pyridyl group mapped with HY 4 also formed hydrophobic interactions with ALA106 and PHE119, while the benzene rings which formed hydrophobic interactions with LEU126, MET123, and PHE62 were also mapped with the corresponding hydrophobic features. In conclusion, the analysis results of pharmacophore modeling and molecular docking were almost consistent, and this indicated the rationality of the result.
In the following study, 961 drug-like compounds obtained from the result of pharmacophore model were docked into the allosteric active pocket. Similar to orthosteric antagonists, the docking score was expressed by -CDOCKER ENERGY. Taking value of -CDOCKER ENERGY and QFIT into consideration, compounds with high-scoring function values were retained. And these compounds were treated as potential P2Y1R allosteric antagonists. Finally, among this potential compounds, desmodilactone was regarded as most promising P2Y1R allosteric antagonist. For instance, desmodilactone was mapped with seven features of Model A-07. And the QFIT value and -CDOCKER ENERGY score were 50.91 and 15.98, respectively. Then from the analysis of the key amino acids, desmodilactone formed hydrogen bond interactions with LEU102 and formed hydrophobic interactions with PHE62 and PRO105. Furthermore, the amino acids residue LEU102 formed hydrogen bonds with the compound, and the oxygen atom of the compound has also formed a hydrogen interaction with the same amino acid ( Figure 10). Moreover, the result of molecular docking and pharmacophore was almost uniform, which indicated that the above results are reasonable.
Previous studies [33] showed that rhinacanthin Q hit by virtual screening has the effect on inhibiting the platelet aggregation. In addition, herbal medicines of other four potential compounds were usually used to treat Blood Stasis Syndrome based on Chinese medicine theory. The compounds and herbal medicines were as follows: 6-methylgingediol (Zingiber officinale), dehydrohirsutanonol (Viscum cruciatum), 5-O-caffeoyl quinic acid butyl ester (Erigeron breviscapus), and Delamide (Aconitum carmichaeli). Furthermore, the former two herbal medicines, Zingiber officinale [34] and dehydrohirsutanonol have proven [35] to be the effective Chinese medicine in inhibiting the aggregation of platelets. On the basis of the above analysis, the results of pharmacophore and molecular docking are reliable.

Conclusions
In recent years, thrombosis has been a serious threat to human health. P2Y1R antagonists, whether recognizing orthosteric or allosteric sites, have a unique advantage in clinical safety of medication. In the paper, series technologies of computer aided drug design were jointly used, including 2D similarity search, pharmacophore modeling, and molecule docking. 2D similarity search firstly attempted to solve the problem of inappropriate classification of P2Y1R antagonists. Then, 55 orthosteric antagonists and 246 allosteric antagonists were obtained. In the following study, GALAHAD models of orthosteric and allosteric antagonists of P2Y1R were generated separately. And Model O-01 and Model A-07 were treated as the optimal pharmacophore model, respectively, to retrieve potential compounds from TCMD. Finally, crystal structures of P2Y1R were used to refine the hit compounds and analyze the ligand interaction mode between hits and P2Y1R.
In conclusion, the study successfully attempted to use the crystal structure of P2Y1R and series of computational approaches to discovery of potential P2Y1R antagonists bound with two distinct sites from Chinese herbs. Finally, three potential P2Y1R antagonist leads were obtained. To be specific, (2R)-hydroxy-4-(9-adenyl)butyric acid and ganoderpurine were regarded as promising P2Y1R orthosteric antagonists, while desmodilactone was treated as the most potential allosteric P2Y1R antagonist. The result was expected to be further employed in the search for potential natural P2Y1R orthosteric and allosteric antagonist. The shortcoming of this study was without combining the biological experiments. In the future, molecular simulation methods and biological experiments will be combination used to discover potential novel antithrombotic drugs.