Discovery of novel and potent P2Y14R antagonists via structure-based virtual screening for the treatment of acute gouty arthritis.

P2Y14 nucleotide receptor is a Gi protein-coupled receptor, which is widely involved in physiological and pathologic events. Although several P2Y14R antagonists have been developed thus far, few have successfully been developed into a therapeutic drug. In this study, on the basis of two P2Y14R homology models, Glide docking-based virtual screening (VS) strategy was employed for finding potent P2Y14R antagonists with novel chemical architectures. A total of 19 structurally diverse compounds identified by VS and drug-like properties testing were set to experimental testing. 10 of them showed good inhibitory effects against the P2Y14R (IC50 < 50 nM), including four compounds (compounds 8, 10, 18 and 19) with IC50 value below 10 nM. The best VS hit, compound 8 exhibited the best antagonistic activity, with IC50 value of 2.46 nM. More importantly, compound 8 restrained monosodium uric acid (MSU)-induced pyroptosis of THP-1 cells through blocking the activation of Nod-like receptor 3 (NLRP3) inflammasome, which was attributed to its inhibitory effects on P2Y14R-cAMP pathways. The key favorable residues uncovered using MM/GBSA binding free energy calculations/decompositions were detected and discussed. These findings suggest that the compound 8 can be used as a good lead compound for further optimization to obtain more promising P2Y14R antagonists for the treatment of acute gouty arthritis.

Introduction

The P2Y14 receptor (P2Y14R) is a member of P2-purigenic receptors, which has been regarded as inhibitory adenylate cyclae G-protein (Gi)-coupled receptor. It inhibits the production of 3′,5′-cyclicadenosine monophosphate (cAMP) through Gi protein, which could be activated by endogenous uidine diphosphate (UDP)-sugars. Activation of P2Y14R has been regarded to be associated with proinflammatory reactions, leading to neutrophil chemotaxis and mast cell degranulation [1], [2], [3], [4]. P2Y14R is distributed among a variety of immune cells and is expressed in extensive tissues [5], [6], [7].

Several animal studies have demonstrated the value of P2Y14R as potential therapeutic target for recruitment of macrophages to liver, induction of insulin resistance in diabetes and local inflammation [8], [9], [10]. However, there are few studies focused on relationship between P2Y14R and acute gouty arthritis, which is a group of characteristic inflammatory reactions caused by innate immune disorders. Acute gouty arthritis is triggered by deposition of monosodium urate crystals (MSU) in the joint, resulting from the activation of Nod-like receptor 3 (NLRP3) inflammasome [11], [12], [13].

Our recent studies have showed that inhibition of NLRP3-mediated pyroptosis is a viable strategy for the prevention and treatment of acute gouty arthritis [14], [15]. Till now, the treatment of gout still lacks the ideal drug. Previous study suggested that MSU can induce high expression of P2Y14R in human keratinocytes [16], offering strong evidence that P2Y14R might play causal role in MSU-related diseases. Meanwhile, the activation of P2Y14R is closely related to the content of intracellular cAMP, which was demonstrated to negatively regulate NLRP3 inflammasome [17], involved in inflammatory, diabetes, immune processes and other related complications [18], [19]. Therefore, P2Y14R is likely to regulate the inflammatory response through NLRP3 inflammasome via cAMP in acute gouty arthritis.

