Structural Modeling and in Silico Screening of Potential Small-Molecule Allosteric Agonists of a Glucagon-like Peptide 1 Receptor.

The glucagon-like peptide 1 receptor (GLP-1R) belongs to the pharmaceutically important class B family of G-protein-coupled receptors (GPCRs), and its incretin peptide ligand GLP-1 analogs are adopted drugs for the treatment of type 2 diabetes. Despite remarkable antidiabetic effects, GLP-1 peptide-based drugs are limited by the need of injection. On the other hand, developing nonpeptidic small-molecule drugs targeting GLP-1R remains elusive. Here, we first constructed a three-dimensional structure model of the transmembrane (TM) domain of human GLP-1R using homology modeling and conformational sampling techniques. Next, a potential allosteric binding site on the TM domain was predicted computationally. In silico screening of druglike compounds against this predicted allosteric site has identified nine compounds as potential GLP-1R agonists. The independent agonistic activity of two compounds was subsequently confirmed using a cAMP response element-based luciferase reporting system. One compound was also shown to stimulate insulin secretion through in vitro assay. In addition, this compound synergized with GLP-1 to activate human GLP-1R. These results demonstrated that allosteric regulation potentially exists in GLP-1R and can be exploited for developing small-molecule agonists. The success of this work will help pave the way for small-molecule drug discovery targeting other class B GPCRs through allosteric regulations.

Introduction

In the past 40 years, the number of overweight people has increased sixfold worldwide with 33.6% of US men and women being obese, and this disturbing trend is projected to continue in the foreseeable future.[1] Obesity can cause a number of health problems including cardiovascular diseases and diabetes, with type 2 diabetes representing 90–95% of diabetes patients. The glucagon-like peptide 1 receptor (GLP-1R), a member of class B family of G-protein coupled receptors (GPCRs), is an effective target for the treatment of type 2 diabetes,[2] and its incretin peptide and varied peptide mimetics are adopted drugs.[3]

Despite remarkable antidiabetic effects, GLP-1 peptide-based agonists have several shortcomings.[4,5] They are available only in a format for injection, lacking effective long-term glucose control capability, and they can cause side effects and result in low quality of life in some patients. Hence, there is significant interest in the development of nonpeptidic small-molecule agonists of GLP-1R with enhanced bioavailability.[6−11] However, such a strategy remains problematic because of the nature of the orthosteric binding site for GLP-1 (Figure A), which is large and relatively shallow.[12] Up to now, no small-molecule drugs acting as GLP-1R agonists are available in the market. Therefore, novel approaches in developing small-molecule drugs targeting GLP-1R are very desirable for the treatment of type 2 diabetes.

Figure 1

Comparison of orthosteric and allosteric binding sites for class B GPCRs. (A) Orthosteric binding site of an endogenous peptide illustrated by the cryo-EM structure of GLP-1R (PDB ID: 5VAI). (B) Allosteric binding site of a small-molecule ligand illustrated by the crystal structure of CRF-1R (PDB ID: 4K5Y). The peptide and small-molecule ligands (colored in pink) were represented in a CPK model. For clarity, only the TM domain of the receptors was shown.

Given the allosteric nature of GPCRs, targeting the allosteric sites on GPCRs for small-molecule therapeutic intervention represents an alternative and promising approach for drug discovery.[13−16] Cinacalcet, a positive allosteric modulator of the calcium sensing receptor[17] and Maraviroc, a negative allosteric modulator of the chemokine receptor CCR5,[18] are the first two allosteric drugs of GPCRs entering the market. Compared to ligands acting at orthosteric sites, allosteric ligands demonstrate several potential benefits, including better subtype selectivity and reduced side effects,[19,20] as well as biased functional selectivity and novel modes of efficacy.[15,21] In addition, allosteric agonists may benefit the development of orally delivered peptide drugs targeting the orthosteric site by augmenting the efficacy of endogenous and exogenous peptide and its analogs. Recently, the existence of the allosteric site for class B GPCRs was confirmed by the reported crystal structure of the corticotropin-releasing factor receptor 1 (CRF-1R),[22] which showed an antagonist binding at an allosteric site deep inside the transmembrane (TM) domain of the receptor (Figure B). Hence, targeting the allosteric sites of GLP-1R for small-molecule drug discovery could represent a promising alternative for overcoming shortcomings related to GLP-1 peptide-based treatment.

