G protein-coupled receptors (GPCRs) comprise the largest family of proteins in humans, which consists of more than 800 receptors. (1) They are membrane proteins that are composed of seven transmembrane helices (TM1-TM7) connected by three extracellular (ECL1-ECL3) and three intracellular (ICL1-ICL3) loops (collectively termed the seven transmembrane domain; 7TMD). (2) Class B1 (also known as Secretin-like) GPCRs possess a large structured extracellular N-terminal domain (ECD) and bind large peptides that are implicated in many physiological and pathophysiological conditions. (3) Binding of the endogenous peptide ligands is believed to be a two-step process whereby the C-terminus of the peptide interacts with the receptor ECD, followed by interaction of the N-terminus of the peptide with the receptor 7TMD. These interactions are associated with conformational changes that are transmitted through the 7TMD to the receptors’ cytoplasmic region. (4) Such structural changes are responsible for interaction of the receptors with, and activation of, G proteins that subsequently elicit multiple intracellular responses. (5) Recent cryogenic electron microscopy (cryo-EM) and crystallographic structural findings for Class B1 GPCRs have rendered them more accessible targets in pharmaceutical research, assisting the design and synthesis of many targeted peptide and nonpeptide analogues. (6) In recent years, advances in GPCR structural biology methods have enabled the characterization of several Class B1 GPCRs in complex with their peptide ligands and/or G proteins and have provided information about the structural mechanisms of ligand recognition and selectivity. (3,7)
The growth hormone-releasing hormone receptor (GHRHR), a Class B1 GPCR, is activated by growth hormone-releasing hormone (GHRH), a 44-residue neuropeptide produced in the hypothalamus that regulates the production and secretion of growth hormone (GH) by pituitary somatotropes. In turn, GH stimulates the secretion of hepatic insulin-like growth factor I (IGF-I), both of which are important for the regulation of growth and metabolism. Both GH and IGF-1 have also been found to act as potent mitogens in cancers. (8) However, the role of these hormones in cancer remains to be fully elucidated and exploited in a therapeutic context, with GH and IGF-1 targeting strategies resulting in side effects such as fluid retention, hypertension, and increased body weight. (9) Overexpression of GHRHR and its most prominent splice variant, SV1, has also been associated with several types of cancer, including breast, (10) prostate, (11) thyroid, (12) pancreatic, (13) and esophageal (14) cancers. SV1, in which the first three exons of GHRHR are replaced with a fragment of intron 3, inducing a novel in-frame start codon, is known to have ligand-independent activity and is believed to confer growth-stimulating effects in many tumor tissues. (15) GHRHR peptide antagonists (which block GHRH-induced and/or ligand-independent receptor activity) have been shown to inhibit growth of various cancer cell lines in addition to showing benefits in models of prostate, (16) breast, (17) pancreatic, (18) gastric, (19) and ovarian (7,20) cancer. However, features of peptide antagonists, such as the need to inject, poor pharmacokinetic profiles, and issues regarding preparation synthesis/storage and uniformity, compromise the application of these peptides in a clinical setting. To date, no nonpeptide (small-molecule) compounds targeting the GHRHR have been reported. Thus, the identification of small-molecule antagonistic compounds targeting the receptor would have a possible clinical application in the treatment of several conditions, including cancer.
We have previously performed computational studies to propose structural elements related to relative ECD movement and activation of the GHRHR (21) and, recently, cryo-EM structures of both the GHRHR and SV1 splice variant in complex with GHRH and the heterotrimeric Gαs (Gs) G protein have been determined. (7,22) The application of cryo-EM to determine the structure of many GPCRs has provided a wealth of information to aid in the understanding of drug-receptor interactions. Furthermore, this has enabled the construction of increasingly accurate in silico homology models, (23,24) also aided by recent advances in protein structure prediction with Alphafold, (24) and has resulted in many successful small-molecule compound identification projects. (25) Purchasable “fragment-like” or “lead-like” compound databases have enriched the available pool of small molecules for in silico virtual screening methodologies and subsequent in vitro testing. Indeed, several million organic molecules are commercially available, via services such as the ZINC (or other) databases. (26) When initial hits toward a target are identified, they are usually used as templates in hit-to-lead optimization schemes through organic chemistry synthesis. Currently, the availability of commercially ready-to-test compounds facilitates a process called “optimization by catalogue” to rapidly explore the chemical space of an active hit without the need of actual chemical synthesis in the hit optimization process. (27)
In this work, we have constructed a receptor-based pharmacophore using a model of the GHRHR 7TMD and performed virtual screening on a commercially available library of ∼9 million compounds for the identification of small molecules targeting this interface. The virtual screening of the compound library yielded forty-four compounds to take forward into in vitro signaling assays. Cyclic adenosine monophosphate (cAMP) enzyme-linked immunosorbent assays (ELISAs) were performed, and a potential inhibitor was identified. This was followed up by exploring its chemical space through compounds within the virtual screening (VS) library to identify compounds with chemical similarity. These were then examined in the cell-based model to establish their ability to modulate GHRHR activity, and three were found to inhibit receptor activity. We therefore report a series of novel compounds that could serve as lead compounds for further refinement with the goal of developing GHRHR small-molecule analogues for use to examine the role of GHRHR in various pathologies in a preclinical setting, in addition to clinical application.
pcDNA3mammalian expression vector encoding the human GHRHR (NM_000823.4) with a C-terminal FLAG epitope tag was a kind gift from Prof Kelly Mayo (Northwestern University, USA). pcDNA3.1(+) mammalian expression vector encoding the human GLP-1 receptor (GLP-1R; NM_002062.5) was a kind gift from Prof Jae Young Seong (Korea University, South Korea). Empty vector (pcDNA3.1(+)) was purchased from Invitrogen (Carlsbad, CA, USA). GHRH (GHRH(1–29); Sermorelin) and GLP-1 were purchased from EZBiolabs (Carmel, IN, USA). The peptide antagonist, JV-1–36, was purchased from Phoenix Pharmaceuticals (Burlingame, CA, USA). Forskolin (FSK) was purchased from Sigma-Aldrich (St Louis, MO, USA). Small-molecule test compounds were purchased from Molport (Riga, Latvia).
Modeler v9.17 (28) was used to build a homology model of the human GHRHR 7TMD in the inactive state (Uniprot ID Q02643). The structural model was built using the cryo-EM structure of the human glucagon receptor (GCGR) 7TMD domain (PDB code 4L6R) as a template. (29) The conserved S152 1.50b in TM1, H177 2.50b in TM2, E245 3.50b in TM3, W272 4.50b in TM4, N318 5.50b in TM5, G359 6.50b in TM6, and G393 7.50b in TM7, were used as reference points in the 7TMD sequence alignment. For GHRHR, corresponding residues are S140 1.50b, H165 2.50b, E223 3.50b, W250 4.50b, N296 5.50b, G339 6.50b, and G369 7.50b. Superscripts denote the single most conserved residue in each TM among Class B1 GPCRs, which is designated as X.50b (where X refers to the corresponding TM). (30) The overall stereochemical quality of the homology model was evaluated by the discrete optimized energy (DOPE) (31) and a thorough visual inspection and subsequently subjected to a 500-step energy minimization using the Amber 99SB-ILDN (32) force field.
The ∼9 million purchasable compounds of the Clean Drug-like subset of the ZINC database (33) were filtered according to molecular weight (230–500 Da), AlogP (0–6), hydrogen bond acceptors (≥2), and aromatic rings (≥1) using Openbabel, (34) finally leading to ∼2.6 million unique compounds. For each compound, 250 conformations were generated with a relative energy difference window of 20 kcal/mol.
The structural model of the GHRHR 7TMD was used to generate a pharmacophore model based on the chemical features of the orthosteric pocket residues, using the Discovery Studio 3.5 software. (35) Eight pharmacophore features were manually added, mapping potential ligand interactions with residues Y133 1.43b, K175 2.60b, V179 2.64b, S209 3.36b, T213 3.40b, D274 ECL2+2, I285 5.39b, Y342 6.53b, F345 6.56b, and L362 7.43b. Thirty-six different pharmacophores, containing all combinations of seven out of eight features, were screened one at a time, as described previously. (36) Excluded features were manually added to represent the van der Waals surface of the transmembrane orthosteric pocket of the GHRHR.
The compounds were fitted to the 36 different pharmacophores, and the fit value scoring function was used to rank the compounds from best to worst fitting. A secondary flexible fitting was performed on compounds that had a fit value of more than 3. The resulting unique compounds were ranked according to their score, and the 4326 compounds with a fit value of more than 3.6 were visually inspected for chemical complementarity. Forty-four of these compounds were selected and purchased for evaluation.
The Alphafold inactive-state specific model of GHRHR from GPCRdb (37) was used for the docking calculations. Auto Dock Tools were used to prepare the corresponding protein file with atom partial charges assigned, and compounds were docked into the binding site of the receptor using AutoDock Vina. The protein was held rigid during the docking process, while the ligands were allowed to be flexible. Docking simulations were performed using a grid box with dimensions of 20 × 20 × 20 Å and a search space of 10 binding modes, and the exhaustiveness parameter was set to 20.
Molecular Dynamics (MD) simulations of GHRHR in complex with the initial hit compound MK04 and one of the best optimized hit MK04 compounds, MK04–6, were undertaken. All simulations were performed using the GROMACS 2023.1 simulation package (38) in an NPT ensemble with Periodic Boundary Conditions (PBC). The simulations were performed in triplicates. The systems were embedded into a lipid bilayer composed of phosphatidylcholine (POPC) using the OPM database (39) and the CHARMM-GUI membrane builder. (40) The simulated systems consisted of 260 POPC molecules and were solvated using the TIP3P water model in a periodic size box of 100 × 100 × 134 Å 3. Sodium and chloride ions were added to neutralize the charge of each system with a concentration of NaCl of 0.15 mol/L. The FF19SB force field (41) was used for the parametrization of the proteins and the POPC lipids, and GAFF2 (42) was used to model the ligands.
Prior to MD simulations, the protein–ligand complexes were subjected to 5,000 steps of energy minimization using the steepest descent algorithm, followed by a standard six-step CHARMM-GUI equilibration protocol. Briefly, the first two equilibration steps were carried out in the NVT ensemble (500 ps for each step with a time step of 1 fs) at a constant temperature of 310 K using the v-rescale thermostat with a coupling constant of 1 ps. The next four steps were performed in the NPT ensemble (125 ps with a time step of 1 fs for the first step and 500 ps with a 2 fs time step for the other three steps) at a constant temperature of 310 K using again the v-rescale thermostat and a constant pressure of 1 bar using the c-rescale barostat with a time constant of 2 ps and compressibility of 4.5 × 10–5 bar–1 with gradually decreased restraint force constants to various components. These components include harmonic restraints on heavy atoms of the protein, planar restraints to hold the position of lipid head groups of membranes along the Z-axis, dihedral restraints to keep fatty acid chain double bonds in the cis conformation, and C2 chirality for each lipid molecule. Once equilibrated at constant pressure, 1 μs unbiased MD simulations were performed with the atomic coordinates of each system saved every 100 ps. Covalent bonds with hydrogen atoms were constrained using the LINear Constraint Solver (LINCS) algorithm allowing for a 2 fs time step. Temperature control was kept at 310 K using the Nosé–Hoover thermostat with a coupling constant of 1 ps, while a target pressure of 1 bar was maintained isotropically using the Parinello-Rahman barostat with a time constant of 5 ps and compressibility of 4.5 × 10–5 bar–1. Long-range electrostatic interactions were treated with the particle-mesh Ewald method with a maximum grid spacing of 1.2 Å. Finally, nonbonded interactions were calculated with a cutoff of 12 Å and a switching distance of 10 Å. The MD trajectories were postprocessed and analyzed using the standard GROMACS tool (gmx_rms) for the calculation of the ligands’ root-mean-square deviation (RMSD) and the protein’s backbone RMSD as a function of time using the starting frame as a reference structure. GetContacts software was used for the calculation of the frequency of the interactions between each ligand and receptor. (43)
Human embryonic kidney cells stably expressing large T antigen (HEK 293T, catalog no. CRL-3216, RRID:CVCL_0063, from ATCC, Manassas, VA, USA) were cultured in Dulbecco’s Modified Eagle Medium containing Glutamax (DMEM, Life Technologies; Carlsbad, CA, USA) supplemented with 10% (v/v) fetal calf serum (Biochrome; Berlin, Germany) in a humidified incubator at 37 °C with 5% CO 2. To aid cell attachment, assay plates were precoated with a 1:30 dilution of Matrigel growth factor-reduced basement membrane matrix (Corning; Corning, NY, USA) prior to cell seeding. Cells were transiently transfected with expression vectors using X-tremeGENE HP DNA transfection reagent (Roche;