To date, the current researches on P2Y14R antagonists only reported three types of compounds including pyrimidine piperidine, 2-naphthoic acid and 3-substituted benzoic acid [7], [9], [20], [21], [22]. Among them, the most active and selective P2Y14R antagonist is (4-(piperidin-4-yl)-phenyl)-7-(4-(trifluoromethyl)-phenyl)-2-naphthoic acid (PPTN, IC50 = 4 nM). However, the currently reported antagonists represented by 2-naphthoic acid suffered from poor solubility, low oral bioavailability, and high difficulty in synthesizing raw materials, bringing greater difficulties to further discussion of structure-activity relationship and biological evaluation [22]. In addition, based on P2Y14R homology models, a novel P2Y14R antagonists with scaffold, 3-(4-phenyl-1H-1,2,3-triazol-1-yl)-5-phenyl substituted benzoic acid was reported by Jacobson and co-workers using molecular docking and molecular dynamics (MD) simulation approaches. The identified P2Y14R antagonists showed quite acceptable binding affinities and the IC50 value of most potent P2Y14R antagonist was 31.7 nM. Based on these observations, there remains ongoing need to explore potent P2Y14R antagonists with novel chemical architectures. Besides, the development of promising P2Y14R antagonists could be a reasonable way for the treatment of gout.

Due to the high cost and time-consuming of high-throughput screening (HTS), virtual screening (VS) has aroused widespread concerns and been widely used in lead compound identifications of drug discovery [23], [24]. For the crystal structures of P2Y14R has not yet been reported, the structure-based virtual screening (SBVS) can be used for finding novel P2Y14R antagonists with diverse chemical scaffolds based on well-established homology modes of P2Y14R [20], [21], [22].

To our knowledge, this is the first case to carry out a molecular docking strategy to massively screen a commercial library for finding novel P2Y14R antagonists based on P2Y14R homology models. Two well-prepared and minimized P2Y14R homology models (HM1 and HM2) [21] were used to screen the ChemDiv database. 19 diverse compounds were selected using drug-likeness properties prediction, REOS filtering, core scaffold clustering and purchased for biological testing. 10 of them (VS hit rate > 50%) exhibited significant antagonistic activity against P2Y14R (IC50 < 50 nM) and the most potent lead, compound 8 displayed a quite satisfactory antagonistic activity with IC50 value of 2.46 nM. Then, the feasibility of compound 8 as a drug candidate for treating gout treatment was investigated through a series of pharmacodynamics and mechanism of action. The results demonstrated that compound 8 restrained MSU-induced pyroptosis of THP-1 cells through blocking the activation of NLRP3 inflammasome, which was attributed to its inhibitory effects on P2Y14R-cAMP pathways. Finally, the Molecular Mechanics/Generalized Born Surface Area (MM/GBSA) binding free energy calculations/decompositions were employed to preliminarily detect the interaction patterns between P2Y14R and two most potent hits (compounds 8 and 18). The key favorable residues for P2Y14R antagonists binding were detected and discussed. These findings may guide us to discovery more promising P2Y14R antagonists for treating acute gouty arthritis in the near future.

Materials and methods

P2Y14R homology models for docking-based virtual screening

The P2Y14R homology models (HM1 and HM2) [21] well established by Trujillo et al. were selected, optimized and applied in the Glide docking-based VS campaign of Schrödinger 9.0 software [25]. By utilizing the Protein Preparation Wizard module of Schrödinger 9.0, all water molecules were removed, the broken side chains were repaired and missing hydrogen atoms were added. Then, using the OPLS2005 force field, the partial charges and protonation states were assigned for each homology model.

Molecular docking-based virtual screening procedure

First of all, the Receptor Grid Generation module of Glide of Schrödinger 9.0 was used to generate binding site/pocket for molecular docking. The binding pocket size was set to 10 Å × 10 Å × 10 Å and centered on the centroid of the ligand in each P2Y14R homology model.

Then, the ChemDiv library including more than 2 million compounds was selected as screening database and screened against two P2Y14R homology models. Using the LigPrep mode of Glide, all compounds in the ChemDiv database were preprocessed carefully. For each compound in ChemDiv, the tautomers were generated at pH = 7.0 ± 2.0 and the different combinations of chiralities were also generated by setting the maximum number of stereoisomers to 32 by using Epik. At last, the final well-prepared ChemDiv database comprising more than 2.6 million compounds was set to Glide docking-based VS pipeline.

P2Y14R inhibitory activities screening