High-throughput screenings have identified a few small-molecule allosteric agonists binding to the TM domain of GLP-1R, with some in covalent modification, for example, compound 2 and BETP.[6−8] However, the lack of a high-quality structure of the TM domain of GLP-1R until very recently[23−25] hindered further development of these lead compounds into drugs. In this work, we first performed homology modeling and conformational sampling to generate a three-dimensional (3D) structure model of the TM domain of GLP-1R. Then, we applied the rational ligand-based and structure-based drug design techniques to screen the ZINC database[26] for identification of druglike small-molecule as allosteric agonists of GLP-1R. Finally, the agonistic and allosteric effects of the top-ranked compounds were studied using a cAMP response element (CRE)-based luciferase reporting system and insulin assay. These results confirmed that allosteric regulation exists in GLP-1R and can be exploited for developing small-molecule agonists of GLP-1R.

Results

3D Model of TM Domain of GLP-1R

Two homology models of the TM domain of GLP-1R were constructed, one based on the template structure of CRF-1R (PDB ID: 4K5Y) and the other based on that of the glucagon receptor (GCGR) (PDB ID: 4L6R). Both models were then subjected to the primary enrichments. Results from preliminary enrichment studies indicated that the GLP-1R model built based on the CRF-1R structure showed better enrichment for the 23 GLP-1R ligands. Therefore, this GLP-1R model was selected for subsequent conformational sampling calculations.

Hundred conformations of the TM domain of GLP-1R were generated through conformational sampling calculations (Figure ). Among them, 10 distinct conformations were identified. For each of the 10 GLP-1R models, the ROC plot was generated. The best enrichment results were obtained for the Conformation #8 (Figure S1 online). About 70% of the 23 active compounds were present in the top 30% of the docking results. Therefore, Conformation #8 was selected for the virtual screening process.

Figure 2

Superimposition of 100 conformations of the TM domain of GLP-1R generated through conformational sampling.

Small-Molecule Agonists of GLP-1R Identified through in Silico Screening

For ligand-based screening, the 23 active compounds from PubChem[30] were used as the query individually to identify similar compounds in the small-molecule library of ZINC database.[26] A total of 5689 compounds that had similar ligand properties (molecular weight, x log P, hydrogen donors, hydrogen acceptor and polar surface area) with the 23 active compounds were identified. Next, those 5689 compounds were docked into the same predicted allosteric site of Conformation #8 respectively using Glide. Glide in the Schrodinger Suite is one of the best docking tools to identify potential protein effectors.[53] Each dock pose was then ranked based on its glide and docking scores. Top-ranked poses were visually inspected and they all bound to the proposed allosteric site in the TM domain (Figure ). Nine top-ranked compounds were identified and eight of them were purchased and experimentally tested for their potential activity against GLP-1R (Table ).

Figure 3

Docking pose based on the homology model of GLP-1R generated above and chemical structure of two small-molecule agonists of GLP-1R. (A) Compound M_4 (colored in blue) interacts with residues L360 and P358, which are indicated by dashed lines; (B) chemical structure of compound M_4; (C) compound M_9 (colored in pink) interacts with residues L244 and R190, which are indicated by dashed lines; and (D) chemical structure of compound M_9.

Table 1

Top-Ranked Compounds in the Predicted Allosteric Binding Site of the GLP-1R Model with Their Docking Scores

no.	ZINC ID	glide score (kcal/mol)	docking score (kcal/mol)	in vitro test
M_1	ZINC01008161	–9.234	–9.234	yes
M_2	ZINC08396882	–9.449	–9.137	yes
M_3	ZINC01017526	–9.094	–9.093	yes
M_4	ZINC00702587	–9.115	–9.088	yes
M_5	ZINC00678328	–8.957	–8.957	yes
M_6	ZINC08397024	–8.935	–8.935	yes
M_7	ZINC00793676	–8.872	–8.872	no
M_8	ZINC02146229	–8.872	–8.869	yes
M_9	ZINC08400241	–8.852	–8.851	yes

In Vitro Activity of the Top-Ranked Compounds

In vitro activity of the eight top-ranked compounds (Table ) from virtual screening was first studied using the GLP-1R-dependent luciferase reporter system. In this screening system, the rat GLP-1R plasmid was transiently transfected in HEK293 cells, stably expressing luciferase reporter gene (HEK-CREB luciferase). The sequence identity between human GLP-1R and rat GLP-1R is ∼90%, and there is no difference in the predicted allosteric site between the two sequences. Hence, the rat GLP-1R was chosen for screening compounds. The activation of the rat GLP-1R was measured as the amount of luminescence in response to cAMP, which in turn was normalized to the amount of protein. The relative luminescence change at each concentration of compound was calculated with respect to the vehicle control. These readings were used to plot the dose–response curves. Three independent experiments were performed for each compound, and negative control was included in all experiments to evaluate nonspecific effect of compounds (if any). From in vitro screening, two compounds (M_4 and M_9) were found to activate GLP-1R. The EC50 value of compounds M_4 and M_9 was 24 and 26 μM, respectively (Figure A,B).

Figure 4

In vitro agonistic activity of small-molecule GLP-1R agonists in HEK293 cells coexpressing GLP-1R or VIPR1 receptor and a 3x-CRE-luciferase reporter. (A) Dose–response curves of compound M_4 in the presence and absence of rat GLP-1R, respectively. (EC50 = 24 μM). HEK293-CREBluciferase cell line transiently expressing rat GLP-1R was treated with different concentrations of M_4. GLP-1R activation was measured as the amount of luminescence produced, which was normalized by protein concentration; (B) dose–response curves of compound M_9 in the presence and absence of rat GLP-1R, respectively. (EC50 = 26 μM). The HEK293-CREB luciferase cell line transiently expressing rat GLP-1R was treated with different concentrations of M_9. GLP-1R activation was measured as the amount of luminescence produced, which was normalized by protein concentration; (C) effect of GLP-1R agonist M_4 in the presence or absence of the VIPR peptide antagonist. The HEK293-CREB luciferase cell line was treated with different concentrations of M_4 in the presence or absence of the VIPR peptide antagonist (5.6 μM). The nonspecific effect of M_4 was measured as the amount of luminescence produced, which was normalized by protein concentration. In all experiments, normalized luminescence was plotted with respect to vehicle control [0.5% dimethyl sulfoxide (DMSO)]. For (A,B), the dose–response curve plotted for negative control showed some nonspecific effect. The difference between dose–response curves for HEK293-CREB luciferase cells with and without GLP-1R was significant (p < 0.0001). The dose–response curves were generated using a sigmoidal dose response (variable slope) from GraphPad Prism 6.0. Data in all three figures are representative of three independent experiments with at least three technical replicates for each treatment conditions, and error bars for each concentration were plotted as standard error of the mean (SEM; n = 3).

Besides GLP-1R, the HEK293 cells are known to express other functional class B GPCRs including vasoactive intestinal peptide receptor 1 (VIPR1).[48] Sequence alignment of residues in the predicted allosteric binding site showed high sequence identity between human GLP-1R and some of these class B receptors (Figure S2 online), implicating that compounds M_4 and M_9 could bind to other class B GPCRs. Therefore, the nonspecific luciferase activity by these two compounds was studied using HEK293-CREB luciferase cells transiently transfected with an empty vector. The results indicated that these compounds only exhibited limited nonspecific effect on HEK-CREB cells (Figure A,B).

Nonspecific stimulation of VIPR1 by compounds M_4 and M_9 could lead to side effects. To further assess this effect, the optimal concentration of a VIPR peptide antagonist that would inhibit VIPR1 activation present on HEK293_CREB cells in the presence of the VIPR peptide agonist was estimated (Figure S3 online). Then, the effect of M_4 on VIPR1 stimulation was studied in the presence or absence of the VIPR peptide antagonist. By comparing the changes of luciferase activity in the absence and presence of the VIPR peptide antagonist, the data indicated that M_4 slightly induced nonspecific VIPR1 activity (Figure C).

Compound M_4 Synergizes with GLP-1 To Activate Human GLP-1R

Low levels[38] and a decreased response of GLP-1 have been observed in some type 2 patients.[39−42] Therefore, it will be of interest to determine whether compounds M_4 and M_9 can act as an allosteric modulator of GLP-1R and enhance the bioactivity of endogenous GLP-1. The activation of GLP-1R by M_4 (24.71 μM) in combination with different concentrations of GLP-1 (0.014–1450 nM) was studied by luciferase activity responding to cAMP production using HEK293-CREB cells transiently expressing human GLP-1R. GLP-1R activity stimulated by GLP-1 in combination with M_4 (24.71 μM) was significantly increased than by using GLP-1 alone, and the synergistic effect was found to be dose-dependent (Figure A). This synergistic activity suggested that M_4 might act as an allosteric modulator of human GLP-1R.

Figure 5