HEK293 cell lines stably expressing the P2Y14R were purchased from Keygen Biotech Co, ltd. Cells were plated in 384-well plates approximately 24 h before the assay at the density of 10,000 cells per well. Before assay, cells were briefly washed with phosphate-buffered saline solution to remove traces of serum and then incubated with 7.5 μL induction buffer contained 30 μM Forskolin (Med Chem Express, Cat. #HY-15371), 10 μM UDP-glucose (Sigma Aldrich, Cat. # U4625) and various concentrations of test compounds (0.01 nM, 0.1 nM, 1 nM, 10 nM, 100 nM) for 30 min at 37 °C, each concentration of 3 repetitions. P2Y14R inhibitory activities at each concentration were evaluated by detecting cAMP levels in order to calculate IC50 values.

Cell culture

THP-1 cell line purchased from American Type Culture Collection (Manassas, VA, USA) was cultured and stimulated with phorbol 12-myristate 13-acetate (PMA) as previous studies. Then, cells were pre-treated with Compound 8 or PPTN for 1 h, followed by the stimulation with MSU (500 μg/ml) for 12 h. Subsequently, the culture supernatants were collected for further investigation.

Measurements of IL-1β and cAMP

IL-1β concentrations in culture supernatants and cAMP levels in cell lysis were detected with ELISA Kit (Neobioscience, Shenzhen, China) or cAMP-Glo™ Assay Kit (Promega, WI, USA).

Pyroptosis assay

For pyroptosis analysis, active Caspase-1 and PI fluorescence of samples were measured using flow cytometry. Active caspase-1 was detected with FLICA 660 Caspase-1 Detection Kit (Immuno Chemistry Technologoes, USA), and propidium iodide (PI) staining was used to assess the integrity of cellular membrane.

Immunofluorescence

After MSU stimulation, the cells were 4% paraformaldehyde fixed for 20–30 min. Permeabilization was performed with 0.3%–0.5% Triton X-100 for 20–30 min. When blocking for 1 h to avoid non-specific protein interactions, the samples were incubated with the primary antibody and secondary antibodies in sequence as previous studies. Fluorescent images were visualized by confocal laser scanning microscope (Fluoview, FV1000, Olympus, Japan).

Western blot

The THP-1 cells collected from each group were lysed in a RIPA buffer (Sigma, St. Louis, MO, USA). Samples containing approximately 50 mg protein was separated by 8–12% SDS-PAGE followed by the transference to polyvinylidene fluoride membranes (Millipore Corporation, MA, USA). Subsequently, PVDF membranes were treated with primary antibodies overnight at 4 °C after being blocked. The membranes were washed three times with Tris buffer saline-Tween20 (TBST), followed by incubation with appropriate horseradish peroxidase-conjugated secondary antibodies for 2 h. Finally, protein bands were visualized with an enhanced chemiluminescence (ECL) system (Keygen Biotech, China) and scanned with a Chemiluminescence gel imaging system (Tanon-5200Multi, China).

Statistical analysis

The data are expressed as mean values ± SDs. Data analyses were performed by one-way ANOVA with Tukey multiple comparison test (Graphpad Prism 7.0a), with p < 0.05 considered as significant.

Results and discussion

Molecular docking-based virtual screening pipeline

Two well-established P2Y14R homology models (HM1 and HM2) [21] proposed by Trujillo et al. were selected and minimized for the following docking-based virtual screening pipeline (Fig. 1). Three scoring functions of Glide docking (HTVS, SP XP) were applied to perform the sequential VS strategy [25]. The 50,000 highest -ranked compounds of the prepared ChemDiv database predicted by HTVS were re-docked using SP scoring mode. Then, the 5000 highest-ranked compounds of SP were re-calculated using the XP scoring function. At last, 1000 highest -ranked compounds were obtained for each P2Y14R homology model. Followed by removing duplicates, Lipinski “Rule-of Five” filter [26] and drug-likeness models built in our previous studies [27], [28], [29], [30], the compounds with reactive, undesirable functional groups or toxic were also deleted by applying REOS criterion [31]. Then, the compounds with less than two chiral centers were retained and then the remaining compounds were clustered using the Tanimoto coefficient evaluated based on MACCS structural keys (Tanimoto coefficient cut off value = 0.7). At last, 19 compounds were selected from ChemDiv database and purchased for experimental testing (Table 1).

Fig. 1

The predicted binding poses and interaction patterns of (a) homology model 1 (HM1) and (b) homology model 2 (HM2) of P2Y14R. (The co-ligands in HM1 and HM2 are UDP-[1] glucose and UDP, respectively).

Table 1

Biological activities, representative molecular properties and key parameters identified in docking-based VS of the 19 purchased compounds from ChemDiv database.

Compd	ID_numbera	IC50(nM)	docking scoreb	MW	logP	logS	model
1	1278–0075	ND	−14.97	299.30	1.20	−2.29	HM1
2	1683–7093	35.4	−12.76	351.38	3.11	−4.46	HM1
3	2372–3546	ND	−16.34	303.29	2.18	−3.68	HM2
4	3473–2589	28.7	−14.39	349.43	3.54	−4.83	HM2
5	3975–0036	ND	−15.82	344.11	2.23	−4.00	HM2
6	4393–0019	45.5	−14.34	451.49	5.37	−6.67	HM1
7	5369–0063	16.5	−15.10	329.31	0.74	−2.13	HM2
8	6521–0066	2.46	−13.06	305.33	2.81	−3.79	HM1
9	7244–0067	18.6	−14.72	273.29	3.07	−2.40	HM2
10	8011–4760	5.35	−14.22	362.43	2.87	−5.92	HM1
11	8012–2120	ND	−12.98	389.42	3.77	−4.01	HM1
12	8013–6020	ND	−12.35	344.36	3.89	−4.31	HM1
13	8020–2337	ND	−14.06	231.17	1.21	−2.16	HM2
14	C301-4660	12.39	−15.15	387.46	2.50	−5.52	HM2
15	F293-0086	ND	−15.51	404.51	4.12	−4.81	HM2
16	F293-0205	ND	−14.64	438.53	4.51	−5.53	HM2
17	K783-4166	ND	−12.64	339.73	4.30	−5.62	HM1
18	Y040-3078	5.12	−13.60	314.34	3.90	−5.10	HM2
19	Y041-2308	7.71	−15.10	396.35	3.55	−5.82	HM1
20	PPTNc	2.74		475.50	3.80	−10.42

The compound number labeled in the ChemDiv database. According to the purity statements, the purity of all compounds purchased from the ChemDiv database is higher than 95%.

The predicted binding affinity for compounds using the XP function based on HM1 or HM2 homology models.

Positive control.

In vitro P2Y14R inhibitory activities screening

P2Y14R inhibitory activities of testing compounds were determined based on production of cAMP in a HEK293 cell line stably expressing P2Y14R. The results were listed in Table 1. As can be seen in Table 1, 10 of 19 purchased compounds (VS hit rate = 52.63%) showed quite acceptable inhibitory activity (IC50 < 50 nM) for P2Y14R. The chemical structures of 10 identified P2Y14R antagonists with IC50 value below 50 nM are shown in Fig. 2 and those of the remaining compounds were shown in the Fig. S1 in the Supporting Information. Among them, four compounds (compound 8, 10, 18 and 19) exhibited satisfactory antagonistic activity of P2Y14R (IC50 value below 10 nM, Fig. 3), and compound 8 showed the most potent antagonistic activity (IC50 = 2.46 nM). The schematic representations of the predicted binding poses and interaction patterns between the P2Y14R and the two most potent identified antagonists (compounds 8 and 18) are depicted in Fig. 4.

Fig. 2

Structures of 10 potent antagonists of P2Y14R (IC50 < 50 nM) identified from Glide docking-based VS.

Fig. 3

Fluorescent assay of P2Y14R binding affinities (IC50 curves) of four identified P2Y14R antagonists (compounds 8, 10, 18 and 19) with IC50 value below 10 nM, PPTN was run as positive control.

Fig. 4

(a) The predicted conformations of compounds 8 and 18 derived from Glide docking (the complexes of compound 8-P2Y14R was colored in golden and compound 18-P2Y14R was colored in green) and (b) predicted interaction patterns for compounds 8 and 18 in the binding pocket by applying HM1 and HM2 as docking structure, respectively. (For interpretation of the references to colour in this figure legend, the reader is referred to the web version of this article.)

In vitro anti-inflammatory effects of compound 8 through regulation of cAMP and NLRP3 inflammasome

As shown in Fig. 5, cAMP concentrations were significantly decreased after MSU stimulation, which was reversed by pre-treatment of compound 8 and PPTN. More importantly, MSU administration led to a significant increase in the proportion of pyroptotic cells characterized by Caspase-1/PI double positive staining analyzed by flow cytometry. As expected, this alternation was also improved in compound 8 and PPTN treated cells (Fig. 6a and b). Consistently, IL-1β levels in the supernatant of THP-1 cell culture medium were obviously increased in model group. Both compound 8 and PPTN interventions apparently inhibited the release of IL-1β, reflecting the mitigation of inflammation caused by MSU (Fig. 6c). As shown in Fig. 7, protein expressions of NLRP3, ASC (apoptosis-associated speck-like protein containing a CARD) and Caspase-1 p20 were apparently increased in THP-1 cells with MSU stimulation. And aforementioned alterations were reversed by pre-treatment of compound 8 and PPTN.

Fig. 5

Effects of compound 8 and PPTN on levels of cAMP in MSU-treated THP-1 cells. Compared with Control group: ###P < 0.001. Compared with Model group: *P < 0.05, **P < 0.01, ***P < 0.001. Each group (n = 4).

Fig. 6

Effects of compound 8 and PPTN on proportions of Caspase-1 single positive and Caspase-1/PI double positive cells (a and b), as well as levels of IL-1β (c) in cell culture supernatants of MSU-treated THP-1 cells. Compared with Control group: ###P < 0.001. Compared with Model group: *P < 0.05, **P < 0.01, ***P < 0.001. Each group (n = 4).

Fig. 7

Effects of compound 8 and PPTN on protein expressions of NLRP3, ASC and Caspase-1 (p20) (a and b) in MSU-treated THP-1 cells. Compared with Control group: ###P < 0.001. Compared with Model group: *P < 0.05, **P < 0.01, ***P < 0.001. Each group (n = 4). Representative confocal microscopy photographs of THP-1 cells with immunofluorescence changes are presented (c).

On the other hand, inhibitory effect of compound 8 on NLRP3 inflammasome was also confirmed by immunofluorescence data (Fig. 7c). When compared to control cells, model cells apparently showed higher fluorescence intensity in NLRP3 and ASC staining without observed difference in DAPI (4′,6-diamidino-2-phenylindole) intensity.

Primary structure-activity relationship discussions using MM/GBSA free energy decompositions

For exploring the detected antagonistic activity differences, the most potent VS hits (compounds 8 and 18) of P2Y14R were selected and docked into the respectively binding pocket of P2Y14R homology models (HM1 and HM2) using Glide XP scoring mode. In order to investigate the interaction pattern between PPTN and P2Y14 receptor, the PPTN was docked into the binding pocket of HM1 and HM2 using SP and XP scoring modes of Glide docking. The docking results demonstrated that PPTN cannot produce acceptable docking poses against P2Y14 receptor. Considering the higher protein flexibility of P2Y14 receptor, the PPTN may adopt quite distinct binding mode with P2Y14 receptor, compared with assayed compounds in our study.

By employing the MM/GBSA approach [32], [33], [34], the predicted binding poses of compounds 8 and 18 interacting with P2Y14R were optimized and rescored. The predicted total binding free energies using MM/GBSA rescoring of compounds 8 and 18 were −42.85 and −61.20 kcal/mol, respectively (Table 2). Then, for quantitatively discerning the contribution of each key residues of P2Y14R binding with compounds 8 and 18, the antagonist-residues interaction spectra were depicted and analyzed. As can be seen in Fig. 8a, two most potent antagonists of P2Y14R have quite distinct binding sites in the binding pocket of P2Y14R. For example, the residues ofVal99, Asn156, Cys172, Lys176, Arg253 and Gln260 play as the key residues for the compound 8 binding with P2Y14R, and their favorable contributions to the total binding energy (ΔGpred) are all lower than −1.5 kcal/mol. Compared with compound 8, compound 18 has quite different favorable binding residues. The dominant residues of compound 18 interacting with P2Y14R are Lys77, Ala98, Phe101, Arg253, Gln260 and Lys277. The same key residues for compounds 8 and 18 binding with P2Y14R are Arg253 and Gln260. The energy contributions of Arg253 and Gln260 for compounds 8 were −11.07 and −4.20 kcal/mol (Fig. 8a), and those for compound 18 were −10.30 and −3.82 kcal/mol (Fig. 8b), respectively. Considering inherent high flexibility of P2Y14R structure, we found that maintaining stable/strong interactions with these favorable residues (Lys77, Ala98, Val99, Phe101, Asn156, Cys172, Lys176, Arg253, Gln260 and Lys277) are the requirements for obtaining promising P2Y14R antagonists. This finding will provide some clues to design/develop more optimal antagonists of P2Y14R in the lead optimization stage.

Table 2

The predicted binding free energies using MM/GBSA rescoring of compounds 8 and 18.

Compd	Polar contributions		Nonpolar contributions		ΔGprede
	ΔEelea	ΔGGB1b	ΔEvdwc	ΔGSAd
8	−698.72	699.19	−36.92	−6.40	−42.85
18	−388.98	373.62	−39.38	−6.46	−61.20

Electrostatic contribution.

Polar part of desolvation.

Van der Waals contribution.

Non-polar part of desolvation.

The predicted total binding energies using MM/GBSA calculations.

Fig. 8

(a) The binding poses of compounds 8 and 18 optimized from the MM/GBSA calculations (the favorable residues for compounds 8 and 18 binding with P2Y14R are colored in golden and green, respectively. The same key residues for two compounds are colored in red), (b) the antagonist-residues interaction spectra of compounds 8 and 18.

Conclusions

In the current work, we adopted Glide docking-based virtual screening strategy for finding potent P2Y14R antagonists using two well-established P2Y14R homology models. 19 potential hits with quite novel chemical scaffolds were set to antagonistic activity testing. 10 of them revealed significant antagonistic activity against P2Y14R. The IC50 of the most potent identified P2Y14R antagonist (compound 8) can reach 2 nM, which was higher than the previously reported 2-naphthoic acid compound PPTN. To further confirm its feasibility as a drug for the prevention and treatment of acute gouty arthritis, we established a THP-1 cell model exposed to MSU to simulate acute gouty arthritis. The results demonstrated that compound 8 can significantly restore cAMP production and reduce IL-1β secretion. More importantly, compound 8 blocked the pyroptosis of THP-1 cells and inhibited the activation of NLRP3 inflamasome. These findings indicate that the compound 8 might be applied as a good lead compound for further modification/optimization for the treatment of acute gouty arthritis.

Compliance with ethics requirements

This article does not contain any studies with human or animal subjects.

Declaration of Competing Interest

The authors declared that they have no conflicts of interest to this work.

We declare that we do not have any commercial or associative interest that represents a conflict of interest in connection with the work submitted.