Synergistic effect of agonist M_4 on GLP-1R and VIPR1 receptor. (A) Synergistic effect of M_4 in the presence of GLP-1 on GLP-1R; (B) synergistic effect of M_4 in the presence of VIP on VIPR1; GLP-1R and VIPR1 activation were assessed as luminescence normalized to protein concentration and plotted as luminescence fold change with respect to vehicle control (0.5% DMSO). Data are the average of three independent experiments with at least three technical replicates for each conditions, and error bars for each concentration were plotted as SEM (n = 3). Statistical analysis was done using 2-way ANOVA (****p < 0.0001; **p < 0.001).

Given the nonspecific effect of M_4 on VIPR1 present on the HEK293-CREB cell line (Figure C), it was necessary to examine the behavior of M_4 (24.71 μM) on VIPR1 in the presence of VIP (7.5–7510 nM), using the same cell line with overexpressed VIPR1. Only at the high VIP concentration, VIPR1 activity in response to VIP in combination with M_4 was higher than using VIP alone, indicating the potential nonspecific activity of M_4 on VIPR1 (Figure B). However, the amplitude of such a synergistic effect of M_4 was significantly lower than when M_4 was used with a broad range of GLP-1, suggesting that M_4 significantly improved the GLP-1R-mediated cAMP production compared with the VIPR1-mediated cAMP production. Overall, M_4 synergizes well with GLP-1 to stimulate GLP-1R activity, potentially by acting as an allosteric modulator of GLP-1R. Further work is certainly required to prove that M_4 acts on the allosteric site of GLP-1R.

M_4 Stimulates Insulin Secretion

The goal of this work is to develop small-molecule agonists of GLP-1R that will stimulate insulin production in pancreatic β cells. The insulin production activity of M_4 was assessed by in vitro insulin secretion assay in INS-1832/13 cells. For the assay, INS-1 832/13 cells were first starved with Krebs Ringer Bicarbonate (KRB) buffer. After starvation, cells were treated with KRB containing glucose (16.7 mM) and GLP-1 (181 nM) or M_4 (20 μM) for 10–20 min. After the treatment, the amount of insulin production was assessed using the insulin detection ELISA kit. The results indicated that like GLP-1, GLP-1R agonist M_4 can stimulate insulin secretion in the presence of 16.7 mM of glucose, and the insulin production by GLP-1 and M_4 was more than twofold compared to vehicle control at both time points (Figure ). In addition, no significant difference was observed between the amount of insulin produced by GLP-1 and M_4. These data indicated that M_4 can induce glucose-dependent insulin production in GLP-1R expressed cells and might act as a real GLP-1R agonist with potential clinical application.

Figure 6

Insulin production induced by GLP-1 and M_4 in INS-1 832/13 cells. INS-1 832/13 cells were treated with GLP-1 (181 nM) and M_4 (20 μM) in the presence of 16.7 mM glucose after 2 h of starvation with KRB buffer. Data are the average of three independent experiments, and error bars for each concentration were plotted as SEM (n = 3). Statistical analysis was done using two-way ANOVA (**p < 0.01). The comparison was done with respect to the amount of insulin produced in the corresponding vehicle control (DMSO 0.125%).

Homology Model of GLP-1R Showed Characteristics of an Active Conformation

When this project was first started several years ago, there was no experimentally determined structure available for the TM domain of GLP-1R. In 2017, a cryo-EM structure of the rabbit GLP-1R in its active conformation (PDB ID: 5VAI) and a crystal structure of human GLP-1R in its inactive conformation (PDB ID: 5VEX) were reported.[24,25] Comparison of both the backbone atoms and the predicted allosteric binding site among these two structures and the homology model generated above suggested that the homology model used for in silico screening is different from both the inactive and the active structure reported (Table ). Consistently, docking compounds M_4 and M_9 into the same predicted allosteric site of the reported active structure (PDB ID: 5VAI) respectively showed worse binding for either compound with the best docking score of about −4.8 kcal/mol.

Table 2

Structural Comparison Among the Inactive Crystal Structure (PDB ID: 5VEX), the Active Cryo-EM Structure (PDB ID: 5VAI), and the Homology Model of GLP-1R

	TM backbone atom rmsd (Å)			binding site atom rmsd (Å)
criteria	5VEX	5VAI	model	5VEX	5VAI	model
5VEX		3.58	2.45		3.67	2.59
5VAI	3.58		3.99	3.67		3.29
model	2.45	3.99		2.59	3.29

Interestingly, although the homology model is more similar to the inactive structure overall, it showed certain characteristics of an active conformation. The root-mean-square deviation (rmsd) of the predicted binding site between the homology model and the reported active structure is less than that between the active structure and the inactive structure. Further, a highly conserved interhelix interaction network (H1802.50 E2473.50 T3536.42 Y4027.57)[44] in class B GPCRs, that includes the conserved ionic lock between H1802.50 and E2473.50, and which is widely regarded to help maintain the inactive conformation of GLP-1R[25] (Figure A), was disrupted in both the active structure (Figure B) and the homology model (Figure C). In the active structure, the ionic lock between H1802.50 and E2473.50 was broken and T3536.42 moved far away from the rest of the residues in the network; while in the model, the same ionic lock was partially broken and T3536.42 also moved far away.

Figure 7

Comparison of an interaction network in the TM domain of the experimental structures and homology model of GLP-1R. (A) Crystal structure of an inactive GLP-1R; (B) cryo-EM structure of an active GLP-1R; (C) homology model.

Discussion

Type 2 diabetes and the underlying obesity is becoming a worldwide threat to human health. GLP-1R is an effective target for treatment of type 2 diabetes, and the development of small-molecule agonists will offer several potential benefits and help overcome problems associated with GLP-1 peptide drugs.[15,19−21] However, targeting the orthosteric site in GLP-1R for small-molecule discovery is ineffective because of the nature of the orthosteric binding site, which is large and relatively shallow.[12] An alternative and potentially promising approach is to target the allosteric sites on GLP-1R instead.

A major hurdle in following this approach is the lack of 3D structure information for the TM domain of GLP-1R until very recently.[23−25] Past small-molecule drug discovery efforts were often carried out by high-throughput screening.[6−8] In this work, we attempted to take the rational design approach by constructing a 3D model of the TM domain of GLP-1R in its active conformation through the combination of homology modeling, conformation sampling, and enrichment studies. First, the standard homology models of the TM domain of GLP-1R were constructed using Modeller;[29] second, enrichment studies were carried out to identify the better homology model; third, conformation sampling was performed using ProDy;[34] and finally, further enrichment studies were done to identify the best conformation for ligand screening. This approach is similar to the previously reported approach.[43]

The construction of the structural model of the TM domain of GLP-1R enabled us to apply the rational structure-based drug design techniques along with ligand-based design techniques to the discovery of small-molecule, allosteric agonists of GLP-1R. Employing the 3D model generated and performing in silico ligand-based and structure-based screening, we have identified nine compounds as potential agonists of GLP-1R for the experimental test. Two of them were confirmed as intrinsic agonists of GLP-1R by in vitro activity assay. Further studies showed that one of the two active compounds (M_4) exerted a synergistic effect on the activity of GLP-1 against GLP-1R in a dose-dependent manner (Figure A), suggesting that it does not compete with GLP-1 in binding to the orthosteric site on GLP-1R. This is consistent with our docking results that indicated that both bound at the allosteric site on the TM domain (Figure ). Together, these data showed that compound M_4 likely functions as a positive allosteric modulator agonist, although further experimental confirmation is needed.

These small-molecule agonists we identified are structurally and chemically different from those reported in the literature.[6−11] Hence, they represent novel allosteric agonists of GLP-1R. The two ligands (M_4 and M_9) were also found to induce agonistic activities in the absence of GLP-1R expression (Figure A,B). Given the high sequence similarity between GLP-1R and other class B GPCRs, for example, VIPR1 (Figure S2 online), this nonspecific effect is not unexpected. Because HEK293 cells are known to express the functional VIPR1 and undesired stimulation of this receptor could lead to numerous side effects, the nonspecific effect of M_4 on VIPR1 was studied. The data indicated that M_4 slightly induced nonspecific VIPR1 activity (Figure C) and when combined with high concentration of VIPR peptide agonist, M_4 also demonstrated a synergistic effect (Figure B), though the amplitude of such an effect was much lower than when M_4 was used with a broad range of GLP-1 (Figure A). Overall, despite the limited activity of VIPR1 induced by M_4, given the high sequence identity of the proposed binding site between GLP-1R and some other Class B GPCRs (Figure S2 online), work is in place to chemically modify compound M_4 in order to increase its potency and improve its specificity to GLP-1R. On the other hand, dual and triple agonists targeting GLP-1R along with other class B GPCRs have been demonstrated with increased therapeutic benefit recently.[45,46] It could also be of interest to explore the potential of compound M_4 as a dual or triple agonist.

The experimental structures of the TM domain of GLP-1R were finally reported in the year of 2017.[23−25] Comparison of the experimental structures with the homology model generated in this work indicated that the homology model is quite different from the reported active structure and is more similar to the inactive structure overall (Table ). On the other hand, using the homology model for in silico screening, small-molecule agonists were indeed identified and one of them can induce insulin secretion in INS-1 832/13 cells.

To some extent, this apparent discrepancy is understandable. The backbone atom rmsd between the reported inactive and active structures is ∼3.6 Å, suggesting quite a significant difference between them. On the other hand, the homology model was constructed based on another inactive structure. Although conformational sampling was carried out using the ProDy software package, it was based on the simple anisotropic network model (ANM) and only 20 ANM models were generated. With such limited sampling, only a small fraction of the entire conformational space could be explored. Given the significant difference between the inactive and active structures, the resulted model would likely remain more similar to the inactive structure of GLP-1R. This is indeed the case observed here. The rmsd of the predicted binding site between the homology model and the reported active structure is ∼3.3 Å, greater than that between the homology model and the inactive structure.

On the other hand, the homology model has shown certain characteristics of an active conformation such as partial disruption of the conserved interaction network (H1802.50 E2473.50 T3536.42 Y4027.57).[44] Consistently, the rmsd of the predicted binding site between the homology model and the reported active structure is less than that between the active structure and the inactive structure. Considering the fact that a GLP-1R can adopt multiple conformations,[6] the homology model reported here, though different from the reported active structure overall, could be much similar to another active conformation, and this may help explain its success in identifying agonists reported in this work.

The development and characterization of small-molecule allosteric agonists of GLP-1R using the rational structural model-based approach as proposed here will also help with the small-molecule drug discovery of other members of the pharmaceutical important class B family of GPCRs, in particular those whose structure is not available yet. All of the natural ligands for the class B GPCRs are moderately long linear peptide hormones whose binding site on the receptors is similar to that of GLP-1R. Because of the nature of these binding sites, small-molecule discovery targeting them remains a general challenge.[47] The existence of the allosteric site for class B GPCRs was confirmed by the reported crystal structure of CRF-1R, which showed an antagonist binding at an allosteric site deep inside the TM domain of the receptor.[22] By targeting the same allosteric site on GLP-1R using our modeling approach and structure-based drug design techniques for the development of small-molecule agonists, this work demonstrates the feasibility for small-molecule drug discovery targeting other class B GPCRs through allosteric regulations.

Methods

The proposed method includes several steps (Figure ): (i) structural modeling of GLP-1R; (ii) in silico ligand-based and structure-based ligand screening; and (iii) experimental validation of the agonistic effects of the identified small-molecule compounds.

Figure 8

Flowchart of the homology modeling and molecule screening procedure.

Structural Modeling of GLP-1R

The 3D structure model of the TM domain of GLP-1R was generated in several steps (Figure ). In 2013, two crystal structures of class B GPCRs were reported, one for the CRF-1R (PDB ID: 4K5Y)[22] and the other for the glucagon receptor (GCGR) (PDB ID: 4L6R).[27] GLP-1R has the sequence identity of ∼33% with CRF-1R and ∼54% with GCGR. Hence, both were suitable as templates for the homology modeling of GLP-1R and both were adopted for the initial homology modeling. To construct the homology model based on either template, the GLP-1R sequence and the template structure were imported in MOE (Molecular Computing Group Inc., version 2011.10), and all hetero atoms, lysozyme, and water molecules were deleted. Structure-based alignment was then carried out in MOE. The alignment results were confirmed with the alignments provided by the GPCRdb database.[28] The homology model based on either template was generated using Modeller (version 9.14)[29] based on the MOE alignment. Models of the GLP-1R based on either template were then subjected to the primary enrichment studies.

To carry out primary enrichment studies, 23 active compounds identified as GLP-1R agonists were downloaded from the PubChem BioAssay database[30] (PubChem AID: 624172). For each of the 23 compounds, 36 structurally similar, inactive compounds (decoys) were also identified using the DecoyFinder tool (http://urvnutrigenomica-ctns.github.io/DecoyFinder/) with druglike settings.

These active and decoy compounds were then prepared for docking using Ligand Preparation wizard of the Schrodinger Suite (version 2014-2).[31] Homology models generated above were prepared in Protein Preparation wizard with default settings. For the GLP-1R model based on the CRF-1R structure, the compound binding site was predicted using the Sitemap tool[32] and further confirmed based on the co-crystalized antagonist binding site on the CRF-1R structure, while for the model based on the GCGR structure, the binding site was predicted using the Sitemap tool.[32] Docking calculations were carried out using the Glide SP settings.[33] Knime (version 2.9.2) was used to automate ensemble docking calculations. The enrichment calculator script provided by Schrodinger was used to calculate the enrichment factor and to generate the ROC plots. Results from preliminary enrichment studies was used to identify the better model for further conformational sampling.

Conformational sampling was carried out using the ProDy library.[34] Twenty ANM models based on the coarse-grained C-alpha atom were first generated starting with the better homology model identified above, and then extended to the full-atom model. From these 20 models, 100 conformations were generated and optimized using NAMD.[35] From these 100 conformations, 10 conformations differing by an rmsd of at least 1.5 Å from the average rmsd model were selected for further enrichment studies using the same enrichments protocol described above. For each of the 10 GLP-1R models, ROC plots were again generated. The conformation with the best enrichment results was selected for the virtual screening process.

In Silico Screening by Ligand-Based and Structure-Based Approaches

For the initial screening, the 23 active compounds from PubChem[30] were used as the query individually to identify similar compounds in the small-molecule library of ZINC database[26] through the Shape Signatures approach.[36] Using this approach, a chemical library can be rapidly scanned for likely matches to a compound. The ZINC database of commercially available druglike molecules is currently in excess of 21 × 106 diverse molecules. Top-ranked compounds that had similar ligand properties (molecular weight, x log P, hydrogen donors, hydrogen acceptor and polar surface area) with any of the 23 active compounds were identified.

Next, those similar compounds were prepared for docking using the Ligand Preparation module of the Schrodinger Suite (version 2014-2), and the best conformation of the GLP-1R identified in the previous section was prepared in the Protein Preparation wizard with default settings. Those compounds were then docked into the same predicted allosteric site on the GLP-1R model. The Glide SP protocol with default settings was used for the docking experiments.[33] Each dock pose was then ranked based on its glide and docking scores, and top ranked poses were visually inspected. After analyzing the docking results, top-ranked compounds were identified for the experimental evaluation.

In Vitro Testing of Potential GLP-1R Agonists by CRE Luciferase Reporter Assay

Materials

HEK293 cells stably expressing the CRE/CREB luciferase reporter gene (BPS Bioscience #60515), RPMI medium (Corning #10-040), KRB (Amsbio #KRB-1000), l-glutamine (Gibco #25030-081), 4-(2-hydroxyethyl)-1-piperazineethanesulfonic acid (HEPES; Gibco #15630-080), sodium pyruvate (Gibco #11360-070), β-mercaptoethanol (MP #806444), d-glucose (#G-7528), fetal bovine serum (Fisher #03600511), penicillin/streptomycin (Corning #30-002-Cl), hygromycin B (Alfa Aesar #J60681), lipofectamine 2000 (Invitrogen #11668027), VIPR peptide agonist (Sigma #V3628), VIPR peptide antagonist (Sigma #SCP0260), GLP-1R peptide agonist (Sigma #9416), 6-well cell culture plates (Ultra Cruz #sc-204443), 96-well cell culture plates (Sigma #CLS9102), luciferase cell culture lysis reagent (Promega #E1531), luciferase assay reagent (Promega #E1501), and Ultra-Sensitive Rat Insulin Kit (Crystal Chem #90060) were purchased from vendors. Rat GLP-1R plasmid (#14944)[51] and VIPR1 plasmid (#51865)[52] were purchased from Addgene, Flag-tagged Human GLP-1R and Flag-tagged pCMV-N-Flag negative control vector were purchased from vendors (Sino Biological Inc. #HG13944-NF and #CV061), and pcDNA3.1 vector and INS-1 832/13 cells were kindly provided by Dr. Xianxin Hua (University of Pennsylvania).

Transfection and Cell Culture

HEK293 cells stably expressing CRE/CREB Reporter (luciferase) were cultured in RPMI medium supplemented with 8% (v/v) fetal bovine serum, 2% (v/v) penicillin/streptomycin, and 100 μg/mL of hygromycin B. Cells were maintained in an incubator at 37 °C with 5% CO2. Cells (1.8 × 106 per well) were seeded into 6-well cell culture plates 1 day before transfection. After overnight incubation, one well of cells was transiently transfected with 3.4 μg of rat GLP1R or empty vector using lipofectamine 2000. After 4 h of transfection, transfection medium was replaced by RPMI medium supplemented with 5% (v/v) fetal bovine serum and 2% (v/v) penicillin/streptomycin. After 24 h of incubation, cells were trypsinized and seeded into 96-well cell culture plates (5.5 × 104 cells per well) and maintained at 37 °C in a 5% CO2 incubator for 24 h. After 24 h of incubation, the transfected cells were treated with compounds as indicated.

Luciferase Assay

Compounds dissolved in 100% DMSO were diluted to indicated concentration in RPMI medium (0.5% DMSO included for all cell culture). After 4 h of treatment, cells were washed with cold 1× phosphate-buffered saline and harvested by cold luciferase cell culture lysis reagent. Luciferase activity was measured by a Wallac 1420 multiplate reader. Luciferase activity of HEK293 reporting cells cultured using 0.5% DMSO and full RPMI medium was used as a vehicle control. The protein concentration of each well was determined by Bradford assay.[37]

Nonspecific Effect of GLP-1R Agonist M_4 on VIPR1

HEK293 cells stably expressing CRE/CREB Reporter (luciferase) were cultured in RPMI medium supplemented with 8% (v/v) fetal bovine serum, 2% (v/v) penicillin/streptomycin, and 100 μg/mL of hygromycin B. Cells were maintained at 37 °C with 5% CO2. Cells (1.8 × 106 per well) were seeded into 6-well cell culture plates 1 day before transfection. After overnight incubation, cells were transiently transfected with 3.4 μg of VIPR1 or N-flag tagged pCMV3 vector using lipofectamine 2000. After 4 h of transfection, medium was replaced by RPMI supplemented with 5% (v/v) fetal bovine serum. After 24 h of incubation, cells were trypsinized and seeded into 96-well cell culture plates (5.5 × 104 cells per well) and maintained at 37 °C in 5% CO2 for 24 h. After 24 h of incubation, the transfected cells were treated with compounds as indicated.

Glucose Stimulated Insulin Production in INS-1 832/13 Cells

INS-1 832/13 cells were cultured in RPMI supplemented with 2 mMl-glutamine, 1 mM sodium pyruvate, 10% FBS, 10 mM HEPES, 100 units/mL penicillin, 100 μg/mL streptomycin, and 50 μM β-mercaptoethanol. Cells were maintained in an incubator at 37 °C with 5% CO2. To determine the effect of GLP-1R agonist M_4 on insulin production, INS-1 cells were seeded onto 6-well plates at the cell density of 1 × 106 per well. After 48 h of incubation, cells were washed twice with 1 mL of KRB and starved for 2 h in fresh KRB supplemented with 0.1% serum. After 2 h of starvation, the buffer was replaced with 1 mL of KRB containing 0.1% serum, 16.7 mM glucose, and 20 μM of M_4 or 181 nM of GLP-1 with 0.125% DMSO or 0.125% DMSO alone (vehicle control) and incubated at 37 °C with 5% CO2. After 10 and 20 min, the supernatant was collected, centrifuged at 1000 rpm for 5 min at 4 °C, and aliquoted and stored at −20 °C. These samples were used to determine the insulin concentration using the insulin detection kit following the manual.

Data Analysis

The concentration-dependent dose–response curve was generated using Graph Pad Prism 6.0 for Mac (GraphPad Software Inc., San Diego, CA). The curves were fitted based on the equation, Y = bottom + (top – bottom)/(1 + 10((log EC) or sigmoidal dose response (variable slope). The EC50 value was calculated from Prism.

Multiple Sequence Alignment

The sequence of the human GLP-1R was taken as the query sequence for similarity search against the SwissProt database[49] using BLAST[50] at www.ncbi.nlm.nih.gov to identify homologous proteins of high sequence identity. Top hit human sequences which have at least 40% sequence identity with GLP-1R, along with human VIPR1 and VIPR2 sequences, were extracted from the SwissProt database and multiple sequence alignment were performed using MOE (version 2018.0101) with the structure-based alignment function.

Comparison of Experimental Structures and Structural Model of GLP-1R

In 2017, a cryo-EM structure of the rabbit GLP-1R in its active conformation (PDB ID: 5VAI), which has >94% sequence identity with human GLP-1R, and a crystal structure of human GLP-1R in its inactive conformation (PDB ID: 5VEX) were reported.[24,25] These two experimental structures were compared with the homology model generated above and with each other. For comparison, all three structures were imported in MOE (version 2018.0101). Then, using the seven TM helix sequences of the homology model only as the standard, all other residues and molecules in the three structures were deleted. The remaining structures were superimposed and pairwise backbone atom rmsd and binding site all-atom rmsd were calculated. Next, attempts were made to dock compounds M_4 and M_9 into the same predicted allosteric site on the active structure of the rabbit GLP_1R using the Glide XP protocol with default settings.