Close
About
FAQ
Home
Collections
Login
USC Login
Register
0
Selected
Invert selection
Deselect all
Deselect all
Click here to refresh results
Click here to refresh results
USC
/
Digital Library
/
University of Southern California Dissertations and Theses
/
Structure – dynamics – function analysis of class A GPCRs and exploration of chemical space using integrative computational approaches
(USC Thesis Other)
Structure – dynamics – function analysis of class A GPCRs and exploration of chemical space using integrative computational approaches
PDF
Download
Share
Open document
Flip pages
Contact Us
Contact Us
Copy asset link
Request this asset
Transcript (if available)
Content
Structure – Dynamics – Function Analysis of Class A
GPCRs and Exploration of Chemical Space Using
Integrative Computational Approaches
by
Nilkanth Patel
A Dissertation Presented to the
FACULTY OF THE USC GRADUATE SCHOOL
UNIVERSITY OF SOUTHERN CALIFORNIA
In Partial Fulfillment of the
Requirements for the Degree
DOCTOR OF PHILOSOPHY
(MOLECULAR BIOLOGY)
May 2020
ii
Acknowledgements
First and foremost, I would express my deepest gratitude and sincerest acknowledgment
to my advisor, Professor Vsevolod Katritch. During my research and study, he gave me
a lot of guidance, support and nurtured my academic life. His generosity of sharing
thoughts and knowledge, his patience to give supervision, his attitude to seek perfection
is, and will always be, the lighthouse to me.
I would like to thank Professor Vadim Cherezov, Professor Raymond Stevens, Professor
Valery Fokin, Professor Aiichiro Nakano and Professor Frank Alber for providing me
guidance, for shaping my research work, and serving on my Ph.D. dissertation
committees.
Also, I am very grateful to my colleagues in the lab, Dr. Saheem Zaidi, Dr. Arman
Sadybekov, Anastasiia Sadybekov, Dr. Jessica Grandner, Dr. Barbara Zarzycka and Dr.
Petr Popov for stimulating scientific discussions providing precious scientific insights, and
for their support and cooperation in the past five years. I also want to thank my colleagues
from the experimental team, Dr. Kyle McClary, Shuming Hao, Dr. Ming Lee, Dr. Benjamin
Stauch, Dr. Linda Johansson, Dr. Kate White, Dr. Rebecca Miller, Dr. Haitao Zhang, Dr.
Matthew Eddy, Dr. Martin Audet and Dr. Gye Won Han, for allowing me to be a part of
the exciting journey of understanding GPCRs. I am thankful to the Bridge administrative
team - Angela, Elizabeth, Kathlynn, Katie, and Allie, for their support throughout the years.
I wish to express my gratitude to my friends and classmates, who have helped me in my
Ph.D. life: Saket Choudhary, Rishvanth Prabakar, Jitin Singla, Parth Patel, Dr. Nan Hua,
Archith Nirmalchandar, Deepak Verma, Satish Sahoo, and many others.
iii
Finally, I want to thank my family, my parents, this became possible only with their
everlasting support.
iv
Table of Contents
Acknowledgements ....................................................................................................... ii
List of Figures.............................................................................................................. vii
List of Tables ................................................................................................................. x
1 Introduction............................................................................................................. 1
1.1 Rapid Growth of Class A GPCRs Structural Repertoire ....................................... 4
1.2 Computational Methods for Analyzing class A GPCRs ........................................ 6
1.2.1 Molecular docking ......................................................................................... 6
1.2.2 Molecular Dynamics Simulations ................................................................. 10
2 Computational Analysis of class A GPCRs and Their Chemical Space .......... 11
2.1 Case Study 1: Melatonin Receptors ................................................................... 11
2.1.1 Benchmarking Receptor Models ................................................................. 14
2.1.2 Experimental hit identification and validation ............................................... 19
2.1.3 Chemical & Conformational Diversity of Hits ............................................... 21
2.1.4 Structural basis of Subtype Selectivity of the Hits ....................................... 24
2.1.5 Functional selectivity of the hit ligands ........................................................ 25
2.1.6 Summary and future directions .................................................................... 29
2.2 Case Study 2: Opioid Receptors ........................................................................ 31
2.2.1 Activation-related changes in the DOP ........................................................ 31
2.2.2 A common denominator for opioid receptor activation................................. 34
2.2.3 Differences between peptide and small molecule recognition at the DOP .. 38
2.2.4 Structural basis for the selectivity of DOP agonists ..................................... 40
v
2.2.5 Structure-activity relationships of benzamide DOP agonists ....................... 44
2.2.6 Summary and future directions .................................................................... 51
3 Computational Analysis of class A GPCR Structures using MD
simulations............................................................................................................ 52
3.1 Case study 1: Angiotensin-II Receptor ............................................................... 53
3.1.1 Active-like conformation of 7TM bundle....................................................... 54
3.1.2 Helix VIII blocks G protein and β-arrestin binding........................................ 59
3.2 Case study 2: Cysteinyl Leukotriene Receptors ................................................. 64
3.2.1 An unusual pattern of microswitches ........................................................... 66
3.2.2 CysLT1R ligand-binding pocket ................................................................... 67
3.2.3 The lateral entrance of ligands directly from membrane .............................. 70
3.2.4 Molecular dynamics reveal the flexibility of the ligand access gate ............. 70
3.3 Summary and future directions ........................................................................... 74
4 Mapping the Link Between the Structure – Dynamics – Function
Relationship of class A GPCRs ........................................................................... 75
4.1 Development of the GAUGE Tool ....................................................................... 76
4.1.1 Structures of Class A GPCRs ...................................................................... 77
4.1.2 Feature calculation and selection ................................................................ 78
4.1.3 Feature-based clustering and reference index generation .......................... 81
4.1.4 Model training and performance .................................................................. 82
4.2 Static Structural Analysis .................................................................................... 83
4.2.1 Case study: Structural analysis of DOR ...................................................... 89
4.3 Evaluation of MD simulation trajectories ............................................................. 90
vi
4.3.1 Case study: MD analysis of β2AR active and inactive states ....................... 91
4.4 Summary and future directions ........................................................................... 94
5 Conclusion ............................................................................................................ 95
References ................................................................................................................... 97
Appendix .................................................................................................................... 115
A1. List of publications resulted from this doctoral research work ........................... 115
A2. GAUGE Tool for activation related conformation analysis ................................. 120
vii
List of Figures
Figure 1 Human GPCR family. Dendrogram showing the superfamily ............................ 2
Figure 2 Importance of GPCRs in physiological processes in humans ........................... 3
Figure 3 GPCR proteins as a target for drugs (red) or clinical candidates (green). ......... 4
Figure 4 GPCR structural biology – the field has shown rapid growth
8
........................... 5
Figure 5 Ligand guided receptor pocket optimization protocol ........................................ 8
Figure 6 Molecular modeling techniques employed along with docking. ......................... 9
Figure 7 Predicted binding modes of selected known MT ligands ................................. 15
Figure 8 Structural features in selected hit candidate compounds. ............................... 16
Figure 9 Functional characterization of selected hits..................................................... 22
Figure 10 Predicted binding poses for top six new chemotypes discovered with VLS .. 23
Figure 11 G-Protein vs. Arrestin pEC50 plots for MT1 and MT2 receptors.................... 26
Figure 12 Predicted binding poses for compounds 28 (a, b), and 37 (c, d) in MTR ...... 27
Figure 13. Activation-related changes in the DOP......................................................... 32
Figure 14 Polar network around D128
3.32
and basic amine ........................................... 33
Figure 15 Activation-related changes in the ECL3 of the DOP ...................................... 36
Figure 16 Water-mediated interactions during DOP activation. ..................................... 39
Figure 17 Docking pose of DPI-287-related DOP agonists ........................................... 42
Figure 18 Ligand interactions near the sodium ion binding site at DOP ........................ 48
viii
Figure 19 Overview of the MD simulation protocol ........................................................ 53
Figure 20 Different orientations of Helix-VIII in class A GPCRs .................................... 54
Figure 21 Active-like conformation of AT2R .................................................................. 56
Figure 22 Conserved L[M]
3.46
-I[A]
6.37
-Y[Y]
7.53
microswitch and sodium-binding pocket in
AT1R and AT2R ............................................................................................................ 58
Figure 23 Helix VIII blocks the putative G protein/β-arrestin-binding site of AT2R ........ 60
Figure 24 Molecular dynamics simulations of AT2R...................................................... 61
Figure 25 MD simulations of Helix-VIII and its interactions with the membrane ............ 62
Figure 26 Overall structure of CysLT1R and its comparison with other receptors ......... 65
Figure 27 Functional motifs of CysLT1R show unusual inactive-state features ............ 66
Figure 28 Orthosteric ligand-binding pocket in CysLT1R. ............................................. 69
Figure 29 MD simulation frames showing the lipid access between TM4-TM5 ............. 71
Figure 30 MD simulations of the TM4-TM5 gate closure ............................................... 72
Figure 31 Overview of GAUGE Tool ............................................................................. 77
Figure 32 Workflow for the prediction model generation and analysis .......................... 78
Figure 33 Comparison of conserved activation microswitches of the active state ......... 79
Figure 34 Transmembrane region features clustering ................................................... 85
Figure 35 PIF motif region features clustering ............................................................... 86
Figure 36 Distance matrices for (a) DRY motif and (b)the sodium ion binding site ....... 87
ix
Figure 37 Distance matrices for (a) NPxxY motif and (b) feature space including all
microswitches and activation related structural features ............................................... 88
Figure 38 Activation-related changes in the DOP.......................................................... 89
Figure 39 B2AR MD simulation analysis using the GAUGE analysis tool ..................... 92
Figure 40 Histogram of frames with activation scores predicted using GAUGE ............ 93
x
List of Tables
Table 1 Hit compounds from VLS with potency EC50 < 1 µM ....................................... 20
Table 2 Docking results for selected small molecule and peptide DOP agonists. ......... 44
Table 3 Performance of SVM based prediction models in GAUGE ............................... 82
Table 4 Activation state related score predictions of DOR structures using GAUGE .... 90
Table 5 MD simulation analysis using β2AR trajectories ............................................... 91
1
1 Introduction
G-protein coupled receptors are the most abundant membrane protein family in the
human genome, comprising 826 transmembrane protein receptors
1
. Based on their
amino acid sequences, they are further classified into five groups. These receptors
present at the lipidic bilayer of the cell membrane and share a common 7-transmembrane
topology. These receptors form the physiological basis of the most signal transduction
processes, and in humans, they are responsible for recognizing a variety of stimuli – as
small as a photon particle, to small molecular chemicals such as dopamine, adrenaline,
caffeine, and opioids, to large macromolecules such as peptides and antibodies
2
. Most
recognition interactions take place near the extracellular region of the 7-TM (TM) domain.
Information received from the extracellular environment through these recognitions can
propagate inside the cell by a concerted cascade of conformational changes leading to
the coupling of the TM domains with G-proteins, arrestins, or other scaffolding proteins
3
.
Thus, the information encoded in the conformational space gets through across the
membrane bilayer.
The class A members of the GPCR family are the largest and most diverse subfamily with
719 members, including 422 olfactory receptors responsible for taste and smell
sensations, and 287 non-olfactory members with a critical role in diverse physiological
functions(Figure 1, 2)
3
. This research work is focused on non-olfactory class A GPCR
members’ computational, structural, and chemical biology.
2
Figure 1 Human GPCR family. Dendrogram showing the superfamily and subfamilies with
Class A GPCR members shown in blue
3
.
The class A GPCRs are involved in a variety of physiological functions – including vision,
cardiovascular, endocrine, sleep and memory, pain and pleasure, reproduction, etc.
Thus, they are also at the center of current drug discovery and development efforts with
more than 34% (475 in total) presently USFDA approved drugs target one or more
proteins across the GPCR family, with more than 320 drug candidates in the clinical trials
3
(Figure 3)
4
. In humans, out of 286 class A GPCRs, 93 are established as clinical targets
and 51 are potential targets under clinical investigation.
Figure 2 Importance of GPCRs in physiological processes in humans
5
.
The immense importance of GPCRs in the treatment of a variety of diseases and
disorders, as well as in the improvement of the quality of human life makes them the most
intensively studied drug targets.
4
Figure 3 GPCR proteins as a target for drugs (red) or clinical candidates (green)
4
.
1.1 Rapid Growth of Class A GPCRs Structural Repertoire
In order to find better drugs and tool compounds to precisely control the functioning of a
GPCR, it becomes critical to understand the structure-function relationship of these
proteins. As GPCRs are membrane proteins, solving their 3D structure using
experimental techniques such as X-ray crystallography, NMR, or cryo-EM was
considered extremely challenging, until technological breakthroughs for crystallizing
5
these receptors in a membrane-like environment happened almost 15 years ago which
led to one of the first GPCR high-resolution crystal structures in 2007
6
.
Figure 4 GPCR structural biology – the field has shown rapid growth
8
.
Since then, the GPCR structural biology field has witnessed almost exponential growth,
with more than 50 unique class A GPCR structures published by the year 2019
7,8
. This
trend is more likely to continue with the breakthroughs in cryo-EM techniques enabling
the resolution of large complexes/assemblies, and the advancements in the field of X-ray
Free Electron Lasers (XFELs) for making the diffraction possible for small crystals even
at room temperatures
5
. Also, the progress in the field of Nuclear Magnetic Resonance
(NMR) spectroscopy for understanding GPCR dynamics has made significant strides in
6
the recent years
9
. Other biophysical techniques such as small-angle X-ray scattering
(SAXS) and electron-electron double resonance spectroscopies have also been valuable
tools in understanding the signaling complex formations in GPCRs
10,11
. Further, other
orthogonal techniques such as photo-crosslinking help understanding the interaction
maps of GPCR protein interactions and proven to help generate 3-D structural models of
complexes
12
. Collectively, the technological progress witnessed in the last 15 years has
a promise to uncover the structural landscape of GPCRs.
1.2 Computational Methods for Analyzing class A GPCRs
The computational molecular modeling field has made huge strides in different aspects
of aiding the structural and chemical biology for drug discovery. From screening ultra-
large chemical libraries for structure-based drug discovery, to performing MD simulations
of large systems (Millions of atoms) for longer simulation time, to predict the protein folds
using deep neural networks with unprecedented accuracies
13
, computational approaches
have become an essential component in current drug design and discovery efforts.
Furthermore, the computational field has witnessed even more dramatic progress in the
past decade – with the advent of highly scalable compute technologies such as Graphical
Processing Units (GPUs), and algorithmic advances and democratization of Machine
learning tools
14,15
. This introduction is aimed to provide a brief overview of computational
modeling tools and techniques applied to the GPCR proteins and related research.
1.2.1 Molecular docking
Docking is one of the first computational molecular modeling techniques developed for
finding new chemical hits in the structure-based drug discovery campaigns. Recent
7
studies have shown the impact of molecular docking coupled with ultra-large-scale
chemical library screening, which discovered high affinity drug-like small molecules
16
.
Molecular docking of a ligand into a protein involves conformational sampling of the ligand
at the ligand binding pocket or probable ligand interacting regions of the protein. Several
different algorithms are available to perform the ligand sampling and energy
minimizations. Some examples of such approaches include biased probability Monte
Carlo based sampling, and genetic algorithm based sampling, as implemented in MolSoft
ICM and GOLD docking tools, respectively
17
. Similarly, different scoring functions for
docking solutions have been developed based on ab initio or empirically derived terms
18
.
There are other popular and open-source docking tools available, such as DOCK and
Autodock Vina
19,20
. Currently, docking has been widely employed as shown in Figure 6,
along with other molecular modeling techniques for understanding drug function and
selectivity, and in drug design and discovery projects which complements docking and
improves the outcome of the process
21
. For example, ligand guided receptor optimization
methods such as LiBERO (Figure 5), have been proven to improve the VLS outcomes.
In such optimizations, the model with the best scores and discrimination power between
actives or true ligands against the decoys or non-binders is selected for further VLS
campaigns. Several publicly available databases can be used to generate actives set and
decoy set for a target protein
22
. The docking procedure requires a 3D model or structure
of target protein and a small molecule compound(s) to be docked as an input. If the
experimental structure is not available, one can perform homology modeling using a
closely related receptor structure as a template. Homology modeling can be done using
several publicly available tools such as Swiss-Model and Modeller
23
, or commercially
8
available tools such as QuickModel (in ICM) or Prime (in Schrodinger)
24
. In the case of
class A GPCR proteins, the 3-D structure of a GPCR or its model is validated by docking
the cognate co-crystallized ligand at the binding pocket. This is further expanded by
docking other related and high-affinity ligands to validate the consistency of interactions
and docking poses. Several commercial and open source tools available for docking,
such as ICM-Pro (Molsoft), Glide(Schrodinger), AutoDock, GOLD, etc
25,26
. .
Figure 5 Ligand guided receptor pocket optimization protocol
22
.
Once the ligand-binding pocket is characterized, and fundamental interactions are
identified, the model or structure can be applied for large scale VLS simulations. The
input chemical libraries can be generated using combinatorial chemistry or obtained
9
from the publicly accessible VLS library repositories such as ZINC database, or from
chemical suppliers such as Enamine, ChemBridge, or MolPort
27
.
Figure 6 Molecular modeling techniques employed along with docking
21
.
If there are several conformations of binding pocket or closely related homologous
proteins used for VLS, one may also use 4-D docking protocols available in ICM-Pro,
which performs docking into several binding pocket models without linearly increasing the
10
computational cost
24
. In some cases, waters can also play a critical role in mediating
interactions between the ligand and the receptor, and waters may be absent in the
experimental input structures because of insufficient resolutions
28
. In such cases,
modeling waters at the binding pocket can improve the performance of docking
predictions as reported in recent studies targeting class A GPCRs
29
.
1.2.2 Molecular Dynamics Simulations
Molecular dynamics simulations in the field of GPCRs have provided valuable insights
into the functioning of the receptors at the atomic level detail. MD simulation studies have
provided insights into the ligand biding at the receptor, interaction mapping of the ligand-
receptor complex, transition from active to inactive state, coupling of TM domains with G-
proteins, and β-arrestins, etc
30,31
. These simulations can also provide insights into the
receptor conformations relevant to ligand binding, yet unobserved in the experimental
structures
32
. MD simulations are also being utilized as a tool for refining cryo-EM maps
and model fitting to improve the resolution
33
. Additionally, several methods are available
for performing enhanced conformational sampling in MD simulations. Some examples
include application of biased potential through reaction coordinates, biasing the
simulation parameters such as temperature (Replica Exchange MD), or by selecting
conformational representatives by clustering and starting new MD simulations from these
conformations (Markov models based MD simulations)
34,35
. These techniques proved
useful and provided insights into the dynamic processes happening at longer timescales,
which may not be feasible to perform either in an unbiased MD setup or are
computationally too expensive. These computational approaches are proven to be
indispensable in understanding the GPCR structure-function-dynamics relationship.
11
2 Computational Analysis of class A GPCRs and Their
Chemical Space
Understanding the molecular recognition by a receptor provides precious insights into the
receptor function, and forms the rational basis of functional modulation by chemical
means. GPCRs interact and recognize a wide array of small-molecules. 3-D structural
information of these receptors is precious in rational design of better drugs. New
structure-based drug design studies have found promising small molecules to tackle
several outstanding challenges, for example, the opioid epidemic
36
. Recent expansions
in chemical search space with increasing computational power and available structural
information provides a unique opportunity to design and discover better therapeutics
16
.
In this chapter, two such case studies are discussed of analyzing structures of class A
GPCRs, their interactions with ligands, and virtual ligand screening for finding novel
ligands for melatonin receptors type 1 and 2 and the delta-opioid receptor – analyzing the
3-D structures of these receptors, their interactions with known ligands, and VLS
campaigns to find new ligands
1
.
2.1 Case Study 1: Melatonin Receptors
The type 1A and 1B melatonin receptors (MT1 and MT2) are G protein-coupled receptors
(GPCRs) that respond to the neurohormone melatonin (N-acetyl-5-
1
During the course of my research work, I have applied similar approaches of structural analysis
and ultra-large scale VLS campaigns for targets including Melatonin, CysLT receptor, Opioid
receptors, Angiotensin-II receptors, APJ receptor, BLT receptors, PTH1 receptor, FZD4 receptor
- using the cloud computing resources. A case study report is available on the webpage :
https://edu.google.com/why-google/case-studies/university-southern-
california/?modal_active=none
12
methoxytryptamine)
37,38
. Melatonin is found in all mammals, including humans, where it
regulates sleep and helps to synchronize the circadian rhythm with natural light-dark
cycles
39,40
. Chemically, melatonin is synthesized from serotonin in the pineal gland of the
brain during darkness
41
. Both MT1 and MT2 share canonical helical 7-transmembrane (7-
TM) topology
42,43
, although they are differentially expressed and implicated in diverse
biological functions and pathologies
44
. While exogenous melatonin has been commonly
used for the treatment of insomnia and jetlag, more effective and long-lasting MT agonists
such as ramelteon have been approved for primary chronic insomnia treatment because
of their low side-effect profile as compared to other sleeping aids such as
benzodiazepines
45,46
. Other MT agonists such as tasimelteon and agomelatine, are used
for non-24-hour sleep-wake disorders in blind individuals and as an atypical anti-
depressant for major depressive disorders, respectively
47,48
. Recent studies also suggest
MT receptors play an essential role in learning, memory, and neuroprotection
49
and these
studies illustrate the potential utility of partial and selective MT 2 receptor agonists as
antinociceptive drugs
50
. Moreover, MT2 single nucleotide polymorphisms (SNPs) are
implicated in type-2 diabetes
51
(T2D), emphasizing the importance of MT receptors in a
wide variety of functions relevant to human health and the quality of life
52
.
Although MT1 and MT2 receptors have distinctive in vivo functions, most of the currently
available drugs non-selectively activate both MT1 and MT2 receptors
53
. The discovery of
novel and selective MT ligands may, therefore, lead to useful tool compounds for the
dissection of the melatonin system, and to accelerate the development of alternatives to
existing drugs targeting MT receptors
54,55
.
13
Recently, the three-dimensional structures of MT1 and MT2 were determined using an X-
ray free-electron laser (XFEL), providing atomic-level details of receptor-ligand
interactions
42,43
. Although both receptors were resolved in complexes with agonists –
agomelatine, 2-phenylmelatonin, 2-iodomelatonin, and ramelteon, thermostabilizing
mutations that were necessary for crystallization rendered these receptors functionally
inactive. Therefore, the accuracy of these agonist-bound inactive structures in
reproducing the active-state conformation of the orthosteric pocket, and their utility in the
prospective discovery of new agonists requires further validation.
Here, we utilized the MT structural information to perform VLS on both MT 1 and MT2
receptors, using libraries of 8.4 million available-for-purchase lead-like and drug-like
compounds
27
. Subsequent experimental testing of 62 compounds selected from the top-
scoring molecules led to the discovery of ten new agonist chemotypes with sub-
micromolar potencies, with one of them, compound 21, displaying pM agonist potency
(EC50 = 0.36 nM). Six of these hits, including the most potent one, demonstrated
selectivity for MT2 . while five hits were MT2 partial agonists. Moreover, the two most
potent MT2 compounds, 21 and a close derivative of melatonin, 28, show reduced arrestin
signaling, thus resulting in substantial bias towards G-protein signaling (bias factors 15.5
and 33.9). Our results demonstrate that structure-based VLS approach can yield novel,
highly potent, and selective ligand chemotypes with potential utility as chemical probes
with distinct properties and potential lead for the treatment of circadian rhythm related
sleep and mood disorders.
14
2.1.1 Benchmarking Receptor Models
To evaluate the ability of the structure-based receptor models to recognize high-affinity
melatonin receptor ligands, we performed extensive docking of a subset of known ligands
of MT1 and MT2 receptors into (i) the unmodified 3D structures obtained from X-ray
crystallography (MT1_Xtal, MT2_Xtal), as well as (ii) into the receptor models where the
ligand-binding pocket was optimized by conformational modeling (MT 1_Opt, MT2_Opt).
Analysis of the docking poses for the known MT ligands in both crystal structures and
optimized MT receptor models showed favorable binding scores with docking poses
consistent with the orientation and binding modes of crystallized ligands (see Fig 7a-d).
The major interactions include aromatic stacking of the heterocyclic core with ECL2
hydrophobic residue F179/192
ECL2
(the residue numbers for MT1 and MT2 listed for
UniProt
56
ids: P48039 and P49286, respectively, followed by superscripted Ballesteros –
Weinstein numbering scheme
57
), as well as hydrogen bonding interactions with
N162/175
4.60
and Q181/194
ECL2
. The performance of each model was then evaluated as
the area under the corresponding receiver operator characteristic (ROC) curve (AUC),
benchmarking the ability of these models to detect ligands among decoys. The AUC
values for the optimized models of MT receptors showed substantial improvement over
AUC values for MT crystal structures (MT1_Opt=87 vs. MT1_Xtal=69; and MT2_Opt=82
vs. MT2_Xtal=70) (see Fig. 7e). Overall, these results validated the improved VLS
performance of the optimized models MT1 and MT2 receptors, which were then used for
large-scale prospective VLS.
15
Figure 7 Predicted binding modes of selected known MT ligands (a) 5-HEAT and (b) S-
26131 docked into MT 1_Opt model, whereas (c) IIK-7 and (d) 4P-PDOT docked into
MT2_Opt model. (e) ROC plots for MT receptor crystal structures and optimized models
2.1.1.1 Prospective VLS and candidate selection
A library of 8.4 million commercially available compounds was docked into the optimized
MT1_Opt and MT2_Opt structural models (see Methods), and for every compound,
docking scores and binding interactions were predicted. The top 5,000 scoring
compounds were selected from these VLS runs for each receptor, which was further
evaluated by redocking into both MT1 and MT2 receptors with increased computational
sampling. The initial hit list contained 700 compounds predicted to bind to both receptors.
To evaluate these top docking solutions, we created additional models of MT receptors
with thermostabilizing mutations restored in the proximity of the orthosteric site and
conformationally optimized to wild-type residues (A104
3.29
G and W251
6.48
F in MT1;
W264
6.48
F in MT2). The dock scores for selected MT ligands were better than -32.0 kJ/mol
in these receptor models with WT orthosteric pocket residues.
16
To capture chemotype diversity, we selected the top 500 compounds for each receptor
using chemical clustering in combination with docking scores. A final set of 62 compounds
(23 from MT1 VLS; 25 from MT2 VLS; 14 from both MT1 and MT2 VLS) were selected for
purchase based on a multidimensional composite criterion involving compound novelty,
chemical diversity, well-defined interaction patterns at the binding site interactions with
residues N162/175
4.60
and Q181/194
ECL2
, and interaction similarities with ligands
observed in the crystal structures (Fig. 8).
Figure 8 Structural features in selected hit candidate compounds. (a) 2-phenylmelatonin in
complex with MT1 receptor with the topology of chemical features shown as colored spheres
indicating R1(orange) = 5-methoxy, R2 (green)= alkylamido chain, and R3 (blue)=2-phenyl
substitutions, whereas (b) Venn diagram summarizing the topologically equivalent chemical
features in selected 62 hit candidate compounds from MT 1 and MT2 VLS
Most of the compounds represented new chemotypes with Tanimoto chemical distance
24
values > 0.22, separating them from known high-affinity MT ligands available in
CHEMBL24
58
. We also chose a close analog of melatonin – compound 28 (Tanimoto
17
distance=0.05), which to our knowledge, has not yet been characterized as a ligand for
MT receptors
58,59
. Compound 28 served as an additional positive control, which also
helped us to evaluate the effect of a single chemical group substitution in melatonin on
the binding and function at the MT receptors.
2.1.1.2 Receptor Model Preparation & Optimization
X-ray crystal structures of MT1
42
and MT2
43
receptors in complex with 2-phenylmelatonin
(PDB IDs
60
: 6ME3, 6ME6) were used to generate virtual screening models. Both
structures were converted from PDB coordinates to ICM objects using the object
conversion protocol implemented in ICM-Pro
24
. This process includes the addition of
hydrogens and assignments of secondary structures, the energetically favorable
protonation states to His, Asn and Gln side chains, and of formal charges to the ligand in
a complex with the receptor, followed by local minimization of polar hydrogens using
energy minimization protocols in ICM-Pro. The orthosteric ligand-binding pocket was
further optimized with energy-based global optimization protocol in ICM using Biased
Probability Monte-Carlo (BPMC), where the orthosteric ligand and amino acid side chains
within 5 Å radius were kept flexible and co-optimized
17
. To validate the models, a set of
20 known MT receptor ligands were selected from ChEMBL database
58
along with 780
MT receptor decoys selected for each MT1 and MT2 receptor from GPCR decoy database
(GDD)
61
and docked into crystal structures, and optimized ligand models of MT receptors.
Following the previously described ligand guided receptor binding pocket optimization
protocol, the Receiver Operator Characteristic curves (ROC) were plotted based on the
True Positive Rates (TPR) and False Positive Rates (FPR)
22
to evaluate the model
optimization. The AUC values were calculated as the areas under these ROC curves and
18
used as a model selection criterion for prospective VLS runs. The RMSD values of ligand
binding pocket side chain heavy atoms for MT1 and MT2 were 0.51 Å and 0.76 Å,
respectively, compared to their corresponding crystal structures.
To perform additional evaluation of screening results with the thermostabilizing mutants
in the proximity of orthosteric site restored to wild-type (WT) residues, the Phe residue at
F251/264
6.48
was mutated to Trp, followed by local minimization of side-chain
conformations using energy-based sampling and minimization protocols
17
. Similarly,
A104
3.29
G was also restored in the MT1 receptor model. These models were used to
evaluate the binding of selected hit candidates into the WT orthosteric pocket.
2.1.1.3 Screening Library
We selected a subset of commercially available (in-stock and on-demand) fragment-like
compounds from the ZINC database with physicochemical properties similar to already
reported melatonin receptor ligands
58, 62
. We considered compounds with molecular
weight ≤ 250 Da and logP values ranging from 1 to 5 based on the logP data of high-
affinity MT ligands. The initial dataset comprised of ~10 million compounds were
converted from SMILES to 3D format, and formal charges were assigned. This set was
further reduced to 8.4 million compounds after applying additional filters for net charges
(between -1 to 1) and removing compounds with highly reactive functional groups and
promiscuous PAINS chemotypes (“molPAINS” and “bad groups” in ICM-Pro v.3.8-6)
63
.
2.1.1.4 Virtual Ligand Screening
The VLS of 8.4 million compounds library for MT1 and MT2 models were performed using
the VLS protocol in ICM-Pro
64
. The receptor energy potential maps were calculated using
19
a fine potential grid (0.5 Å). Several energy terms, including van der Waals, hydrophobic,
electrostatic and hydrogen bonding interactions were considered for map calculations.
Full torsional flexibility of ligands was allowed, and their internal conformational strain was
considered while the receptor atoms were assigned rigid for docking. The docking was
performed using BPMC conformational sampling and energy minimization protocol
implemented in ICM-Pro for scoring and finding the best docking solutions at the default
effort level 1. These top compounds were further docked into corresponding MT receptor
models with an increased sampling effort value of 3. Each VLS run for the 8.4 million
compound library used 32,000 CPU core hours on 3 Linux workstations with a total of
192 CPU cores.
2.1.2 Experimental hit identification and validation
The selected 62 candidate compounds from VLS were tested experimentally for binding
and functional profiles in the MT receptors. Eleven compounds (see Table 1, Fig. 9)
demonstrated sub-micromolar potencies in Gi/o mediated signaling assays (18% hit rate).
Ten of these eleven compounds also showed binding affinities Ki < 10 µM in competition
binding assay. The melatonin derivative 28 identified by VLS was as potent as melatonin
itself in MT2 (EC50=0.04 nM) and had the same potency (EC50=0.04 nM) at MT1. The most
potent new chemotype, 21, displayed an EC50 = 0.36 nM for MT2, with a 30-fold selectivity
over MT1 receptor (MT1 EC50=12 nM). Overall, seven hits had EC50 < 100 nM for at least
one of the melatonin receptors. Similar to other low molecular weight MT ligands, most of
the hits identified to belong to a library of fragment-like compounds with molecular weights
less than 250 Da, and they have exceptionally high ligand efficiency (LE), exceeding the
~0.3 value considered as a standard for a promising lead. For example, compound 21
20
(Mol. Wt. = 224 Da) had the highest LE values of 0.83 kcal/mol per non-hydrogen atom
for MT2 and 0.69 kcal/mol per non-hydrogen atom for MT1 receptors
65,66
. The excellent
LE of these molecules suggests their possibility for further optimization of their drug-like
properties.
Table 1 Hit compounds from VLS with potency EC50 < 1 µM for at least one MT receptor.
MT 1 MT 2
Compound
pKi ± SEM
a
pEC
50
±SEM EC
50
(nM)
Emax
b
±SEM LE
c
pKi±SEM pEC
50
±SEM EC
50
(nM) Emax±SEM LE Selectivity
d
Tanimoto
e
21 6.31 ± 0.11 7.91 ± 0.05 12.0 93.8 ± 2.5 0.69 6.91 ± 0.05 9.44 ± 0.08 0.36 86.1 ± 3.2 0.83 30.6 0.50
23 5.42 ± 0.03 7.16 ± 0.09 57.5 96.9 ± 5.3 0.56 5.56 ± 0.13 7.69 ± 0.08 20.42 91.7 ± 3.0 0.60 2.7 0.22
28 7.78 ± 0.10 10.39 ± 0.04 0.04 95.3 ± 2.6 0.86 7.63 ± 0.08 10.35 ± 0.10 0.04 69.4 ± 4.0 0.85 0.7 0.05
29 5.22 ± 0.07 6.83 ± 0.06 144.5 87.5 ± 4.5 0.53 5.61 ± 0.05 7.46 ± 0.10 34.67 69.4 ± 8.0 0.58 3.3 0.43
37 5.07 ± 0.13 ND > 30000 ND ND 5.45 ± 0.10 6.85 ± 0.19 141.25 61.1 ± 9.1 0.53 >1000.0 0.57
44 4.19 ± 0.36 3.33 ± 0.36 57544.0 72.8 ± 4.7 0.33 4.95 ± 0.30 6.58 ± 0.13 263.03 88.9 ± 6.3 0.51 267.2 0.59
45 4.54 ± 0.15 5.06 ± 0.12 8709.6 90.6 ± 14.3 0.44 5.26 ± 0.19 6.37 ± 0.13 426.58 75.0 ± 7.4 0.56 16.9 0.59
47 4.58 ± 0.07 5.25 ± 0.16 2344.2 112.4 ± 5.2 0.46 5.91 ± 0.12 7.99 ± 0.10 10.23 91.7 ± 3.0 0.66 186.9 0.60
54 5.03 ± 0.06 6.06 ± 0.07 741.3 82.8 ± 4.3 0.54 5.56 ± 0.10 7.74 ± 0.10 18.20 75.0 ± 3.7 0.68 36.9 0.43
57 4.84 ± 0.03 5.72 ± 0.11 1778.3 87.5 ± 9.1 0.47 5.37 ± 0.04 6.88 ± 0.15 131.83 66.7 ± 8.3 0.57 10.3 0.53
62 4.32 ± 0.11 4.39 ± 0.42 42658.0 54.1 ± 10.0 0.36 5.49 ± 0.33 7.28 ± 0.14 52.48 58.3 ± 4.8 0.60 875.9 0.64
Melatonin 9.06 ± 0.14 11.38 ± 0.06 0.004 100.0 ± 5.6 0.93 9.27 ± 0.14 10.30 ± 0.14 0.05 100.0 ± 5.6 0.84 0.1 0.00
a
Standard error of the mean, N = 3.
b
Activation compared to melatonin.
c
Ligand efficiency (based on EC50).
d
selectivity in folds (calculated as: Antilog (log(Emax/EC50) MT2- log (Emax/EC50) MT1)). MT1 selectivity is
shown as underlined values.
e
Tanimoto distance to closest MT receptor ligands in ChEMBL database with
pAct>6. Hits with EC50<100 nM are displayed in bold, with Emax< 70% in italic.
21
Chart 1 Chemical structures of hit compounds with EC50 < 1 μM at MT receptors
2.1.3 Chemical & Conformational Diversity of Hits
Most of the hit compounds, as shown in Figure 9, are novel and display diverse
chemotypes distinct from known high-affinity MT ligands (ChEMBL, pAct > 6.0), with
Tanimoto distance exceeding 0.4 for all but two ligands (28 and 23). While the majority of
known MT agonists reported in ChEMBL have either substituted indene or naphthalene
core scaffold made of fused heterocyclic rings, only two of the eleven hits reported here
have fused heterocycles in place, whereas several others have two substituted aromatic
rings connected by a flexible chain. Most compounds have diverse substitutions at
positions topologically equivalent to the methoxy, acetylamido and C2 substitution at
melatonin. Three compounds-- 21, 29 and 62-- have a pyrimidine scaffold, whereas four
compounds-- 23, 37, 44, and 57-- have a methoxyphenyl group in place of the 5-methoxy
indole scaffold in melatonin. Compound 47 has substituted benzyl moieties, while
22
compound 45 has a fused tricyclic scaffold of pyridine, thiophene, and cyclopentane. Only
2 compounds, 28 and 54, have substituted indoles similar to melatonin.
Figure 9 Functional characterization of selected hits for agonist activity at MT1 and MT2
receptors in Gi-mediated signaling assays, with melatonin as reference. Results were
normalized to the Emax value (%) of receptor activation by melatonin.
The predicted binding poses of the selected hit compounds in their docking models of MT
receptors are shown in Fig. 10. Nine out of eleven hits have methoxy or a similar group
predicted to make hydrogen bonding interactions with N162/175
4.60
, which was found to
be a critical residue for receptor activation, though it plays a limited role in ligand affinity
or structural stability of the receptor
42
.
23
Figure 10 Predicted binding poses for top six new chemotypes discovered with VLS (a)
21, (b) 23, (c) 62, (d) 29, (e) 47 and (f) 54 into MT 2 receptor (purple). The center panel
shows a canonical 7-TM receptor structure of MT 2 receptor (blue helices; part of TM-V is
not displayed for clarity) in complex with 2-phenylmelatonin shown as green spheres (PDB
id: 6ME6).
Furthermore, seven of the hits were predicted to form hydrogen bonding interactions with
Q181/194
ECL2
similar to alkylamide tail of melatonin, whereas five hits were predicted to
occupy a significant space in the pocket flanked by TMs-II, III, and VII forming
hydrophobic interactions, especially with residues Y281/294
7.38
and Y285/298
7.43
, as
previously found in the MT receptor structures
42,43
. These hydrophobic interactions are
similar to those formed by the phenyl moiety of 2-phenylmelatonin. Both types of
hydrogen bonding and hydrophobic interactions were found to be critical for a ligand’s
24
steric fit into the MT receptor binding pocket and are the primary determinants of ligand
affinity.
2.1.4 Structural basis of Subtype Selectivity of the Hits
Six of the identified hits were found to be at least 30-fold more potent at MT2 as compared
to MT1. Among the hits reported here with novel scaffolds, compound 21 has the highest
potency for both MT2 (EC50=0.36 nM) and MT1 (EC50=12 nM). Compound 21 was
predicted to bind both the MT receptors in a similar orientation by forming hydrogen
bonding interactions with N162/175
4.60
and Q181/194
ECL2
with its methoxy anchor and
acetylamido tail, respectively. These interactions had been reported to be critical for
ligand affinity and potency at the MT receptors
42,43
.
Other compounds also possess remarkable MT2 selectivity. For example, compound 47
is 187-fold selective towards MT2 (EC50=10 nM for MT2, and 2.34 µM for MT1,
respectively). The pyrrole ring mimics the indole ring of melatonin, the amide group forms
hydrogen bonding with N162/175
4.60
and the chlorophenyl group forms hydrophobic
interactions with ECL2 and TM-II, III, and VII residues (Fig. 10). Despite the lack of polar
interactions with Q181/194
ECL2
, the compound displays sub-micromolar potency for MT2.
Similarly, compound 45 also lacks a substitution topologically equivalent to acetylamido
tail of melatonin (R2 feature) and yet has a sub-micromolar potency and 17-fold selectivity
for MT2 (EC50 = 427 nM). In contrast, compound 44 was predicted to form interactions
with Q181/194
ECL2
, but it lacks an R3 equivalent substitution, which still makes it 267-fold
selective for MT2 (EC50 = 263 nM). These findings suggest that either R2 or R3 could be
sufficient in maintaining the potency and selectivity at MT 2.
25
2.1.5 Functional selectivity of the hit ligands
All the discovered hits show activity as agonists in Gi-protein signaling assays at both
MT1 and MT2 receptors. This is remarkable, considering that most GPCR structures have
a limited capacity to distinguish agonists vs. antagonists
67,68
and prospective VLS
campaigns often result in antagonists even when an agonist-bound VLS model is
used
16,69
. This result may also reflect the unusual propensity of MT ligands to show
agonist activity
54
. At the same time, some compounds show functional profiles notably
distinct from full and balanced agonism, especially at MT2. Thus, four of the hits, 28, 29,
57, and 62 had their efficacy (Emax) reduced to less than 70% in MT2, and thus can be
considered partial agonists
70
. The identified hits were also evaluated for their β-arrestin
recruitment. In the case of the MT1 receptor, there are no significant deviations from the
overall balanced G-protein / Arrestin signaling profiles for most compounds (Fig. 11a).
One exception is compound 37, which completely lacks G-protein signaling, though it still
binds to MT1 and displays modest β-arrestin activity. In case of MT2, however, there are
several compounds that show a marked reduction in β-arrestin signaling compared to G-
protein, especially compounds 21 and 28, which show bias factors of 15.5 and 33.9,
respectively (Fig. 11b). These results suggest that MT ligands may show rather distinct
functional bias profiles in G-protein and β-arrestin signaling, as observed in many other
GPCRs
71,69
, though the biological importance of this bias in MT remains to be
investigated
72
.
26
Figure 11 G-Protein vs. Arrestin pEC50 plots for MT 1 and MT2 receptors.
To gain more insights into these variations in ligand activity at MT receptors, we analyzed
conformational differences among ligand-receptor pairs for these compounds. As the hit
compounds are fragment – like and may attain multiple energetically-favorable poses at
the orthosteric site upon docking, the specific conformational features driving partial
agonism remain unclear. However, analysis of compounds 28 and 37 with the most
pronounced bias to G-protein in MT2 and -arrestin in MT1 suggests potential
explanations for these phenomena.
27
Figure 12 Predicted binding poses for compounds 28 (a, b), and 37 (c, d) in MT 1 and MT 2
receptors, respectively.
Compound 28 is highly similar to melatonin, except for an amino substitution in place of
the terminal acetyl at the acetylamido tail group. This substitution renders compound 28
as a partial agonist at MT2 (Emax=69.4%) while maintaining full agonism at MT1
(Emax=95.3%). Docking predictions suggest that compound 28 assumes an orientation at
the binding pocket similar to 2-phenylmelatonin with subtle differences, as shown in Fig.
12. In MT2 orthosteric site, the ureido group of compound 28 forms hydrogen bonding
interactions with the side chains of polar residues Q194
ECL2
and N268
6.52
owing to its
28
additional amino nitrogen. Such interactions, however, are energetically unfavorable in
the case of MT1 with possible steric clashes
42,43
. Instead, the interactions of acetylamido
group of melatonin with Q181
ECL2
are replaced by a hydrogen bond between Y281
7.38
and
oxygen from the ureido group in compound 28. These interactions become favorable in
MT1 as the residue Y281
7.38
is rotated towards TM-VI placing it 4 Å away from T178
ECL2
.
In case of MT2, however, residues Y294
7.38
– T191
ECL2
are 3 Å apart, forming an
intermolecular hydrogen bond with Y281
7.38
oriented away from TM-VI allowing favorable
orientation of Q194
ECL2
to form a hydrogen bond with compound 28.
A similar pattern of ligand-receptor interaction is observed from the docking of the most
selective compound 37 into MT receptors. Compound 37 has a distinct and much bulkier
substitution with a 3-cyclopropyl-1,2,4-oxadiazol group (R2 feature). In MT1, this
oxadiazol group is predicted to form hydrogen bonds with Q181
ECL2
, T178
ECL2,
and
Y281
7.38
, while the cyclopropyl group is predicted to fit in the sub pocket formed by ECL2,
TM-V, and VII. These interaction pattern changes in MT 2 where the prominent hydrogen
bonding between side chains of T191
ECL2
and Y294
7.38
is formed, precluding hydrogen
bonding of these residues to compound 37. Moreover, the methoxy group (R1 feature)
of the compound, which maintain a hydrogen bond with N175
4.60
in MT2 is lost with
N162
4.60
side chain in MT1, due to a subtle shift of the compound. Indeed, the methoxy –
N162
4.60
interaction is found to be critical in receptor activation
42
, and loss of this
interaction is likely to explain the lack of activity of compound 37 in Gi-mediated signaling
at MT1. This interaction difference, however, does not seem to affect the β-arrestin
mediated signaling by compound 37 in MT1. Taken together, these results affirm the
critical role of N162/175
4.60
anchoring interactions in Gi-mediated receptor activation and
29
also suggest a distinct role of residue Y281/294
7.38
in governing ligand bias at MT
receptors.
2.1.6 Summary and future directions
The discovery of potent and selective MT ligands with novel chemotypes holds promise
for developing next-generation drugs to treat circadian rhythm and mood disorders, pain,
insomnia, type-2 diabetes, and cancer
49,52
. Here we utilized the three-dimensional
structures of melatonin receptors in complex with the agonist 2-phenylmelatonin
42,43
to
prepare optimized receptor models and used these models to perform prospective virtual
ligand screening of large fragment-like and lead-like compound libraries. This approach
resulted in the discovery of ten new chemotypes of potent agonists, both full and partial,
for MT1 and MT2. The number of sub-micromolar hits and potency for the best of them is
among the highest reported for the VLS campaigns in class A GPCRs
16
. Two of the hit
compounds, 37 and 62, are MT2-selective partial agonists, which may have a desirable
profile for targeting eliciting antinociceptive effects mediated by melatonin receptors and
may be potentially useful in developing novel analgesics for pain management with
reduced side effects
50
.
To our knowledge, this study represents one of the few cases of successful VLS to identify
agonists with novel chemotypes, sub-nanomolar potencies, and a high degree of receptor
subtype selectivity for a class A GPCR
16,73
-. This study also represents a successful
implementation of molecular modeling and structure-based virtual screening aimed at the
melatonin receptors over previous attempts and validates the impact of the availability of
high-quality structures in capturing vital ligand-receptor interactions
74,42,43
. Additionally,
this study verifies the importance of activated orthosteric site conformations of recent MT 1
30
and MT2 receptor crystal structures. Even though the receptors were thermostabilized by
9 and 8 point mutations (MT1/MT2, respectively) to aid crystallization while rendering them
inactive, the agonist-bound orthosteric pockets remain relevant for structure-based drug
discovery applications. Our benchmarking also underscored the critical role of structural
optimization in improving the outcomes of a structure-based VLS campaign. Another
critical aspect of this successful VLS is the discovery of novel chemotypes with reliable
docking poses. With our screening library assembled to be fragment-like in terms of
molecular weights, our obtained hits are diverse and amenable for chemical optimization
to improve their pharmacological profiles. Thus, it also proves the utility of fragment-like
compounds in the early stages of drug discovery.
The chemical diversity, selectivity, high potency, and agonist activities of the identified
hits serve as a valuable starting point for the development of tool compounds to explore
the biology of melatonin receptors. These novel chemotypes could provide the
foundations of developing next-generation treatments for insomnia, pain, sleep, and
mood-related disorders, T2D and cancer.
31
2.2 Case Study 2: Opioid Receptors
Selective activation of the -opioid receptor (DOP) has excellent potential for the
treatment of chronic pain, benefitting from ancillary anxiolytic and antidepressant-like
effects. Moreover, DOP agonists show reduced adverse effects as compared to -opioid
receptor (MOP) agonists, which are blamed for the current “opioid crisis”. Herein, we
performed computational modeling and docking of ligand interactions and receptor
activation for the first crystal structures of the DOP in an active-like state, in complex with
the peptide agonist H-Dmt-D-Arg-Phe-Sar-N-Me-3’,5’-(CF3)2-Bn (KGCHM07), and with
the small molecule DPI-287, a non-morphinan-derived DOP agonist. Our study identifies
key determinants for DOP activation and agonist recognition, with distinct differences
between peptide and small molecule agonists.
2.2.1 Activation-related changes in the DOP
Both agonist-bound structures showed all the major conformational characteristics of an
active-like state receptor. Unless otherwise indicated, we will use the higher resolution
BRIL-DOP-KGCHM07 structure for comparison with previously published inactive state
antagonist-bound DOP structures (PDB 4N6H and 4RWD)
75,76
, and with active state
structures of the MOP (PDB 5C1M and 6DDF)
77,78
and KOP (PDB 6B73)
79
. First of all,
the agonist-bound DOP structures display large outward movements of the intracellular
parts of helices V (4.5 Å) and VI (9.4-11.2 Å), and a 3.9 Å inward movement of helix VII
(Fig. 13a), which is a common feature of active conformational states of GPCRs
80
.
32
Figure 13. Activation-related changes in the DOP. Comparison of conserved activation
microswitches of the active-like DOP-KGCHM07 (orange) and DOP-DPI-287 (blue)
structures with the inactive DOP-naltrindole structure (yellow, PDB 4N6H). Structural
superposition of the (A) overall architecture, (B) PIF motif, (C) NPxxY and DRY motifs,
and (D) CWxP motif.
The shift of helix VII at the level of residue N314
7.49
(Fig. 13a), which leads to a collapse
of the allosteric sodium binding pocket in the active state GPCR structures
81
, is even more
pronounced in the determined DOP structures as compared to the active MOP and KOP
(Fig. 13b).
Interestingly, the activation-related helix VI outward movements in the agonist-bound
DOP at the level of residue F270
6.44
is also somewhat more pronounced than in MOP and
KOP. However, the very tips of helix VI (at position 6.28 as a reference) in the active state
MOP and KOP structures are more tilted by 4-6 Å, likely due to the bound G protein or
nanobody, respectively, pushing the helix VI tips further outwards
3
.
33
Figure 14 Polar network around D128
3.32
and basic amine positioning as a potential
hallmark for opioid receptor activation.BRIL-DOP-KGCHM07, orange; BRIL-DOP-DPI-
287, blue; naltrindole DOP antagonist structure (PDB 4N6H), yellow; DIPP-NH2 DOP
antagonist structure (PDB 4RWD), cyan; DAMGO MOP agonist structure (PDB 6DDF),
red.
The rearrangements in the transmembrane helices are accompanied by several changes
in conserved microswitches that are typical for GPCR activation
82
. This includes
concerted changes in the so-called P-I-F motif, where P225
5.50
moves inward by ~1.6 Å,
using I136
3.40
to change its side chain rotamer state, and facilitating the significant
movement of the bulky side chain of F270
6.44
(Fig. 13c) towards helix V by as much as
~3.5 Å at the C atoms. For comparison, this movement is ~2.6 Å in the active state MOP
and KOP structures. The changes in the P-I-F motif are coupled with changes in the
34
NP
7.50
xxY motif, collapsing the sodium binding pocket, with a ~3.5 Å inward shift of
N314
7.49
. Another residue of the sodium pocket and the NP
7.50
xxY motif, Y318
7.53
,
switches its side chain rotamer to a downward orientation, opening the intracellular
entrance to the sodium pocket (Fig. 13b). The overall conformation of the conserved
DR
3.50
Y motif remains largely unaltered between the active-like agonist-bound and the
inactive state DOP structures (Fig. 13b). Notably, the importance of the DRY motif in
maintaining the inactive state in most GPCRs is attributable to a strong salt bridge
between D
3.49
and R
3.50
residues. However, in all inactive-state structures of the DOP and
other opioid receptors, this salt bridge is already disrupted, displaying distances of >3.5Å.
Moreover, the activation changes in the corresponding MOP and KOP structures manifest
only in the side chain reorientation of R
3.50
that directly interacts with the G protein or
nanobody, which is lacking in the DOP structures. Furthermore, we evaluated the two
new DOP structures along with previously solved opioid receptor structures for their
activation-related conformations of microswitches using a GAUGE machine learning-
based prediction tool as described in Chapter 4. All the microswitches were predicted to
be in the active-like or fully activated state.
2.2.2 A common denominator for opioid receptor activation
The new DOP active-like state structures provide atomic-level insights into the critical
components of molecular recognition for small molecule and peptide agonists. The
majority of opioid receptor ligands share a basic, protonated nitrogen atom forming a salt-
bridge interaction to D128
3.32
, which itself is connected to a polar network with potential
35
involvement of Y308
7.43
, T101
2.56
and Q105
2.60
linking helices II, III, and VII. In inactive
state DOP structures, this polar network can extend to Y109
2.64
, involving water-mediated
interactions. However, in the active-like state structures, the Y109
2.64
side-chain shows a
large rotation towards helix I, uncoupling Y109
2.64
from the polar network, a mechanism
that appears to be essential for DOP activation. In case of the DOP-DPI-287 complex,
the distance for anchoring interactions between the protonated amine and D128
3.32
is ~3.0
Å. It is even higher (~3.5 Å) for the peptide KGCHM07, which shows multiple, but
comparably weak, direct or potentially water-mediated interactions with D128
3.32
(Fig. 14a
and c). Another residue of the polar network, Y308
7.43
, forms a direct hydrogen bond to
the primary amine of the peptide KGCHM07, while Y308
7.43
does not interact directly with
DPI-287’s protonated amine (N4). In both structures, Y308
7.43
positioning is preserved by
hydrogen bonds to D128
3.32
, and in DPI-287 by additional π-π-stacking interactions with
the unsubstituted benzyl moiety (Fig. 14b-c). At the same time, T101
2.56
helps in
maintaining the polar network in the DOP-DPI-287-complex by forming hydrogen bonds
with both Y308
7.43
and Q105
2.60
, while the T101
2.56
side-chain loses this interaction in the
peptide-bound DOP-KGCHM07-complex, which uncouples it from the polar network (Fig.
14a and c).
36
Figure 15 Activation-related changes in the ECL3 of the DOP and structural basis for DPI-
287 selectivity. Comparison of ECL3 conformations between (A) inactive (naltrindole,
yellow, PDB 4N6H and DIPP-NH2, cyan, PDB 4RWD) and (B) active DOP binding pockets
(DOP-KGCHM07, orange, and DOP-DPI-287, blue). (C) Alignment of agonist-bound
opioid receptor binding pockets.
The T101
2.56
side chain is also uncoupled from ligand interactions in the recently solved
peptide agonist MOP-DAMGO-complex (PDB 6DDF), whereas T101
2.56
maintains the
polar network in the peptide antagonist DOP-DIPP-NH2-structure (PDB 4RWD). D128
3.32
mutations to N or A virtually abolished KGCHM07 activity, while the potency of the small
molecule DPI-287 was reduced 10-fold for D128
3.32
N (EC50: 0.060 nM vs. 0.61 nM) and
30-fold for D128
3.32
A (EC50: 0.060 nM vs. 1.39 nM) (Fig. 14d). Similarly, previous studies
on opioid peptides reported that modifications of the N-terminal amine, including its
acetylation or substitution by a methyl group, abolished agonistic activities while retaining
low nanomolar affinity
39,40
. However, we were unable to determine the affinity of
KGCHM07 and DPI-287 for the D128
3.32
mutants, because no [
125
I]-deltorphin-I specific
37
binding could be observed. The N-terminal amine of KGCHM07 is buried 1.9 Å deeper
into the binding pocket when compared to the equivalent amine of the peptide antagonist
DIPP-NH2 (Fig. 14e).
Furthermore, the basic amine (N4) of the non-morphinan-based agonist DPI-287 is
shifted 1.4 Å deeper than the tertiary amine of the morphinan antagonist naltrindole (Fig.
14f). The MOP agonist BU72 shows the same 1.4 Å amine shift down in the binding
pocket when compared to the antagonist naloxone (PDB 4DKL). Next, we docked ten
peptide agonists and seven small molecule agonists into our new active-like DOP
structures. Accordingly, the docking poses revealed that all respective amines were
located deeper in the binding pocket when compared to antagonists (Fig. 14e-f).
Therefore, we argue that the polar network around D
3.32
plays an essential role in agonist-
induced activation at the DOP and proposes that the amine positioning is a hallmark for
opioid agonist activity.
2.2.2.1 Protocol for molecular modeling of water molecules in the binding pocket
To further evaluate water-mediated interactions in the binding pocket of DOP-KGCHM07
and DOP-DPI-287 structures, we used the energy-based water prediction tool Sample
Flood available in ICM-Pro v. 3.8.7a (MolSoft)
63
. Water predictions obtained from this
procedure were further evaluated for stability in the given space by performing energy-
based conformational minimization and sampling using water molecules and side chains
of amino acid residues located within 4 Å of predicted water molecules for at least 100,000
Monte Carlo steps. Water molecules showing consistent conformations were further
evaluated by comparison with electron density maps and considered further for docking
and ligand interaction analysis.
38
2.2.3 Differences between peptide and small molecule recognition at the
DOP
Besides the abovementioned prevalent salt bridge formation, another important anchor
of ligand interaction in opioid receptors is the phenol moiety that is conserved in many
peptides and small-molecule ligands. Accordingly, the pentapeptide-like agonist
KGCHM07 contains the N-terminal tyrosine derivative 2,6-dimethyl-L-tyrosine (Dmt),
which was shown to improve binding affinity and activity of peptidic ligands at opioid
receptors
41,42
. Its position is similar to the one observed for tyrosine in the active MOP-
DAMGO structure
32
but shows tighter binding into the pocket of the DOP. Three water
molecules were found in the KGCHM07 binding pocket, supported by weak electron
densities. These are involved in connecting Dmt1 to helices III, V, and VI via a polar
interaction network with K214
5.39
, H278
6.52
, and Y129
3.33
, which is also effectively
connected to D-Arg2 of KGCHM07 (Fig. 16a). Furthermore, our analysis suggests that
the positively-charged D-Arg2 can form a water-mediated salt bridge interaction to
D210
5.35
(Fig. 16a) supported by a 17-fold reduction in KGCHM07-binding to the
D210
5.35
N mutant. The lower resolution of the DPI-287 structure precluded robust
identification of structural water molecules in this case. However, the energy-based
prediction of water molecules suggested three tightly bound water molecules at residues
H278
6.52
and Y129
3.33
, linking DPI-287 to helices III, V, and VI as likewise observed in the
KGCHM07-bound structure (Fig. 16).
39
Figure 16 Water-mediated interactions during DOP activation. Water interactions for the
(a) KGCHM07; b DPI-287; c DIPP-NH2 (PDB 4RWD); and d naltrindole (PDB 4N6H)
binding pockets. Structural water molecules are shown as blue spheres, while additional
highly stable water molecules were predicted for the KGCHM07 and DPI-287 binding
pockets and are shown as red spheres.
The mutagenesis studies showed a ~3-fold and ~5-fold decrease of DPI-287 binding to
the Y129
3.33
F and the Y129
3.33
A mutant, respectively, and ~7-fold and ~38-fold decrease
in KGCHM07 binding for the same mutants. This finding confirms the involvement of
Y129
3.33
in water-mediated polar networks in both structures. Similarly, the EC50 in a
40
H278
6.52
A mutant was reduced ~10-fold for DPI-287 and ~50-fold for KGCHM07. The
Phe3 side chain of KGCHM07 extends towards helices II and III, extracellular loop (ECL)
1, ECL2, and is positioned in a partially hydrophobic pocket formed by Q105
2.60
,
W114
ECL1
, V124
3.28
, L125
3.29
and C198
ECL2
(Fig. 15a). In order to stabilize the DOP fusion
protein, a K108
2.63
D mutation was introduced, located at the extracellular entrance of the
binding pocket. Interestingly, our mutagenesis studies revealed that KGCHM07 binds to
the K108
2.63
D mutant with ~4-fold higher affinity as compared to the WT receptor. Docking
of the KGCHM07 peptide into a DOP model with the K108
2.63
residue as in the WT
receptor suggests that KGCHM07 binds in the same pose as in the crystal structure.
However, the ligand pushes the K108
2.63
side chain away from its optimal extended
conformation. This suboptimal interaction with the K108
2.63
side chain may explain why
KGCHM07 prefers the D108
2.63
mutant. The flexible Sar4 residue of KGCHM07 adopts
an energetically less favorable cis-amide bond to Phe3, while all remaining amide bonds
are found to be in the trans-conformation. This enables the C-terminal
bistrifluoromethylated benzyl moiety to address the ECL3 and the extracellular ends of
helices VI and VII. A large side chain-rotation of W284
6.58
by approximately 125°,
compared to other DOP structures (Fig. 15a-b), opens a hydrophobic pocket consisting
of I277
6.51
, F280
6.54
, V281
6.55
, W284
6.58
, I289
ECL3
, R291
ECL3,
and L300
7.35
, harboring the
benzyl moiety. This moiety is further stabilized by - stacking interactions and a
hydrogen bond to W284
6.58
(Fig. 14a).
2.2.4 Structural basis for the selectivity of DOP agonists
Like in MOP and KOP, the active state conformation of the DOP reveals contractions of
the orthosteric binding pocket around the agonist. Helix VI moves into the agonist-binding
41
pocket by 1.6 Å, while helix VII undergoes a 2.5 Å sideways movement. These helix
movements close to the binding pocket result in conformational changes in the ECL3 as
compared to the antagonist binding pocket. In the inactive state, R291
ECL3
stabilizes the
ECL3 by forming hydrogen bonds with the carbonyl functions of V287
6.61
and W284
6.58
(Fig. 15a). In the KGCHM07-bound structure, a large movement (10.0 Å based on the
guanidine carbon) of R291
ECL3
into the binding pocket can be observed, resulting in the
disruption of this hydrogen bond network (Fig. 15b). The side chain of R291
ECL3
is,
therefore, more flexible in the agonist-bound DOP, and its electron densities only allowed
us to model the full R291
ECL3
residue in chain A of the BRIL-DOP-KGCHM07 structure,
where it forms a lid over the hydrophobic pocket harboring KGCHM07’s
bistrifluoromethylated benzyl moiety. Although KGCHM07 is not DOP-selective
83
, our
BRIL-DOP-KGCHM07 structure reveals that R291
ECL3
is accessible to the agonist binding
pocket, and can play a role in the selectivity of DOP-binding peptides as the MOP has a
glutamic acid and the KOP a histidine in the same position (Fig. 15d).
42
Figure 17 Docking pose of DPI-287-related DOP agonists. (A) Alignment of the docking pose of
the selected DPI-287 analogs BW373U86, SNC-80, and SNC-162 (all gray) with DPI-287 (blue).
The blue box indicates the moiety with differences between these three docked analogs. (B)
Docking pose of a DPI-287 analog with N-3,4-(methylenedioxy)benzyl substitution (green) and
lacking the phenolic hydroxy function into a DOP model derived from the DOP-DPI-287 structure
with G95
2.50
D, S131
3.35
N and D108
2.63
K reversed to WT, superimposed with DPI-287 (blue). The
surface of the derivative is shown in green, and the black arrow indicates that the ligand can
penetrate deeper into the entrance of the former sodium binding pocket. (C) Superposition of the
docking poses of DPI-130 (brown) and DPI-3290 (yellow) with DPI-287 suggests that the rotated
W284
6.58
is essential for DOP binding. (D) Chemical structures and DOP binding properties
(human opioid receptors) of (+)-BW373U86, SNC-80 and SNC-162 (E) Chemical structures and
binding properties (rat opioid receptors) of DPI-130 and DPI-3290
43
The DOP-selective small molecule, DPI-287, interacts with the extracellular ends of
helices VI and VII. This finding is in agreement with mutagenesis studies that report DOP-
selective small molecules to lose their affinity upon replacement of amino acids 291-300
of the DOP by the respective MOP residues
43
, and that ECL3 mutants alone do not affect
the binding affinity of small-molecule DOP agonists
44,45
. The N,N-diethylbenzamide forms
multiple hydrophobic contacts within a pocket consisting of V281
6.55
, F280
6.54
, W284
6.58
,
and L300
7.35
(Fig. 15b and Fig. 15c). Structural comparison with other opioid receptors
reveals that DPI-287’s N,N-diethylbenzamide cannot occupy the same receptor space in
the MOP and KOP as in the DOP, due to steric interactions in positions 6.58 (charged in
the case of MOP (K305) and KOP (E297)) and 7.35 (W320 in the MOP and Y312 in the
KOP) (Fig. 15c-d). In most cases studied, L300
7.35
forms optimal hydrophobic interactions
with DOP-selective ligands, and any larger substitution would pre-empt such interactions
due to steric clashes. On the other hand, charged side chains in helix VI to replace
W284
6.58
would also make the sub-pocket less favorable for forming hydrophobic
interactions. Similarly, the indole moiety of the DOP-selective antagonist naltrindole is
known to be responsible for DOP selectivity
23,28
. Due to ligand geometry, the N,N-
diethylbenzamide of DPI-287 is shifted by approximately 2.5 Å towards the center of the
binding pocket as compared to the indole moiety of naltrindole. The N1-nitrogen of DPI-
287’s piperazine ring is not involved in specific receptor-ligand interactions. However, the
stereochemistry of the adjacent benzylic carbon is essential for DOP selectivity in DPI-
related compounds
46
. Indeed, the R-stereoisomer in this position, but not the S-isomer,
enables the N,N-diethylbenzamide of DPI-287 to address ideally the non-conserved
region close to the extracellular ends of helices VI and VII (Fig. 15b and Fig. 15c).
44
2.2.5 Structure-activity relationships of benzamide DOP agonists
The two new structures of the DOP bound to a small molecule and a peptide agonist
provide the structural basis to evaluate the key fingerprints that determine DOP
selectivity. We performed molecular docking of several small-molecule analogs of DPI-
287 at DOP, MOP, and KOP (Table 2).
Table 2 Docking results for selected small molecule and peptide DOP agonists.
PubChem
ID
Chemical Structure
Ki human
DOP, nM
Ki human
MOP, nM
Ki human
KOP, nM
Receptor
Docking
Model
Dock
Score
119029
0.32 260 130-3,400 DOP-DPI -37.56
123924
1.2-1.7 352-1,300
3,535-
4,169
DOP-DPI -25.1
9891642
2.0 (rat) 1800 (rat) 2300 (rat) DOP-DPI -31.56
45
10196358
0.4 (rat) 1.58 (rat) 21.8 (rat)
DOP-
KGCHM07
-28.3
9826770
0.18-11.39
(rat)
0.18-7.56
(rat)
0.46 (rat)
DOP-
KGCHM07
-27.7
6604878
0.5 (rat) 3,000 3,000 DOP-DPI -27.68
9871102
0.54 (IC 50)
702 (IC 50,
rat)
ND DOP-DPI
#
-33.9
11873111
8
1.7 (rat) 0.62 (rat) ND
DOP-
KGCHM07
*
-34.12
46
52937467
10.4 (rat) 0.42 (rat) ND
DOP-
KGCHM07
*
-28.23
52937468
0.6 (rat) 0.15 (rat)
118
(guinea)
DOP-
KGCHM07
*
-24.97
56834761
130 60 >1,000,000
DOP-
KGCHM07
*
-32.15
54669894
5 14.8 >1,000,000
DOP-
KGCHM07
*
-26.85
44299404
0.5-3.98 >10,000 >10,000
DOP-
KGCHM07
*
-33.59
44374513
1.8 (rat) 58.3 (rat) ND
DOP-
KGCHM07
*
-22.15
47
10918703
0.44 (rat) 53.9 (rat) ND
DOP-
KGCHM07
*
-33
10055958
1.73
147-3,930
(rat)
>10,000
(guinea)
DOP-
KGCHM07
*
-35.32
11873093
3
0.28 (rat) 0.08 (rat) ND
DOP-
KGCHM07
*
-22.72
11873093
5
0.48 (rat) 0.28 (rat) ND
DOP-
KGCHM07
*
-29.71
Docking of the selected DPI-287 analogs (+)-BW373U86, SNC-80, and SNC-162 showed
that these ligands assume a similar orientation as that of DPI-287 with comparable
docking scores at the DOP, whereas they exhibited significantly weaker docking scores
at the MOP and KOP. Within this series of compounds, the phenolic hydroxy function of
(+)-BW373U86 was either methylated (SNC-80) or removed (SNC-162), which interferes
with their ability to form polar interactions. Previous work reported a loss of affinity by ~2-
fold and ~7-fold for the DOP, respectively, while increased DOP selectivity was observed.
However, the DOP docking poses of the respective benzyl moieties of SNC-80 and SNC-
48
162 are virtually overlapping with the phenol ring of DPI-287 in the new crystal structure
(Fig. 17a) indicating that the water-mediated phenol interactions are not as important in
DOP as in MOP. (+)-BW373U86 differs from DPI-287 only by its N4-allyl moiety, but
occupies the same position as DPI-287, while the allyl group overlaps with DPI-287’s N4-
benzyl moiety. In contrast, the N4-benzyl group of DPI-287 extends further into the
entrance of the sodium-binding pocket, completely interrupting water-mediated
interactions as observed in the 1.8 Å antagonist naltrindole structure
28
(Fig. 18).
Figure 18 Ligand interactions near the sodium ion binding site at DOP. (a) Binding pocket
of the Phe3 residue of KGCHM07 and superposition with the docking pose of KGCHM07
into a DOP model with mutation K108
2.63
D restored to WT. The KGCHM07 crystal
structure is shown in orange, and KGCHM07 docked into the DOP model with mutation
K108
2.63
D restored to WT is shown in green. (b) The N4-benzyl moiety of DPI-287
penetrates the entrance of the former sodium binding pocket. DOP-DPI-287 complex
(blue); DOP-naltrindole complex (yellow, PDB 4N6H). The binding poses of DPI-287 show
its N4-benzyl moiety close to the former sodium binding pocket. The sodium ion in the
DOP-naltrindole structure is represented with a yellow sphere, and water molecules are
shown in blue spheres.
49
Conformational changes of W274
6.48
of the CW
6.48
xP motif are essential for opening up
the required space for the benzyl moiety (Fig. 16d). Moreover, it has been shown that
substitution of the benzyl group with larger residues like the N-3,4-
(methylenedioxy)benzyl moiety can further increase DOP affinity
84
. The docking poses of
this analog reveal that it can indeed penetrate further into the entrance of the sodium-
binding pocket with only minor adjustment in the pocket-lining side chains, stabilized by
a hydrogen bond to S311
7.46
(Fig. 16b). These findings indicate that the sodium-binding
pocket can be targeted by ligand interactions in the DOP, as suggested for other
GPCRs
85
. However, the functional activity of the ligands can be affected by further
intrusion into the sodium pocket, as recently shown for the leukotriene B4 receptor in
complex with a bitopic ligand protruding deep into the sodium pocket. That ligand no
longer activates the receptor but acts as an antagonist with inverse agonistic activity
86
.
The piperazine ring of DPI-287 is represented in the energy-minimized chair conformation
with axial methyl groups. Moreover, a conformational energy assessment predicted the
axial methyl conformations as more favorable than the equatorial ones (-153.43 kJ/mol
vs. -109.81 kJ/mol). Methyl groups in the axial position can form hydrophobic contacts to
Y129
3.33
, M132
3.36
, I304
7.39
, and Y308
7.43
, thereby ideally occupying the additional binding
pocket space. Furthermore, all docked analogs with the same trans-dimethyl substitutions
showed axial methyl conformations (Fig. 14f and Fig. 15a-c). The two N-(3-fluorophenyl)-
N-methylbenzamide derivatives DPI-130 and DPI-3290 differ from DPI-287, mostly in the
bulkier benzamide moiety in the meta-position of the phenyl ring (Fig. 17e). Our docking
studies show that a rotation of the W284
6.58
side chain, as observed in the DOP-
KGCHM07-complex, can open up space for the 3-fluorophenyl moiety and stabilize it via
50
- stacking interactions (Fig. 17c). However, while DPI-130 and DPI-3290 exhibit high
DOP affinity, they were shown to be non-selective and also bind MOP with low nanomolar
affinities (Fig. 17e)
87
.
2.2.5.1 Protocol for molecular docking of ligands at the opioid receptors
Selected opioid small-molecular ligands and peptide structures were obtained from the
ChEMBL or PubChem databases
58,59
. These ligands were further docked into the DOP-
DPI-287 and DOP-KGCHM07 crystal structures using the docking platform available in
ICM-Pro v. 3.8.7a. The receptor structure complexes were pre-processed for docking to
remove non-receptor fusions. Water molecules inferred from water modeling using
Sample Flood procedure implemented in ICM-Pro, and from the electron density data,
were maintained in place while preparing receptor potential grid maps for docking. Ligand
geometry was pre-optimized, and charge assignments were made using MMFF force
field
88
. Conformational sampling and optimization in the docking process were performed
using the Biased Probability Monte Carlo (BPMC) method with a sampling thoroughness
of 50, beginning with a random seed conformation every time in at least three independent
runs
17
. The resulting docking poses were analyzed for their pose and interaction
consistencies, and selected results were further analyzed. There were no distance
restraints used to bias docking scores.
Reverting the thermostabilizing mutations to WT residues in the docking models for a
DPI-287 analog, the receptor coordinates from DOP-DPI-287 complex structure were
taken and the mutations around the ligand-binding pocket D95
2.50
, N131
3.35
and K108
2.63
were restored, and the side chains were locally optimized using energy minimization-
based refinement protocol in ICM-Pro
24
. To account for the conformational changes, we
51
further performed extensive energy minimization using BPMC protocol in ICM-Pro with at
least 1,000,000 steps of conformational sampling and energy refinements. At the
beginning of the refinement process, heavy atoms were restrained by soft harmonic
tethers to the starting conformation to avoid large structural deviations, and then the
tethers gradually removed for the final refinement by reducing their weights. The refined
structure had an RMSD value of 0.19 Å for all Cα atoms compared to the DOP-DPI-287-
complex and was further used for docking of DPI-287 analog with N-3,4-
(methylenedioxy)benzyl substitution using the docking mentioned in the above protocol.
2.2.6 Summary and future directions
In this study, we implemented computational tools such as molecular modeling and
docking to understand the agonist interactions at the DOP and the role of water molecules
in mediating these interactions. Also, the changes in the ligand-binding pocket and at the
ECL3 during the receptor activation has shed more light on the peptide recognition at the
DOP. On the other hand, DOP-selective small molecules interact with the non-conserved
extracellular ends of helices VI and VII. These results provide insights into the DOP
activation and modulation processes at the atomic level details. These findings will be
fundamental for DOP-centered, structure-based drug discovery programs in a time where
opioid addiction-related deaths are dramatically increasing, and safer opioid analgesics
are urgently needed.
52
3 Computational Analysis of class A GPCR Structures using
MD simulations
Molecular dynamics simulations act as a computational microscope to provide crucial
insights into the dynamic molecular processes which may not be readily apparent from
the static experimental structures. Furthermore, the experimental structural results
obtained under non-physiological conditions, for instance, GPCR structures with fusion
partners and thermostabilizing mutations, may also affect the conformation of the
receptor. In such cases, MD simulations could provide insights into the physiologically
relevant conformational ensemble of such proteins. The time series analysis performed
on the snapshots of simulation trajectories can also provide additional information, with
even more in-depth insights into dynamic processes which may look surprising,
biologically relevant, yet unobvious to infer from the analysis of static experimental
structures. A brief overview of the MD simulation protocol in the context of membrane
protein simulations, is provided in Figure 19.
53
Figure 19 Overview of the MD simulation protocol
In this chapter, I will discuss two case studies where I have applied biased and unbiased
MD simulations to study the receptor stability and functioning, its interactions with the
ligand, and the relevance of experimental features which are unique or identified for the
first time in the X-ray crystallographic studies for Angiotensin-II Receptor Type-2 (AT2R),
and Cysteinyl Leukotriene Receptor Type-1
32,89
.
3.1 Case study 1: Angiotensin-II Receptor
Angiotensin II receptors, AT1R and AT2R, serve as critical components of the renin-
angiotensin–aldosterone system. AT1R has a central role in the regulation of blood
pressure, but the function of AT2R is enigmatic with a variety of reported effects
90
. To
determine the mechanisms of the functional diversity and ligand selectivity between these
receptors, we performed MD simulations using the crystal structures of human AT2R
54
bound to an AT2R-selective ligand and an AT1R/AT2R dual ligand, respectively,
crystallized in an active-like conformation. Unexpectedly, in these structures, the helix VIII
was found in a non-canonical position, stabilizing the active-like state, but at the same
time preventing the recruitment of G proteins or -arrestins, in agreement with the
reported lack of signaling responses in standard cellular assays. Results from the MD
simulations provide insights into the structural basis of distinct functions of the angiotensin
receptors, especially the unique feature of Helix-8.
Figure 20 Different orientations of Helix-VIII in class A GPCRs
91
3.1.1 Active-like conformation of 7TM bundle
The overall architecture of AT2R is comprised of a 7TM bundle (helices I–VII) and an
intracellular amphipathic helix VIII (Fig. 21a). Similar to AT1R and other peptide-binding
GPCRs
92
, AT2R exhibits a β-hairpin conformation of the extracellular loop 2 (ECL2) and
two pairs of disulfide bonds linking the N terminus with ECL3 (Cys35
N-term
–Cys290
ECL3
),
55
and helix III with ECL2 (Cys117
3.25
–Cys195
ECL2
; superscripts indicate residue numbers
as per the Ballesteros–Weinstein nomenclature
35
). A closer comparison of the AT1R and
AT2R structures, however, revealed substantially different conformations (Fig. 21b–e).
Although the previously determined AT1R structures captured the receptor in a classical
inactive state, all three AT2R structures, studied in this work, display an active-like
conformation. Specifically, the intracellular end of helix VI shows an outward
displacement by approximately 11.5 Å, whereas the intracellular end of helix VII exhibits
an inward displacement by 4.9 Å, as compared to the structures of inactive AT 1R (Fig.
21e). Similar large-scale shifts of helices VI and VII that open an intracellular cleft for the
recruitment of G proteins and β-arrestins have been implicated in the activation of all class
A GPCRs
82,93
. Another significant conformational rearrangement occurs at the
extracellular end of helix V, which is shifted towards the ligand-binding pocket by around
4.8 Å, compared to AT1R (Fig. 21d). Similar shifts of the extracellular part of helix V,
although typically with a smaller 2 Å amplitude, have been observed in structures of
several activated GPCRs, such as the β2-adrenergic (β2AR) and serotonin 5-HT2B
receptors
6,94
. Such agonist-stabilized displacement of helix V is essential for triggering re-
arrangements in the P
5.50
-I
3.40
-F
6.44
motif, leading to activation of these receptors.
56
Figure 21 Active-like conformation of AT2R. (a) Overall receptor architecture. Compound
1 is shown as spheres with carbon atoms in magenta, nitrogen in blue, oxygen in red, and
sulfur in yellow. Membrane boundaries, as defined by the OPM database
(http://www.opm.phar. umich.edu), are shown as red lines. (b–e), Structural comparison
of the active-like AT2R (cyan) with the inactive AT1R (green, PDB code 4YAY), and active
G-protein-bound β2AR (yellow, PDB code 3SN6) viewed from within the membrane (b),
the extracellular side (d), and the intracellular side with helix VIII removed (e). The
magnified region in c shows details of microswitches in the conserved motifs. The
backbone of AT2R is shown in the cyan cartoon, and the distinct inactive position of AT1R
57
helix VII is shown in green. Red arrows indicate conformational changes associated with
receptor activation
The large-scale rearrangements of helices during activation are accompanied by
conformational changes in the conserved microswitches. We, therefore, compared the
DR
3.50
Y, NP
7.50
XXY, and P
5.50
-I
3.40
-F
6.44
motifs of AT2R with those in the structures of
inactive AT1R and of fully activated GPCRs, such as Gs-protein-bound β2AR and arrestin-
bound rhodopsin
6,95
. The large-scale movement of helix VI upon activation is enabled by
re-arrangements in the P
5.50
-I
3.40
-F
6.44
motif (Pro223
5.50
-Ile132
3.40
-Phe265
6.44
in AT2R), as
it was previously demonstrated in other structures of activated GPCRs (Fig. 21c). In the
DR
3.50
Y motif, the Arg142
3.50
side chain of AT2R is rotated by around 90° compared to
Arg126
3.50
in AT1R, adopting a similar conformation to other fully activated GPCRs.
Additionally, in the NP
7.50
XXY motif, the side chain of Tyr318
7.53
is shifted by 6.5 Å and
rotated by 45° from the corresponding position of Tyr302
7.53
in AT1R, following the inward
movement of helix VII, as observed in other receptors upon activation. Finally, analysis
of the conserved residue switch
96
in the G-protein-binding pocket L[M]
3.46
–I[A]
6.37
–Y[Y]
7.53
(AT2R residues are shown in brackets) confirms the rearrangement of interactions in this
residue triad consistent with an active-like state of AT2R (Fig. 22a, b). Notably, instead of
a highly conserved large hydrophobic residue in position 6.37, AT 2R has a rare alanine
residue (~3% of class A GPCRs), which markedly reduces the hydrophobic contact area
between helices III and VI in the inactive state, and probably facilitates the activation-
related changes in AT2R. Therefore, all of the major conformational features indicate that
the 7TM bundle of AT2R adopts an active-like conformation, similar to that observed in
58
the crystal structures of other fully activated class A GPCRs bound to signaling partners
or their mimics.
Figure 22 Conserved L[M]
3.46
-I[A]
6.37
-Y[Y]
7.53
microswitch and sodium-binding pocket in
AT1R and AT2R. (a) Comparison of the conserved residue triad between the AT1R
(green, PDB code 4YAY) and AT2R (cyan) structures shows a rearrangement of
interactions consistent with AT2R activation. (b), Modelling of the AT2R in a hypothetical
inactive state (cyan) based on the AT1R crystal structure template (green) shows that
replacement of a large hydrophobic residue in position 6.37, which is conserved in most
class A GPCRs, to a rare small Ala258
6.37
in AT2R markedly reduces the hydrophobic
contact in this region between helices III and VI in the inactive state. c, d, Sodium binding
pocket in AT2R (c), and AT1R (PDB code 4YAY) (d) is shown as a surface with hydrogen
bonds between Asn7.46 and Asn3.35 as orange spheres. Putative sodium ion in the AT1R
59
structure (d) is shown as a solid magenta sphere, while the same position in the AT2R
structure (c) is marked as a dotted sphere. Potential sodium-coordinating residues are
shown as sticks.
In addition to the conserved microswitches, a sodium-binding site consisting of 16 highly
conserved residues is expected to undergo large-scale conformational changes upon
activation of class A GPCRs
81
. The superposition of AT2R with AT1R revealed that the
putative sodium-binding pocket in AT2R is collapsed and rearranged, hindering sodium
ion binding (Fig. 22c, d), mainly owing to the inward shift of helix VII, which is consistent
with the structures of other activated GPCRs. Only 2 out of 16 residues in the putative
sodium pocket (Ile135
3.43
and Ser311
7.46
) are different from their counterparts in AT1R
(Leu119
3.43
and Asn295
7.46
).
3.1.2 Helix VIII blocks G protein and β-arrestin binding
In most GPCR structures, helix VIII lies parallel to the membrane pointing outside of the
7TM bundle
92
, regardless of the activation state of the receptor (Fig. 23a, b). Surprisingly,
in the AT2R structures, helix VIII adopts a very different conformation by flipping over to
interact with the intracellular ends of helices III, V, and VI (Fig. 23c). This non-canonical
conformation of helix VIII was found in all AT2R structures determined in this study,
regardless of the ligand identity and different crystal packing environments, suggesting
that it is likely to be a genuine feature of AT 2R, rather than an artefact of crystallization.
Helix VIII forms a highly complementary interface with the intracellular cavity of the 7TM
bundle and is stabilized by extensive hydrophobic interactions mediated by Phe325
8.50
,
Leu329
8.54
, Val332
8.57
and Phe333
8.58
, as well as by polar interactions between
Arg324
8.49
, Gln326
8.51,
and Lys328
8.53
, and helices III, V and VI (Fig. 23d). Docking and
60
molecular dynamics simulations suggest that helix VIII stabilizes the active-like
conformation of the 7TM bundle of AT2R, while sterically blocking the binding of G
proteins and β-arrestins (Fig. 23c), which is consistent with the lack of robust downstream
signaling by AT2R as assessed by traditional G protein and β-arrestin assays
97
.
Figure 23 Helix VIII blocks the putative G protein/β-arrestin-binding site of AT2R. (a, b),
Varied positions of helix VIII in different GPCRs shown as cartoons. AT2R is in cyan, with
helix VIII colored magenta, AT1R (PDB code 4YAY) in green, NTSR1 (PDB code 4GRV)
in red, G-protein bound β2AR (PDB code 3SN6) in yellow, and arrestin-bound rhodopsin
in blue (PDB code 4ZWJ), viewed from within the membrane (a) and from the intracellular
61
side (c), Shared interaction sites of AT2R 7TM domain (grey surface) for the helix VIII of
AT2R (magenta), the C terminus of G protein (red), and the finger loop of arrestin (yellow).
(d), Interactions between helix VIII (magenta) and helices III, V, and VI (cyan) in AT2R.
Sidechains of helix VIII not resolved in the crystal structure are shown by grey carbon
atoms
In molecular dynamics simulations, helix VIII not only remained in the range of positions
within r.m.s.d. values of less than 4 Å to the crystallographic structure for a total of 4 μs
of unbiased simulation (Fig. 24a–c), but also quickly (<200 ns) and reproducibly (n = 4)
returned to this conformation after consequent perturbations of its position (Fig. 24d, e).
Figure 24 Molecular dynamics simulations. (a–c), Conformational stability of the AT2R
structure is illustrated by representative conformations (c) from a total of 4 μs of molecular
62
dynamics simulations (8 independent 500 ns runs), clustered by r.m.s.d. Traces of
distances measured between different helices are shown for apo AT2R (a) and the AT2R–
compound 1 complex (b). Distances were calculated between the centers of mass of
residues Ser792.39-Ile832.43 for helix II, Arg1423.50-Val1463.54 for helix III, Gln2536.32-
Met2576.36 for helix VI, and Phe325-Lys328 for helix VIII. (d, e), Conformational stability
of helix VIII upon perturbations, using eight starting conformations of helix VIII (d) is
revealed by r.m.s.d. traces (e), which all converge by ~250 ns of simulations. r.m.s.d.
values are calculated for the center of mass of Cα atoms of residues Phe325-Lys328
compared to the crystal structure of AT2R. Tick marks on the y-axis show the starting
frame r.m.s.d. values. Colored lines are plotted using values averaged over a 500 ps
window.
Figure 25 MD simulations of Helix-VIII and its interactions with the membrane bilayer
Results of molecular dynamics simulations for a modified AT2R model with the backbone of helix
63
VIII aligned with helix VIII from AT1R structure (PDB code 4YAY). Conformational snapshots of
the AT2R model (a) are shown for every 100 ns (blue to red spectrum) from one of the six
independent 700 ns molecular dynamics simulation runs (simulation 5). The green cartoon shows
inactive-state conformation of CCR5 (PDB code 4MBS), helix VIII of which was found to be the
closest to the final conformations of AT2R helix VIII in molecular dynamics simulations.
Intracellular view (b) of snapshots from the same MD simulation is shown but at t=0 and t=700
ns. Traces of the distance between helices VI and II (c, top curves), calculated between the
centers of mass of Cα atoms of residues Gln253
6.32
-Met257
6.36
in helix VI and residues Ser79
2.39
-
Ile83
2.43
in helix II, show a change from 21Å (active state) to under 16Å (inactive state). Traces of
the distance between helix VIII and the membrane (c, bottom curves), calculated between the
center of mass of Cα atoms of residues Arg330-Val332 and the closest phosphate atoms of lipid
molecules, indicate a gradual shift of helix VIII towards the lipid bilayer, with the distance
decreasing from ~10Å to under 3Å.
Alternatively, when helix VIII was relocated to the position observed in the AT 1R
structures, it relaxed into a canonical membrane-bound conformation, and this motion
was accompanied by an inward shift of the intracellular tip of helix VI towards its position
in the inactive state AT1R structure in three out of six independent molecular dynamics
runs (Fig. 25a–c)
64
3.2 Case study 2: Cysteinyl Leukotriene Receptors
G protein-coupled cysteinyl leukotriene receptor CysLT1R mediates inflammatory
processes and plays a significant role in numerous disorders, including asthma, allergic
rhinitis, cardiovascular disease, and cancer
98
. Selective CysLT1R antagonists are widely
prescribed as anti-asthmatic drugs; however, they demonstrate low efficiency in some
patients and exhibit variety of side effects. To gain a more in-depth understanding of
functional mechanisms of CysLTRs, we determined crystal structures of CysLT 1R bound
to chemically diverse antagonists, zafirlukast and pranlukast. The structures reveal
unique ligand binding modes and signaling mechanisms, including lateral ligand access
to the orthosteric pocket between transmembrane helices TM4 and TM5, an atypical
pattern of microswitches, and a distinct four-residue-coordinated sodium site. These
results provide essential insights and structural templates for rational discovery of more
effective and safe drugs.
65
Figure 26 Overall structure of CysLT1R and its comparison with other receptors. (A) Side
and (B) top views of CysLT 1R-pran (receptor, orange; ligand, blue) and CysLT1R-zafir
(receptor, green; ligand, yellow). (C) Superposition of CysLT 1R-pran with BLT1-BIIL260
[Protein Data Bank (PDB) ID 5X33, light blue], and P2Y1R-BPTU (PDB ID 4XNV, pink).
(D) Distribution of lipid receptors on the human GPCR sequence homology tree. CysLT 1R
and CysLT2R are marked as red dots, and other lipid receptors are marked as blue dots.
The percentage of the lipid receptors located on each branch is shown, with the exact
numbers given in parenthesis. The membrane boundaries are shown as dashed lines in
(A) and (C).
66
3.2.1 An unusual pattern of microswitches
Specific features that set CysLT1R apart from other class A receptors are evident in its
unique combination of functional motifs or “microswitches”
3
(Fig. 27). Thus, the highly
conserved DR
3.50
Y motif is replaced by FR
3.50
C in CysLT1R (Fig. 27B), eliminating the
common for inactive-state structures salt bridge between D
3.49
and R
3.50
and conferring
more flexibility to R
3.50
. Indeed, this residue adopts two different conformations in
CysLT1R-pran and CysLT1R-zafir structures (Fig. 27B). Restoring the motif by mutating
F120
3.49
D substantially decreases receptor expression and renders the protein non-
responsive to LTD4, while C122
3.51
Y mutation does not change the potency of LTD4 in
IP1 assays
99
.
Figure 27 Functional motifs of CysLT1R show unusual inactive-state features. (A)
Superposition of CysLT1R-pran (orange) with β2AR in inactive (PDB ID 2RH1; violet) and
active (PDB ID 3SN6; teal) conformations. The membrane boundaries are shown as
dashed lines. Loops are removed for clarity. (B to E) Zoom in on functional elements: DRY
67
motif (B), intracellular region (C), P-I-F motif (D), and NPxxY motif (E). A different
conformation of R121
3.50
in CysLT1R-zafir (chain A) is shown as green sticks in (B).
A conserved P
5.50
-I
3.40
-F
6.44
motif has been characterized in many class A GPCRs as a
key microswitch, coupling conformational changes in the orthosteric ligand pocket with
transitions in the cytoplasmic G protein or β-arrestin binding site
100
. While the P-I-F
residues are conserved in CysLT1R, they are found in a different conformation than in
most other antagonist-bound GPCR structures, e.g., in β2AR bound to carazolol (PDB ID
2RH1). In both CysLT1R structures, concerted rearrangements of I
3.40
and F
6.44
rotamers
result in a “switched-on” conformation of this motif, reminiscent of an active-like state (Fig.
27D). This change is related to the replacement of the neighboring “toggle switch” W
6.48
by F
6.48
. Unlike W
6.48
in other GPCRs, the side chain of F
6.48
in CysLT1R adopts a
“downward” conformation, where it would clash with the conformation of F
6.44
from the P-
I-F motif found in inactive receptors. On the intracellular side (Fig. 27C), TM6 of CysLT1R
is well-aligned with the inactive conformation of this helix in other receptors. However,
TM7 in the NPxxY
7.53
motif region is shifted (~3 Å) inward compared to the inactive state
β2AR and most other GPCRs. This TM7 shift also appears to be shared for inactive δ-
branch class A GPCR structures (Fig. 27E).
3.2.2 CysLT1R ligand-binding pocket
The ligand-binding pocket in both CysLT1R-zafir and CysLT1R-pran structures stretches
from ECL2 across the receptor toward a gap between TM4 and TM5, deep in the middle
section of the 7TM bundle (Fig. 28), which is different from any previously observed
pockets in GPCR structures. Zafirlukast and pranlukast are very distinct chemically and,
despite overall similar binding positions, show substantial variation in contacts within the
68
pocket (Fig. 28G-H). The EC part of the binding pocket is lined up with polar and charged
residues that interact with the ligands. Thus, Y104
3.33
makes extensive polar interactions
with tetrazole and benzopyran moieties of pranlukast, as well as with sulfonamide of
zafirlukast, while Y249
6.51
directly interacts only with zafirlukast (Fig. 28A,D). The
sidechain of R79
2.60
anchors pranlukast via a salt bridge to its tetrazole moiety.
In contrast, zafirlukast does not directly contact R79
2.60
, and its toluene group points
towards TM6 and TM7, making hydrophobic interactions. The glycine carboxyl of LTD4
forms a salt bridge with R79
2.60
and a hydrogen bond with Y104
3.33
, mimicking the
tetrazole group of pranlukast, while the cysteinyl carboxyl engages in an additional salt
bridge and hydrogen bonds with R22
1.31
, Y26
1.35
, Y83
2.64
, and Q274
7.32
. Other key LTD4
interactions include a hydrogen bond to Y249
6.51
, stacking with R253
6.55
, and hydrophobic
interactions with V186
5.35
and V277
7.35
. Indeed, mutations in most of these residues
dramatically reduce or altogether abolish LTD4 signaling
99
.
69
Figure 28 Orthosteric ligand-binding pocket in CysLT1R. (A and D) Details of ligand-
receptor interactions for pranlukast (A) and zafirlukast (D). (B and E) Pocket shapes for
pranlukast (B) and zafirlukast (E). (C and F) Pocket entrance for pranlukast, closed “gate”
(C), and zafirlukast, open “gate” (F). (G and H) 2D representations of receptor-ligand
70
interactions for pranlukast (G) and zafirlukast (H). Water molecules are shown as red
spheres in (A). Residues engaged in the same type of interactions with both zafirlukast
and pranlukast are colored in light green, and those engaged in different types are colored
in orange (G, H). The membrane boundary is shown as a dashed line in (B, C, E, and F).
3.2.3 The lateral entrance of ligands directly from membrane
While aromatic rings of benzamide groups of zafirlukast and pranlukast overlap in the
structures, the further ligand paths towards the gap between TM4 and TM5 as well as the
receptor conformation in this region diverge dramatically (Fig. 28B, E). In CysLT1R-pran,
the extended phenyl butyl chain of the ligand is enclosed entirely within the 7TM bundle,
with only a small opening between T154
4.56
and V192
5.41
towards the lipid bilayer. In
CysLT1R-zafir, a large-scale (~7 Å) movement of the extracellular tip of TM5 and a
corresponding rearrangement in ECL2 create a wide sideways opening in the 7TM
bundle, which accommodates the indole group of the ligand and is sufficient for passing
the whole ligand into the pocket laterally (Fig. 28C,F). The tilting motion of the extracellular
part of TM5 is facilitated by the highly conserved P201
5.50
hinge, which is made even
more flexible by G197
5.46
as found in CysLTRs, LPA4,5,6, and BLT1. Such TM5 plasticity
suggests that the opening between TM4 and TM5 can serve as an important gate for the
lateral ligand entrance into the orthosteric pocket.
3.2.4 Molecular dynamics reveal the flexibility of the ligand access gate
Molecular dynamics simulations demonstrate that the open gate conformation of the
pocket as found in the CysLT1R-zafir complex is metastable – after zafirlukast is removed,
the open pocket collapses typically within 200 ns into a closed conformation, resembling
the one observed in the CysLT1R-pran structure (Fig. 28). In some cases, however, the
71
open conformation lingers much longer due to a POPC lipid molecule entering the pocket
from the membrane. Moreover, we observed that the gate could close and open
spontaneously within a 1 µs simulation, indicating a relatively low barrier for transitions
between the open and closed states, supporting the lateral ligand access route into the
orthosteric pocket in CysLT1R (Figure 29, 30).
Figure 29 MD simulation frames showing the lipid access into the gap between TM4-TM5
helices near 350 ns, and the closed state after lipid exits
72
Figure 30 MD simulations of the TM4-TM5 gate closure. Distance plots for interactions
between TM4 (Cα atom of P157) and TM5 (Cα atom of V186), as well as between three
TM4 residues (T154
4.56
, P157
4.59
, and F158
4.60
), and the closest heavy atom of TM5 in ten
independent simulation runs of the apo state based on the CysLT 1R-zafirlukast structure
over the course of a cumulative simulation time of 10 µs. The distance sampling rate was
ten frames per ns, and solid lines represent moving average values from 200 frames in all
cases. Dotted lines on the right show the minimum distances for the open state, whereas
solid lines show closed state distances for the respective residues.
73
3.2.4.1 Protocol for MD simulations of CysLT1R
The CysLT1R-pran and CysLT1R-zafir structures were pre-processed to assign their
protonation states and to model missing side chains using ICMFF energy-based
optimization protocols available in ICM-Pro
24
(v3.8-6). Missing loops were modeled using
Loop modeling and Regularization protocols in ICM-Pro. The pre-processed structures
were used for molecular dynamics simulation as previously described
89
, using input files
generated by the CHARMM-GUI web-server
101
. The initial membrane coordinates were
assigned by aligning the receptor models to CB1 receptor coordinates retrieved from the
OPM database
102
. The structure was simulated in a periodic box containing 178 POPC
lipids, 11,898 water molecules, 31 sodium, and 44 chloride ions. After the initial energy
minimization, the system was equilibrated for 10 ns, followed by production runs of up to
1 µs for CysLT1R-zafir and 650 ns for CysLT1R-pran apo state structures. The simulations
were performed using Gromacs (v.2018) simulation package
103
, and the plots were
generated using Matplotlib plotting package available in Python. The simulations were
performed either on NVIDIA P100 GPU enabled nodes made available by the Google
Cloud Platform (GCP) or with NVIDIA K-80, or P100 GPU enabled clusters at the High-
Performance Computing Center at the University of Southern California.
74
3.3 Summary and future directions
MD simulations are a powerful tool for probing molecular processes at the atomic level
details. It also provides insights into the temporal domain of the events, which could shed
light onto structural features that can be used to modulate the function of these systems.
Crucial details linking the class A GPCRs structure and function have been gained from
the MD simulations. In this chapter, MD simulations have provided some useful insights
into unique structural features observed in the X-ray crystal structures. In the case of
AT2R, we have performed the dynamics of an important structural feature - the Helix VIII,
and the structural basis of how it can block the receptor coupling to G-proteins, and as a
result, the functioning of the receptor itself. Similarly, the unique structural feature
observed in the CysLT1 receptor structure – a cleft between TM4 and TM5, and its
potential significance in the receptor interactions with lipids and lipid like ligands.
As the GPCR structural space expanding even more rapidly with cryo-EM techniques,
MD simulation would play an important role in bridging the structural-functional link, as
well as modeling the receptor-ligand interactions, which may not be delineated by
experimental methods due to resolution limitations.
75
4 Mapping the Link Between the Structure – Dynamics –
Function Relationship of class A GPCRs
The class A GPCR receptors possess a broad spectrum of conformational ensemble
characterizing conformational populations for active, inactive, and intermediate states.
Recent findings from the experimental structural studies and computational modeling
have shed more light on the conformational ensemble, and their link with the functioning
of these receptors
104
. The ability to transition within conformational landscape enables
GPCRs to recognize and interact with a variety of entities to mediate information between
the extracellular and intracellular environments. In most class A GPCRs, the orthosteric
ligand interaction site is located at the extracellular side. External signal carriers (light,
small chemicals, biomacromolecules) interact at this external site, which leads to
conformational changes in the transmembrane region of a receptor. These
conformational changes propagate to the intracellular side, where, based on the signal,
the receptor becomes primed to interact with an appropriate coupling protein, which
further propagates the signal inside the cell
11
. Assessment of conformational changes
becomes more challenging as these changes throughout the transmembrane domains
show weak coupling, culminating into significant structural movements near the
intracellular transmembrane region. These assessments become even more challenging
when performed for thousands to millions of conformations coming from MD simulations.
In order to gain better insights into the dynamic information transfer process across the
membrane bilayer, it becomes essential to map the conformational landscape and the
structural features defining it with the function of a receptor. Some recent work on
76
analyzing static structures of class A GPCRs and evaluating their dynamics have shed
some light on these processes. For example, recent studies have reported a distance-
based structural evaluation method to identify distinct key interacting residue pairs from
MD simulations
105
. In other studies, a universal activation index based on the intracellular
structural features involving transmembrane distances is proposed
106
. Despite these
recent efforts of evaluating the conformational landscape of class A GPCRs, it remains
unclear to quantitatively map the key structural features or activation switches to their
functional relevance. In this work, we present a machine learning-based structural
assessment tool – GAUGE for automated assessment of structural features and their
activation related scoring. This tool could help provide a rapid evaluation of an
experimental structure or a model of class A GPCR and could also provide an automated
and holistic assessment of the conformational ensemble obtained from computational MD
simulations or Cryo-EM studies. Furthermore, this would also provide a more detailed yet
simplified assessment of the dynamic conformational transitions relevant to distinct
functions and the underlying sequential patterns leading to such transitions.
4.1 Development of the GAUGE Tool
The experimental structural repertoire of class A GPCRs has expanded rapidly in recent
years, and currently 306 structures of class A GPCRs available in PDB. We utilized this
structural information in order to prepare a Support Vector Machine (SVM) based
prediction models. Figure 31 provides an overview of the development and application of
this tool.
77
Figure 31 Overview of GAUGE Tool
4.1.1 Structures of Class A GPCRs
We considered available 306 experimental structures of class A GPCRs, which represent
77 unique class A GPCR proteins
60
. The steps involved in the preparation of the GAUGE
tool are shown in Figure 32. In order to make these receptor structures consistent, we
removed all non-main chain molecules, including the residues from the main protein
chain, which may belong to the fusion proteins or engineered non-protein molecules. The
remaining main chain amino acid residues were renumbered after aligning with their wild-
type sequence obtained from UniProt
56
. Finally, these residues were mapped to their
structural numbers using BW and GPCRdb numbering schemes. These residue
numbering scheme for class A GPCRs has been derived from the structural positioning
78
of amino-acid residues in the 7-TM domain, which also considers the residue
conservation throughout the class A GPCR subfamily proteins.
Figure 32 Workflow for the prediction model generation and analysis
4.1.2 Feature calculation and selection
Once the receptor structures were preprocessed for consistencies and mapped with their
corresponding GPCRdb numbers
57,107
, we calculated structural features of selected
microswitches and sub-regions of transmembrane domains which are previously
characterized as key components propagating information across the membrane bilayer.
While the class A GPCRs show a significant sequence variation in the majority of
sequence, the residues forming these key structural features show high conservation.
These key features/locations include PIF, DRY, and NPxxY motifs, Na
+
binding site and
transmembrane helices near the intracellular region of membrane.
79
Figure 33 Comparison of conserved activation microswitches of the active state (orange,
PDB: 4LDL) and inactive state and (blue: PDB: 2RH1) structures of β2AR. Overall
structural architecture(center), PIF motif, DRY, and NPxxY motifs(right), TM domains near
intracellular region (left)
A variety of structural measurements were considered while calculating and selecting the
feature space using the MDAnalysis package available in Python v.3.6
108
. For the PIF
motif, the pairwise distances between Ca-Ca atoms and the center of geometry (COG) of
side-chain heavy atoms for P, I, and F were calculated, along with an additional pair of
F
6.44
– S
3.39
distance. This additional feature was given higher relative weight by the SVM
model as the location of Ca in S
3.39
remains less perturbed by the PIF movements, which
captures the sidechain movement of F. Thus, seven features were considered for
80
generating a prediction model for the PIF motif. In the case of DRY motif, we considered
three features. The first is the distance between the side chain C atoms, which are farthest
from the C atom of the corresponding residues D
3.49
and R
3.50
of the DRY motif. This
feature indicates the formation of ionic lock between these two residues, which is a key
activation related feature observed in some of the active state structures. The other
features for the DRY motif include the distance of the side chain C atom (farthest from
the corresponding C atom) to the center of geometry of Ca atoms belong to 2.50 and 5.50
locations in the GPCRdb numbering. This feature gives the side-chain orientation of R
3.50
residue with regards to the axis parallel to the transmembrane helical axis. The third
feature is obtained from a distance between the side chain C atoms (farthest from the
corresponding C atom) from the residues 3.50 and 6.34. This feature captures the local
changes near TM-3 with regards to presence or absence of activation related interactions
with the intracellular coupling partners of class A GPCRs. For the NPxxY motif, we
calculated 12 distance-based features measuring distances between the center of
geometry of Ca atoms of the seven residues located at the intracellular region of TM
domain, to the center of geometries of side-chain heavy atoms of residues at N
7.49
, P
7.50,
and Y
7.53
locations. These features capture the local structural rearrangements near the
NPxxY motif related to the activation of the receptor. The fourth set of features is
calculated for the Na
+
site which is a key feature coupling orthosteric site to the
intracellular region of the TM helices
81
. We selected four features describing the Na+
binding site, which includes distances between the COG of side-chain heavy atoms for
residues at locations S
3.39
-N
7.49
, W
6.48
- N
7.49
, N
3.35
- W
6.48
, and N
3.35
- N
7.49
. These side
chains undergo rearrangements to lose the coordination and the collapse of the Na+
81
binding pocket which is a key activation event. Finally, for the transmembrane domain
conformations near the intracellular region, we calculated three features measuring the
distances between the COG of the bottom four residues from TM2 – TM6, TM3 – TM6,
and TM3 – TM7.
4.1.3 Feature-based clustering and reference index generation
The conformational ensemble of the class A GPCR represents a diverse set of
populations with their distinct preferences towards the receptor functioning. In order to
perform a holistic analysis of activation related conformational spectrum and improve the
robustness of prediction model, we considered all the experimental 3-D structures of class
A GPCRs available in the PDB
60
. For every activation related location considered here
as described earlier, we performed a simple agglomerative hierarchical clustering using
Unweighted Pair Group Method with Arithmetic Mean (UPGMA)
109
. In the initial step of
constructing a UPGMA based clustering, a distance matrix based on the set of descriptors
is generated for an activation related region. For example, we calculated the pairwise
distances for the center of mass of the Cα atoms of selected residues at the intracellular
TM regions (i.e. TM2 – TM6, TM3 – TM6, and TM3 – TM7) for all preprocessed 235
receptor structures. At this point, we excluded several opsin structures in order to prevent
bias in the input dataset due to their higher representation, and we kept only a few opsin
structures with higher resolution and distinct activation states. Thus, obtained descriptor
distances were used for UPGMA based clustering
109
. From the distance matrix obtained
from the UPGMA clustering, we calculated a “distance vector” from an inactive state
reference receptor structure of high-resolution A2a receptor (PDB: 4EIY)
110
to other
structures. These distances give a measure of how different a conformation or local
82
geometry of a given receptor is from the reference structure (inactive state A2a receptor).
Similar feature-based distance calculations and clustering were performed for all the
microswitches and activation related key locations.
4.1.4 Model training and performance
Once the descriptor set and a reference distance index are calculated for a given
microswitch or a substructure, we generated an SVM based regression model to measure
the distance from inactive state conformation for input receptor structure. We used the
SVM implementation from Sci-Kit learn
111
(sklearn v. 0.21.3) available as Python (v. 3.6)
library for machine learning tools, along with other libraries for numerical computations,
including NumPy
112
(v. 1.17.2) and data analysis library Pandas
113
(v. 0.25.2). The input
target column values or the distances calculated as described previously, are used for
training the prediction models. These set of values are normalized between -1 to 1
corresponding to the conformations ranging from inactive states to active states. The
input dataset comprising data from 235 receptors was split into training (80%) and
test(20%). The model training process was performed with 3-fold cross-validation. Further
details regarding the regression model parameters are provided in the appendix. The
model performance matrices for different regression models are shown in Table 3.
Table 3 Performance of SVM based prediction models for conformational analysis of class
A GPCRs. (MAE: Mean absolute error; MSE: Mean squared error; R
2
: coefficient of
determination)
83
The classification models based on SVM were also generated using a dataset of 30
receptors (18 Inactive state + 12 Active state; activation related annotations obtained from
GPCRdb
8
; appendix A2) structures with 80% test and 20% training split. In all cases, the
classification models achieved absolute performance by identifying the active and inactive
training data correctly with a caveat of smaller test and training set. These generated
models for regression and classification of class A GPCR structures were used for
analyzing publicly available class A structures and MD simulation data.
4.2 Static Structural Analysis
The structural studies of class A GPCRs suggest the conformational landscape
associated with the activation mechanism has a variable degree of energy barriers
throughout the TM domain. For instance, the orthosteric ligand binding site shows a small
conformational change upon binding with an agonist or an antagonist. These subtle
changes have weak coupling with W
6.48
and the PIF motif
93
. These changes link with the
Na+ binding site rearrangements, where presence of Na+ ion increases the activation
barrier by acting as a negative allosteric modulator. The magnitude of conformational
changes gets increased as the information reaches near the intracellular region of TM
domain. These rearrangements have a very high activation barrier, which is affected by
the presence of intracellular coupling partners such as G-proteins or Arrestins
114
. Studies
84
have shown that interactions of a GPCR with a G-protein and agonist binding at the
orthosteric site act in a positive and cooperative way.
Clustering, based on UPGMA, shows the receptor conformations cluster based on their
activation states. Some microswitches/locations show higher degree of overlap between
the active, inactive and intermediate states than others. For example, the clustering in the
case of PIF motif has several intermediate states scattered, which is also evident from
the F
6.44
-S
3.39
side-chain COG distance-based histogram. The activation related
assignments are taken from the GPCRdb
8
. Comparing similar plots for the TM region
near the intracellular side shows a more apparent separation of the active, inactive and
intermediate states based on the activation coordinates taken from TM2 – TM6.
Collectively, these plots indicate the relative height of activation barrier, which is lower in
case of PIF motif compared to TM movements and signifies the importance of presence
of G-protein in order for the GPCR to undergo conformational transition into active state.
There is a caveat to such observations as these published structures may also under the
influence of mutations introduced to improve their stability for experimental
(crystallization/cryo-EM) outcomes.
85
Figure 34 Transmembrane region features clustering (a) histogram showing structures
with TM2 – TM6 distances (b) schematic view of the energy barrier requires for activation
and (c) UPGMA based distance matrix plotted for class A GPCR structures
86
Figure 35 PIF motif region features clustering (a) histogram showing structures with F
6.44
-
S
3.39
distances (b) schematic view of the energy barrier requires for activation and (c)
UPGMA based distance matrix plotted for class A GPCR structures
87
Figure 36 Distance matrices for (a) DRY motif and (b)the sodium ion binding site
88
Figure 37 Distance matrices for (a) NPxxY motif and (b) feature space including all
microswitches and activation related structural features
89
4.2.1 Case study: Structural analysis of DOR
We applied the GAUGE tool for prediction of the conformational state of recently resolved
delta-opioid receptor structures (PDB: 6PT2 and 6PT3). As previously discussed in
chapter 2, these structures are in complex with agonist ligands at the orthosteric site and
without any intracellular coupling protein. These structures have several activation related
structural features in the TM domain similar to the active state structures.
Figure 38 Activation-related changes in the DOP. Comparison of conserved activation
microswitches of the active-like DOP-KGCHM07 (orange) and DOP-DPI-287 (blue)
structures with the inactive DOP-naltrindole structure (yellow, PDB 4N6H). Structural
superposition of the (A) overall architecture, (B) PIF motif, (C) NPxxY and DRY motifs,
and (D) CWxP motif. (reprinted)
Predictions using the GAUGE tool assigns activation scores indicating that PIF, NPxxY,
and TM have features similar to that of other active state opioid receptor structures. The
extent of TM6 movement is not as pronounced as observed in the other active state opioid
90
structures in complex with G-protein or nanobodies stabilizing the conformation. Thus,
NPxxY and TM region are assigned to a score of 0.35 and 0.39, respectively. Moreover,
the lack of coupling partners at the intracellular side rendered DRY motif in the inactive
like state. These DOR structures provide an example of how GAUGE tool can be useful
in identification of crucial activation related features.
Table 4 Activation state-related score predictions of DOR structures using GAUGE. The
assessment shows all activation switches “on”, except the DRY motif, which remains in
an inactive conformation. Receptor names and corresponding PDB IDs are colored based
on their activation state (blue= inactive state; orange=active state), whereas the predicted
activation states for the receptor features PIF, DRY, NPxxY, and transmembrane helix
regions are colored based on the predicted activation state values (blue=inactive state;
red=active state). The C_score and R_score show the classification and regression
values, respectively, for a given structural feature (1=active; -1=inactive).
4.3 Evaluation of MD simulation trajectories
In addition to the experimental techniques such as X-ray crystallography, NMR, and Cryo-
EM, the conformational landscape of class A has been analyzed extensively using
91
computational MD simulations. Recent efforts focusing on streamlining the analysis and
making it more readily available to the non-experts in the simulation field have led to the
development of web portals for MD simulation repositories, visualization, and analysis
115
.
The goal of such endeavors is to gain insights into the conformational landscape in the
context of GPCR activation mechanism. The GAUGE tool could provide insights into
these dynamic processes of GPCR activation from the MD simulation trajectories in an
easy to interpret, and yet, robust way. Moreover, this further provides some valuable
insights into the progression and correlation of activation events throughout the key
activation related microswitches and regions.
4.3.1 Case study: MD analysis of β2AR active and inactive states
Beta-2 adrenergic receptor (β2-AR) is a class A GPCR protein, widely studied in terms of
its pharmacology, structural and computational biology. It was one of the first GPCRs to
be resolved experimentally using X-ray crystallography in 2007
6
. Currently, more than
twenty structures are available in PDB for β2-AR. There are several MD simulation studies
reported in the literature aiming to understand the ligand-receptor interactions and
activation mechanism of β2-AR. We analyzed some of the publicly accessible MD
simulation trajectories made available at the GPCRmd repository
115
. The details of
selected trajectories are given in Table 5
116,117
.
Table 5 MD simulation analysis using β2AR trajectories
92
The active state complex of β2-AR has a high-affinity catecholamine agonist HBI (hydroxy
benzyl isoproterenol), whereas the inactive state structure complexed with a high-affinity
antagonist /partial inverse agonist Carazolol
117
.
Figure 39 β2AR MD simulation analysis using the GAUGE analysis tool. Results are
plotted for two separate simulation runs for active state and inactive state β2AR (PDB ids:
4LDL, 2RH1). The number of trajectory frames is plotted on X-axis and features on Y-axis.
93
There are two MD simulations trajectories analyzed here for each state, and their
corresponding plots for activation related assessments by GAUGE are shown in Figure
39 and 40. Each independent simulation run consists 2500 trajectory frames for every 0.2
ns leading to 500 ns of simulation time, resulting a cumulative simulation of 1 µs for each
state.
Figure 40 Histogram of frames with activation scores predicted using GAUGE for β2-AR
active and inactive state simulations. Score densities are plotted on Y-axis, whereas score
range (-1 to 1) on X-axis.
From the GAUGE analysis, it becomes evident that the TM region conformational
ensemble returns to an intermediate state in the absence of coupling G-protein. The PIF
microswitch undergoes frequent transitions from inactive – intermediate – active states
as a significant energetic contribution, perhaps coming from the orthosteric site and a
feature of ligand interactions. Despite the absence of G-protein, the NPxxY motif
predominantly explores the active like conformational space, while the Na
+
binding site
may undergo some local rearrangements leading to transitions in intermediate states.
94
Overall, in the absence of a coupling protein at the intracellular site, the overall receptor
displays conformational features similar to active like intermediate states.
4.4 Summary and future directions
Identifying the link between structure and function is critical to finding better ways to
modulate the functions. With the rapid growth in structural information of class A GPCRs
and their computational simulations, it becomes vital to find techniques to evaluate and
extract meaningful information, which is also accessible to researchers. Also, it is equally
important to establish a uniform platform where the structure evaluation is performed in a
data-driven and automated way. Alongside the tools currently available for structure
visualization and analysis, the GAUGE is intended to provide a novel, uniform, and easily
accessible way to analyze class A GPCR structures.
With the expanding structural repertoire, it would always be possible to incorporate the
new structural information into such a prediction tool in order to make it more robust.
Predictions from GAUGE also shed light on conformations in the MD simulation trajectory,
which are unusual or belong to the conformational space that is less sampled. Thus, it
would also be valuable as a conformation selector to guide conformational sampling in
MD simulations.
95
5 Conclusion
G-protein coupled receptors are the most versatile players in the cell for interacting with
a variety of stimuli or signaling molecules. Their conformational landscape provides
physical means to the information transfer across the cell membrane. In the past decade,
advances in the structural biology of GPCRs have shed light into the atomic-level details
of these protein structures. Integrative computational approaches have become critical
for the integration, processing, and analysis of vast amount of experimental and
computational data. Development of special-purpose computing infrastructures (Anton
supercomputer, NVIDIA GPUs, cloud computing) have provided unprecedented speed of
performing large scale virtual screening and atomistic simulations, and as a result, drug
discovery quests are becoming increasingly powered by computational predictions.
In the work presented here, I applied molecular docking and energy based structural
refinement and modeling to gain insights into the ligand recognition at the melatonin and
the opioid receptors. Further, I performed large-scale virtual screening campaigns which
led to the discovery of novel and potent ligands of class A GPCRs. On the other hand, I
used MD simulations to understand the structure-function relationship of class A GPCRs.
MD simulations of the Angiotensin II receptors have provided key insights into the
receptor’s self-inhibitory mechanism. These insights could also help understand similar
structural features responsible for atypical signaling profiles in other class A GPCRs. In
the case of CysLT receptors, MD simulations shed light on the potential mode of lipid-
receptor interactions and the structural features facilitating such interactions. Finally,
integrating advanced machine learning method with molecular modeling, I prepared a tool
96
– GAUGE, for automated structural analysis and activation related scoring of class A
GPCRs. This tool can provide insights into an experimentally resolved receptor structure
or a model, as well as it can evaluate MD simulation trajectory frames to gain insights into
the conformational transitions in the context of receptor activation. The work presented
here has shown examples of development and applications of computational methods or
tools which can be applied to study 3D structures or models of class A GPCRs, their link
to the receptor activation, and for exploring the chemical space of class A GPCRs to find
drug like candidates or tool compounds.
97
References
(1) Fredriksson, R.; Schiöth, H. B. The Repertoire of G-Protein-Coupled Receptors in
Fully Sequenced Genomes. Mol. Pharmacol. 2005, 67 (5), 1414–1425.
https://doi.org/10.1124/mol.104.009001.
(2) Schiöth, H. B.; Lagerström, M. C. Structural Diversity of g Proteincoupled Receptors
and Significance for Drug Discovery. Nat. Rev. Drug Discov. 2008, 7 (4), 339–357.
https://doi.org/10.1038/nrd2518.
(3) Katritch, V.; Cherezov, V.; Stevens, R. C. Structure-Function of the G Protein–
Coupled Receptor Superfamily. Annu. Rev. Pharmacol. Toxicol. 2013, 53 (1), 531–
556. https://doi.org/10.1146/annurev-pharmtox-032112-135923.
(4) Hauser, A. S.; Attwood, M. M.; Rask-Andersen, M.; Schiöth, H. B.; Gloriam, D. E.
Trends in GPCR Drug Discovery: New Agents, Targets and Indications. Nat. Rev.
Drug Discov. 2017, 16 (12), 829–842. https://doi.org/10.1038/nrd.2017.178.
(5) Stauch, B.; Cherezov, V. Serial Femtosecond Crystallography of G Protein–
Coupled Receptors. Annu. Rev. Biophys. 2018, 47 (1), 377–397.
https://doi.org/10.1146/annurev-biophys-070317-033239.
(6) Cherezov, V.; Rosenbaum, D. M.; Hanson, M. A.; Rasmussen, S. G. F.; Foon, S.
T.; Kobilka, T. S.; Choi, H. J.; Kuhn, P.; Weis, W. I.; Kobilka, B. K.; et al. High-
Resolution Crystal Structure of an Engineered Human Β2-Adrenergic G Protein-
Coupled Receptor. Science (80-. ). 2007, 318 (5854), 1258–1265.
https://doi.org/10.1126/science.1150577.
98
(7) GPCRdb www.gpcrdb.org.
(8) Pándy-Szekeres, G.; Munk, C.; Tsonkov, T. M.; Mordalski, S.; Harpsøe, K.; Hauser,
A. S.; Bojarski, A. J.; Gloriam, D. E. GPCRdb in 2018: Adding GPCR Structure
Models and Ligands. Nucleic Acids Res. 2018, 46 (D1), D440–D446.
https://doi.org/10.1093/nar/gkx1109.
(9) Shimada, I.; Ueda, T.; Kofuku, Y.; Eddy, M. T.; Wüthrich, K. GPCR Drug Discovery:
Integrating Solution NMR Data with Crystal and Cryo-EM Structures. Nat. Rev.
Drug Discov. 2018, 18 (1), 59–82. https://doi.org/10.1038/nrd.2018.180.
(10) Brosey, C. A.; Tainer, J. A. Evolving SAXS Versatility: Solution X-Ray Scattering
for Macromolecular Architecture, Functional Landscapes, and Integrative Structural
Biology. Curr. Opin. Struct. Biol. 2019, 58, 197–213.
https://doi.org/10.1016/j.sbi.2019.04.004.
(11) Weis, W. I.; Kobilka, B. K. The Molecular Basis of G Protein–Coupled Receptor
Activation. Annu. Rev. Biochem. 2018, 87 (1), 897–919.
https://doi.org/10.1146/annurev-biochem-060614-033910.
(12) Seidel, L.; Zarzycka, B.; Zaidi, S. A.; Katritch, V.; Coin, I. Structural Insight into the
Activation of a Class B G-Protein-Coupled Receptor by Peptide Hormones in Live
Human Cells. Elife 2017, 6, 1–25. https://doi.org/10.7554/eLife.27711.
(13) Wei, G.-W. Protein Structure Prediction beyond AlphaFold. Nat. Mach. Intell. 2019,
1 (8), 336–337. https://doi.org/10.1038/s42256-019-0086-4.
(14) Gawehn, E.; Hiss, J. A.; Brown, J. B.; Schneider, G. Advancing Drug Discovery via
GPU-Based Deep Learning. Expert Opin. Drug Discov. 2018, 13 (7), 579–582.
99
https://doi.org/10.1080/17460441.2018.1465407.
(15) Fleming, N. How Artificial Intelligence Is Changing Drug Discovery. Nature 2018,
557 (7707), S55–S57. https://doi.org/10.1038/d41586-018-05267-x.
(16) Lyu, J.; Wang, S.; Balius, T. E.; Singh, I.; Levit, A.; Moroz, Y. S.; Meara, M. J. O.;
Che, T. Ultra-Large Library Docking for Discovering New Chemotypes. Nature.
https://doi.org/10.1038/s41586-019-0917-9.
(17) Abagyan, R.; Totrov, M. Biased Probability Monte Carlo Conformational Searches
and Electrostatic Calculations for Peptides and Proteins. J. Mol. Biol. 1994, 235,
983–1002.
(18) Halgren, T. A.; Murphy, R. B.; Friesner, R. A.; Beard, H. S.; Frye, L. L.; Pollard, W.
T.; Banks, J. L. Glide: A New Approach for Rapid, Accurate Docking and Scoring.
2. Enrichment Factors in Database Screening. J. Med. Chem. 2004, 47 (7), 1750–
1759. https://doi.org/10.1021/jm030644s.
(19) M. Mysinger, M.; K. Shoichet, B. Rapid Context-Dependent Ligand Desolvation in
Molecular Docking. J. Chem. Inf. Model. 2010, 50 (9), 1561–1573.
https://doi.org/10.1021/ci100214a.
(20) Trott, O.; Olson, A. J. AutoDock Vina: Improving the Speed and Accuracy of
Docking with a New Scoring Function, Efficient Optimization, and Multithreading. J.
Comput. Chem. 2010, 31 (2), 455–461. https://doi.org/10.1002/jcc.21334.
(21) Beuming, T.; Lenselink, B.; Pala, D.; McRobb, F.; Repasky, M.; Sherman, W.
Docking and Virtual Screening Strategies for GPCR Drug Discovery. In G Protein-
Coupled Receptors in Drug Discovery: Methods and Protocols; Filizola, M., Ed.;
100
Springer New York: New York, NY, 2015; pp 251–276. https://doi.org/10.1007/978-
1-4939-2914-6_17.
(22) Katritch, V.; Rueda, M.; Abagyan, R. Ligand-Guided Receptor Optimization. In
Homology Modeling: Methods and Protocols; Orry, A. J. W., Abagyan, R., Eds.;
Humana Press: Totowa, NJ, 2012; pp 189–205. https://doi.org/10.1007/978-1-
61779-588-6_8.
(23) Eswar, N.; Webb, B.; Marti-renom, M. A.; Madhusudhan, M. S.; Eramian, D.; Shen,
M.; Pieper, U.; Sali, A. Comparative Protein Structure Modeling Using Modeller;
2006. https://doi.org/10.1002/0471250953.bi0506s15.Comparative.
(24) Abagyan, R.; Orry, A.; Raush, E.; Totrov, M. ICM Manual, 3.8. MolSoft LLC. La
Jolla, CA 2016.
(25) Forli, S.; Huey, R. Computational Protein–Ligand Docking and Virtual Drug
Screening with the AutoDock Suite. Nat. Protoc. 2016, 11 (5), 905–919.
https://doi.org/10.1038/nprot.2016.051.
(26) Jones, G.; Willett, P. Development and Validation of a Genetic Algorithm for
Flexible Docking. J. Mol. Biol. 1997, 268, 727–748.
(27) Sterling, T.; Irwin, J. J. ZINC 15 - Ligand Discovery for Everyone. J. Chem. Inf.
Model. 2015, 55 (11), 2324–2337. https://doi.org/10.1021/acs.jcim.5b00559.
(28) Zheng, Z.; Huang, X.-P.; Mangano, T. J.; Zou, R.; Chen, X.; Zaidi, S.; Roth, B. L.;
Stevens, R. C.; Katritch, V. Structure Based Discovery of New Antagonist and
Biased Agonist Chemotypes for the Kappa Opioid Receptor. J. Med. Chem. 2017,
acs.jmedchem.7b00109. https://doi.org/10.1021/acs.jmedchem.7b00109.
101
(29) Zheng, Z.; Huang, X.; Mangano, T. J.; Zou, R.; Chen, X.; Zaidi, S. A.; Roth, B. L.;
Stevens, R. C.; Katritch, V. Structure-Based Discovery of New Antagonist and
Biased Agonist Chemotypes for the Kappa Opioid Receptor. 2017.
https://doi.org/10.1021/acs.jmedchem.7b00109.
(30) Nygaard, R.; Zou, Y.; Dror, R. O.; Mildorf, T. J.; Arlow, D. H.; Manglik, A.; Pan, A.
C.; Liu, C. W.; Fung, J. J.; Bokoch, M. P.; et al. The Dynamic Process of B2-
Adrenergic Receptor Activation. Cell 2013, 152 (3), 532–542.
https://doi.org/10.1016/j.cell.2013.01.008.
(31) Dror, R. O.; Dirks, R. M.; Grossman, J. P.; Xu, H.; Shaw, D. E. Biomolecular
Simulation: A Computational Microscope for Molecular Biology. Annu. Rev.
Biophys. 2012, 41, 429–452. https://doi.org/10.1146/annurev-biophys-042910-
155245.
(32) Luginina, A.; Gusach, A.; Marin, E.; Mishin, A.; Brouillette, R.; Popov, P.; Shiriaeva,
A.; Besserer-Offroy, É.; Longpré, J. M.; Lyapina, E.; et al. Structure-Based
Mechanism of Cysteinyl Leukotriene Receptor Inhibition by Antiasthmatic Drugs.
Sci. Adv. 2019, 5 (10). https://doi.org/10.1126/sciadv.aax2518.
(33) McGreevy, R.; Teo, I.; Singharoy, A.; Schulten, K. Advances in the Molecular
Dynamics Flexible Fitting Method for Cryo-EM Modeling. Methods 2016, 100, 50–
60. https://doi.org/10.1016/j.ymeth.2016.01.009.
(34) Kohlhoff, K. J.; Shukla, D.; Lawrenz, M.; Bowman, G. R.; Konerding, D. E.; Belov,
D.; Altman, R. B.; Pande, V. S. Cloud-Based Simulations on Google Exacycle
Reveal Ligand Modulation of GPCR Activation Pathways. Nat. Chem. 2014, 6 (1),
102
15–21. https://doi.org/10.1038/nchem.1821.
(35) Caliman, A. D.; Swift, S. E.; Wang, Y.; Miao, Y.; McCammon, J. A. Investigation of
the Conformational Dynamics of the Apo A2A Adenosine Receptor. Protein Sci.
2015, 24 (6), 1004–1012. https://doi.org/10.1002/pro.2681.
(36) Manglik, A.; Lin, H.; Aryal, D. K.; McCorvy, J. D.; Dengler, D.; Corder, G.; Levit, A.;
Kling, R. C.; Bernat, V.; Hübner, H.; et al. Structure-Based Discovery of Opioid
Analgesics with Reduced Side Effects. Nature 2016, 537 (7619), 185–190.
https://doi.org/10.1038/nature19112.
(37) Reppert, S. M.; Weaver, D. R.; Ebisawa, T. Cloning and Characterization of a
Mammalian Melatonin Receptor That Mediates Reproductive and Circadian
Responses. Neuron 1994, 13 (5), 1177–1185. https://doi.org/10.1016/0896-
6273(94)90055-8.
(38) Pévet, P. Melatonin Receptors as Therapeutic Targets in the Suprachiasmatic
Nucleus. Expert Opin. Ther. Targets 2016, 20 (10), 1209–1218.
https://doi.org/10.1080/14728222.2016.1179284.
(39) Brzezinski, A. Melatonin in Humans. N. Engl. J. Med. 1997, 336, 186–195.
https://doi.org/10.1056/NEJM199701163360306.
(40) Xie, Z.; Chen, F.; Li, W. A.; Geng, X.; Li, C.; Meng, X.; Feng, Y.; Liu, W.; Yu, F. A
Review of Sleep Disorders and Melatonin. Neurol. Res. 2017, 39 (6), 559–565.
https://doi.org/10.1080/01616412.2017.1315864.
(41) Ganguly, S.; Coon, S. L.; Klein, D. C. Control of Melatonin Synthesis in the
Mammalian Pineal Gland: The Critical Role of Serotonin Acetylation. Cell Tissue
103
Res. 2002, 309 (1), 127–137. https://doi.org/10.1007/s00441-002-0579-y.
(42) Stauch, B.; Johansson, L. C.; McCorvy, J. D.; Patel, N.; Han, G. W.; Huang, X. P.;
Gati, C.; Batyuk, A.; Slocum, S. T.; Ishchenko, A.; et al. Structural Basis of Ligand
Recognition at the Human MT 1 Melatonin Receptor. Nature 2019, 569 (7755),
284–288. https://doi.org/10.1038/s41586-019-1141-3.
(43) Johansson, L. C.; Stauch, B.; McCorvy, J. D.; Han, G. W.; Patel, N.; Huang, X. P.;
Batyuk, A.; Gati, C.; Slocum, S. T.; Li, C.; et al. XFEL Structures of the Human MT
2 Melatonin Receptor Reveal the Basis of Subtype Selectivity. Nature 2019, 569
(7755), 289–292. https://doi.org/10.1038/s41586-019-1144-0.
(44) Dubocovich, M. L.; Markowska, M. Functional MT 1 and MT 2 Melatonin Receptors
in Mammals. Endocrine 2005, 27 (2), 101–110.
https://doi.org/10.1385/ENDO:27:2:101.
(45) Hardeland, R.; Cardinali, D. P.; Srinivasan, V.; Spence, D. W.; Brown, G. M.; Pandi-
Perumal, S. R. Melatonin-A Pleiotropic, Orchestrating Regulator Molecule. Prog.
Neurobiol. 2011, 93 (3), 350–384. https://doi.org/10.1016/j.pneurobio.2010.12.004.
(46) Erman, M.; Seiden, D.; Zammit, G.; Sainati, S.; Zhang, J. An Efficacy, Safety, and
Dose-Response Study of Ramelteon in Patients with Chronic Primary Insomnia.
Sleep Med. 2006, 7 (1), 17–24. https://doi.org/10.1016/j.sleep.2005.09.004.
(47) Lavedan, C.; Forsberg, M.; Gentile, A. J. Tasimelteon: A Selective and Unique
Receptor Binding Profile. Neuropharmacology 2015, 91, 142–147.
https://doi.org/10.1016/j.neuropharm.2014.12.004.
(48) De Bodinat, C.; Guardiola-Lemaitre, B.; Mocaër, E.; Renard, P.; Muñoz, C.; Millan,
104
M. J. Agomelatine, the First Melatonergic Antidepressant: Discovery,
Characterization and Development. Nat. Rev. Drug Discov. 2010, 9 (8), 628–642.
https://doi.org/10.1038/nrd3140.
(49) Liu, J.; Clough, S. J.; Hutchinson, A. J.; Adamah-Biassi, E. B.; Popovska-Gorevski,
M.; Dubocovich, M. L. MT 1 and MT 2 Melatonin Receptors: A Therapeutic
Perspective . Annu. Rev. Pharmacol. Toxicol. 2015, 56 (1), 361–383.
https://doi.org/10.1146/annurev-pharmtox-010814-124742.
(50) López-Canul, M.; Comai, S.; Domínguez-López, S.; Granados-Soto, V.; Gobbi, G.
Antinociceptive Properties of Selective MT
2
Melatonin Receptor Partial
Agonists. Eur. J. Pharmacol. 2015, 764, 424–432.
https://doi.org/10.1016/j.ejphar.2015.07.010.
(51) Guillaume, J.-L.; Bonnefond, A.; Lukasheva, V.; Canouil, M.; Le Gouill, C.;
Langenberg, C.; Bouvier, M.; Wareham, N. J.; Lichtarge, O.; Plouffe, B.; et al. Type
2 Diabetes–Associated Variants of the MT 2 Melatonin Receptor Affect Distinct
Modes of Signaling . Sci. Signal. 2018, 11 (545), eaan6622.
https://doi.org/10.1126/scisignal.aan6622.
(52) Karamitri, A.; Jockers, R. Melatonin in Type 2 Diabetes Mellitus and Obesity. Nat.
Rev. Endocrinol. 2019, 15 (2), 105–125. https://doi.org/10.1038/s41574-018-0130-
1.
(53) Zlotos, D. P.; Jockers, R.; Cecon, E.; Rivara, S.; Witt-Enderby, P. A. MT1 and MT2
Melatonin Receptors: Ligands, Models, Oligomers, and Therapeutic Potential. J.
Med. Chem. 2014, 57 (8), 3161–3185. https://doi.org/10.1021/jm401343c.
105
(54) Jockers, R.; Delagrange, P.; Dubocovich, M. L.; Markus, R. P.; Renault, N.; Tosini,
G.; Cecon, E.; Zlotos, D. P. Update on Melatonin Receptors: IUPHAR Review 20.
Br. J. Pharmacol. 2016, 2702–2725. https://doi.org/10.1111/bph.13536.
(55) Zlotos, D. P. Recent Progress in the Development of Agonists and Antagonists for
Melatonin Receptors. Curr. Med. Chem. 2012, 19 (21), 3532–3549.
(56) Bateman, A.; Martin, M. J.; O’Donovan, C.; Magrane, M.; Alpi, E.; Antunes, R.; Bely,
B.; Bingley, M.; Bonilla, C.; Britto, R.; et al. UniProt: The Universal Protein
Knowledgebase. Nucleic Acids Res. 2017, 45 (D1), D158–D169.
https://doi.org/10.1093/nar/gkw1099.
(57) Ballesteros, J. A.; Weinstein, H. Integrated Methods for the Construction of Three-
Dimensional Models and Computational Probing of Structure-Function Relations in
G Protein-Coupled Receptors. Methods Neurosci. 1995, 25, 366–428.
https://doi.org/10.1016/S1043-9471(05)80049-7.
(58) Gaulton, A.; Hersey, A.; Nowotka, M. L.; Patricia Bento, A.; Chambers, J.; Mendez,
D.; Mutowo, P.; Atkinson, F.; Bellis, L. J.; Cibrian-Uhalte, E.; et al. The ChEMBL
Database in 2017. Nucleic Acids Res. 2017, 45 (D1), D945–D954.
https://doi.org/10.1093/nar/gkw1074.
(59) Kim, S.; Chen, J.; Cheng, T.; Gindulyte, A.; He, J.; He, S.; Li, Q.; Shoemaker, B. A.;
Thiessen, P. A.; Yu, B.; et al. PubChem 2019 Update: Improved Access to
Chemical Data. Nucleic Acids Res. 2019, 47 (D1), D1102–D1109.
https://doi.org/10.1093/nar/gky1033.
(60) Berman, H. M.; Westbrook, J.; Feng, Z.; Gilliland, G.; Bhat, T. N.; Weissig, H.;
106
Shindyalov, I. N.; Bourne, P. E. The Protein Data Bank Www.Rcsb.Org. Nucleic
Acids Res. 2000, 28 (1), 235–242. https://doi.org/10.1093/nar/28.1.235.
(61) Gatica, E. A.; Cavasotto, C. N. Ligand and Decoy Sets for Docking to G Protein-
Coupled Receptors. J. Chem. Inf. Model. 2012, 52 (1), 1–6.
https://doi.org/10.1021/ci200412p.
(62) Sterling, T.; Irwin, J. J. ZINC 15 - Ligand Discovery for Everyone. J. Chem. Inf.
Model. 2015, 55 (11), 2324–2337. https://doi.org/10.1021/acs.jcim.5b00559.
(63) Baell, J. B.; Holloway, G. a. New Substructure Filters for Removal of Pan Assay
Interference Compounds (PAINS) from Screening Libraries and for Their Exclusion
in Bioassays. J. Med. Chem. 2010, 53 (7), 2719–2740.
https://doi.org/10.1021/jm901137j.
(64) Abagyan, R.; Totrov, M.; Kuznetsov, D. ICM - A New Method for Protein Modeling
and Design: Applications to Docking and Structure Prediction from the Distorted
Native Conformation. J. Comput. Chem. 1994, 15 (5), 488–506.
https://doi.org/10.1002/jcc.540150503.
(65) Hopkins, A. L.; Colin, R. G.; Alexander, A. Ligand Efficiency : A Useful Metric for
Lead Selection. Drug Discov. Today 2004, 9 (10), 430–431.
(66) Hopkins, A. L.; Keserü, G. M.; Leeson, P. D.; Rees, D. C.; Reynolds, C. H. The
Role of Ligand Efficiency Metrics in Drug Discovery. Nat. Rev. Drug Discov. 2014,
13 (2), 105–121. https://doi.org/10.1038/nrd4163.
(67) Costanzi, S.; Vilar, S. In Silico Screening for Agonists and Blockers of the Β2
Adrenergic Receptor: Implications of Inactive and Activated State Structures. J.
107
Comput. Chem. 2012, 33 (5), 561–572. https://doi.org/10.1002/jcc.22893.
(68) Weiss, D. R.; Karpiak, J.; Huang, X. P.; Sassano, M. F.; Lyu, J.; Roth, B. L.;
Shoichet, B. K. Selectivity Challenges in Docking Screens for GPCR Targets and
Antitargets. J. Med. Chem. 2018, 61 (15), 6830–6845.
https://doi.org/10.1021/acs.jmedchem.8b00718.
(69) Roth, B. L. Molecular Pharmacology of Metabotropic Receptors Targeted by
Neuropsychiatric Drugs. Nat. Struct. Mol. Biol. 2019, 26 (7), 535–544.
https://doi.org/10.1038/s41594-019-0252-8.
(70) Audinot, V.; Mailliet, F.; Lahaye-Brasseur, C.; Bonnaud, A.; Le Gall, A.; Amossé,
C.; Dromaint, S.; Rodriguez, M.; Nagel, N.; Galizzi, J. P.; et al. New Selective
Ligands of Human Cloned Melatonin MT1 and MT2 Receptors. Naunyn.
Schmiedebergs. Arch. Pharmacol. 2003, 367 (6), 553–561.
https://doi.org/10.1007/s00210-003-0751-2.
(71) Kenakin, T. Biased Receptor Signaling in Drug Discovery. Pharmacol. Rev. 2019,
71 (2), 267–315. https://doi.org/10.1124/pr.118.016790.
(72) Cecon, E.; Oishi, A.; Jockers, R. Melatonin Receptors: Molecular Pharmacology
and Signalling in the Context of System Bias. Br. J. Pharmacol. 2018, 175 (16),
3263–3280. https://doi.org/10.1111/bph.13950.
(73) Wang, S.; Wacker, D.; Levit, A.; Che, T.; Betz, R. M.; McCorvy, J. D.;
Venkatakrishnan, A. J.; Huang, X. P.; Dror, R. O.; Shoichet, B. K.; et al. D4
Dopamine Receptor High-Resolution Structures Enable the Discovery of Selective
Agonists. Science (80-. ). 2017, 358 (6361), 381–386.
108
https://doi.org/10.1126/science.aan5468.
(74) Alkozi, H. A.; Montero, J. M. S.; Doadrio, A. L.; Pintor, J. Docking Studies for
Melatonin Receptors. Expert Opin. Drug Discov. 2018, 13 (3), 241–248.
https://doi.org/10.1080/17460441.2018.1419184.
(75) Fenalti, G.; Giguere, P. M.; Katritch, V.; Huang, X.-P.; Thompson, A. a; Cherezov,
V.; Roth, B. L.; Stevens, R. C. Molecular Control of $\delta$-Opioid Receptor
Signalling. Nature 2014, 506 (7487), 191–196.
https://doi.org/10.1038/nature12944.
(76) Fenalti, G.; Zatsepin, N. a; Betti, C.; Giguere, P.; Han, G. W.; Ishchenko, A.; Liu,
W.; Guillemyn, K.; Zhang, H.; James, D.; et al. Structural Basis for Bifunctional
Peptide Recognition at Human δ-Opioid Receptor. Nat. Struct. Mol. Biol. 2015, 22
(February), 1–6. https://doi.org/10.1038/nsmb.2965.
(77) Huang, W.; Manglik, A.; Venkatakrishnan, a J.; Laeremans, T.; Feinberg, E. N.;
Sanborn, A. L.; Gmeiner, P.; Kato, H. E.; Livingston, K. E.; Thorsen, T. S.; et al.
Structural Insights into Mu-Opioid Receptor Activation. Nature 2015, 524, 315–321.
https://doi.org/10.1038/nature14886.
(78) Koehl, A.; Hu, H.; Maeda, S.; Zhang, Y.; Qu, Q.; Paggi, J. M.; Latorraca, N. R.;
Hilger, D.; Dawson, R.; Matile, H.; et al. Structure of the μ-Opioid Receptor-G i
Protein Complex. Nature 2018, 558 (7711), 547–552.
https://doi.org/10.1038/s41586-018-0219-7.
(79) Che, T.; Majumdar, S.; Zaidi, S. A.; Ondachi, P.; McCorvy, J. D.; Wang, S.; Mosier,
P. D.; Uprety, R.; Vardy, E.; Krumm, B. E.; et al. Structure of the Nanobody-
109
Stabilized Active State of the Kappa Opioid Receptor. Cell 2018, 172 (1–2), 55-
67.e15. https://doi.org/10.1016/j.cell.2017.12.011.
(80) Liu, J. J.; Horst, R.; Katritch, V.; Stevens, R. C.; Wuthrich, K.; Wüthrich, K. Biased
Signaling Pathways in B2-Adrenergic Receptor Characterized by 19F-NMR.
Science 2012, 1106 (March), 1106–1111.
https://doi.org/10.1126/science.1215802.
(81) Katritch, V.; Fenalti, G.; Abola, E. E.; Roth, B. L.; Cherezov, V.; Stevens, R. C.
Allosteric Sodium in Class A GPCR Signaling. Trends Biochem. Sci. 2014, 39 (5),
233–244. https://doi.org/10.1016/j.tibs.2014.03.002.
(82) Venkatakrishnan, A. J.; Deupi, X.; Lebon, G.; Tate, C. G.; Schertler, G. F.; Madan
Babu, M. Molecular Signatures of G-Protein-Coupled Receptors. Nature 2013, 494
(7436), 185–194. https://doi.org/10.1038/nature11896.
(83) Guillemyn, K.; Kleczkowska, P.; Lesniak, A.; Dyniewicz, J.; Van Der Poorten, O.;
Van Den Eynde, I.; Keresztes, A.; Varga, E.; Lai, J.; Porreca, F.; et al. Synthesis
and Biological Evaluation of Compact, Conformationally Constrained Bifunctional
Opioid Agonist - Neurokinin-1 Antagonist Peptidomimetics. Eur. J. Med. Chem.
2015, 92, 64–77. https://doi.org/10.1016/j.ejmech.2014.12.033.
(84) Cottney, J.; Rankovic, Z.; Morphy, J. R. Synthesis of Novel Analogues of the Delta
Opioid Ligand SNC-80 Using REM Resin. Bioorganic Med. Chem. Lett. 1999, 9 (9),
1323–1328. https://doi.org/10.1016/S0960-894X(99)00173-0.
(85) Gutiérrez-De-Terán, H.; Massink, A.; Rodríguez, D.; Liu, W.; Han, G. W.; Joseph,
J. S.; Katritch, I.; Heitman, L. H.; Xia, L.; Ijzerman, A. P.; et al. The Role of a Sodium
110
Ion Binding Site in the Allosteric Modulation of the A2A Adenosine G Protein-
Coupled Receptor. Structure 2013, 21 (12), 2175–2185.
https://doi.org/10.1016/j.str.2013.09.020.
(86) Hori, T.; Okuno, T.; Hirata, K.; Yamashita, K.; Kawano, Y.; Yamamoto, M.; Hato,
M.; Nakamura, M.; Shimizu, T.; Yokomizo, T.; et al. Na + -Mimicking Ligands
Stabilize the Inactive State of Leukotriene B 4 Receptor BLT1. Nat. Chem. Biol.
2018, 14 (3), 262–269. https://doi.org/10.1038/nchembio.2547.
(87) Yi, S. P.; Kong, Q. H.; Li, Y. L.; Pan, C. L.; Yu, J.; Cui, B. Q.; Wang, Y. F.; Wang,
G. L.; Zhou, P. L.; Wang, L. L.; et al. The Opioid Receptor Triple Agonist DPI-125
Produces Analgesia with Less Respiratory Depression and Reduced Abuse
Liability. Acta Pharmacol. Sin. 2017, 38 (7), 977–989.
https://doi.org/10.1038/aps.2017.14.
(88) Halgren, T. A. Merck Molecular Force Field. I. Basis, Form, Scope,
Parameterization, and Performance of MMFF94. J. Comput. Chem. 1996, 17 (5–
6), 490–519. https://doi.org/10.1002/(SICI)1096-987X(199604)17:5/6<490::AID-
JCC1>3.0.CO;2-P.
(89) Zhang, H.; Han, G. W.; Batyuk, A.; Ishchenko, A.; White, K. L.; Patel, N.;
Sadybekov, A.; Zamlynny, B.; Rudd, M. T.; Hollenstein, K.; et al. Structural Basis
for Selectivity and Diversity in Angiotensin II Receptors. Nature 2017, 544 (7650).
https://doi.org/10.1038/nature22035.
(90) Zaman, M. A.; Oparil, S.; Calhoun, D. A. Drugs Targeting the Renin-Angiotensin-
Aldosterone System. Nat. Rev. Drug Discov. 2002, 1 (8), 621–636.
111
https://doi.org/10.1038/nrd873.
(91) Tate, C. G. A Receptor That Might Block Itself. Nature 2017, 544 (7650), 307–308.
https://doi.org/10.1038/nature21907.
(92) Zhang, H.; Unal, H.; Gati, C.; Han, G. W.; Liu, W.; Zatsepin, N. A.; James, D.; Wang,
D.; Nelson, G.; Weierstall, U.; et al. Structure of the Angiotensin Receptor Revealed
by Serial Femtosecond Crystallography. Cell 2015, 161 (4), 833–844.
https://doi.org/10.1016/j.cell.2015.04.011.
(93) Katritch, V.; Cherezov, V.; Stevens, R. C. Structure-Function of the G Protein-
Coupled Receptor Superfamily. Annu Rev Pharmacol Toxicol 2013, 53, 531–556.
https://doi.org/10.1146/annurev-pharmtox-032112-135923.
(94) Wacker, D.; Stevens, R. C.; Roth, B. L. How Ligands Illuminate GPCR Molecular
Pharmacology. Cell 2017, 170 (3), 414–427.
https://doi.org/10.1016/j.cell.2017.07.009.
(95) Kang, Y.; Zhou, X. E.; Gao, X.; He, Y.; Liu, W.; Ishchenko, A.; Barty, A.; White, T.
a.; Yefanov, O.; Han, G. W.; et al. Crystal Structure of Rhodopsin Bound to Arrestin
by Femtosecond X-Ray Laser. Nature 2015. https://doi.org/10.1038/nature14656.
(96) Venkatakrishnan, A. J.; Deupi, X.; Lebon, G.; Heydenreich, F. M.; Flock, T.; Miljus,
T.; Balaji, S.; Bouvier, M.; Veprintsev, D. B.; Tate, C. G.; et al. Diverse Activation
Pathways in Class A GPCRs Converge near the G-Protein-Coupling Region.
Nature 2016. https://doi.org/10.1038/nature19107.
(97) Guimond, M. O.; Gallo-Payet, N. How Does Angiotensin AT2 Receptor Activation
Help Neuronal Differentiation and Improve Neuronal Pathological Situations? Front.
112
Endocrinol. (Lausanne). 2012, 3 (DEC), 1–12.
https://doi.org/10.3389/fendo.2012.00164.
(98) Bäck, M.; Powell, W. S.; Dahlén, S. E.; Drazen, J. M.; Evans, J. F.; Serhan, C. N.;
Shimizu, T.; Yokomizo, T.; Rovati, G. E. Update on Leukotriene, Lipoxin and
Oxoeicosanoid Receptors: IUPHAR Review 7. Br. J. Pharmacol. 2014, 171 (15),
3551–3574. https://doi.org/10.1111/bph.12665.
(99) Luginina, A.; Gusach, A.; Marin, E.; Mishin, A.; Brouillette, R.; Popov, P.; Shiriaeva,
A.; Besserer-offroy, É.; Longpré, J.; Lyapina, E.; et al. Structure-Based Mechanism
of Cysteinyl Leukotriene Receptor Inhibition by Antiasthmatic Drugs. 2019, No.
October.
(100) Schönegge, A.-M.; Gallion, J.; Picard, L.-P.; Wilkins, A. D.; Le Gouill, C.; Audet, M.;
Stallaert, W.; Lohse, M. J.; Kimmel, M.; Lichtarge, O.; et al. Evolutionary Action and
Structural Basis of the Allosteric Switch Controlling Β2AR Functional Selectivity.
Nat. Commun. 2017, 8 (2169) (1), 1–12. https://doi.org/10.1038/s41467-017-
02257-x.
(101) Jo, S.; Kim, T.; Iyer, V. G.; Im, W. CHARMM-GUI: A Web-Based Graphical User
Interface for CHARMM. J. Comput. Chem. 2008, 29 (11), 1859–1865.
https://doi.org/10.1002/jcc.20945.
(102) Lomize, M. A.; Pogozheva, I. D.; Joo, H.; Mosberg, H. I.; Lomize, A. L. OPM
Database and PPM Web Server: Resources for Positioning of Proteins in
Membranes. Nucleic Acids Res. 2012, 40 (D1). https://doi.org/10.1093/nar/gkr703.
(103) Abraham, M. J.; Murtola, T.; Schulz, R.; Páll, S.; Smith, J. C.; Hess, B.; Lindahl, E.
113
GROMACS: High Performance Molecular Simulations through Multi-Level
Parallelism from Laptops to Supercomputers. SoftwareX 2015, 1–7.
https://doi.org/10.1016/j.softx.2015.06.001.
(104) Latorraca, N. R.; Venkatakrishnan, A. J.; Dror, R. O. GPCR Dynamics: Structures
in Motion. Chem. Rev. 2016, acs.chemrev.6b00177.
https://doi.org/10.1021/acs.chemrev.6b00177.
(105) Venkatakrishnan, A. J.; Fonseca, R.; Ma, A. K.; Hollingsworth, S. A.; Chemparathy,
A.; Hilger, D.; Kooistra, A. J.; Ahmari, R.; Babu, M. M.; Kobilka, B. K.; et al.
Uncovering Patterns of Atomic Interactions in Static and Dynamic Structures of
Proteins. bioRxiv 2019, 1–26. https://doi.org/10.1101/840694.
(106) Ibrahim, P.; Wifling, D.; Clark, T. Universal Activation Index for Class A GPCRs. J.
Chem. Inf. Model. 2019, 59 (9), 3938–3945.
https://doi.org/10.1021/acs.jcim.9b00604.
(107) Isberg, V.; De Graaf, C.; Bortolato, A.; Cherezov, V.; Katritch, V.; Marshall, F. H.;
Mordalski, S.; Pin, J. P.; Stevens, R. C.; Vriend, G.; et al. Generic GPCR Residue
Numbers - Aligning Topology Maps While Minding the Gaps. Trends Pharmacol.
Sci. 2015, 36 (1), 22–31. https://doi.org/10.1016/j.tips.2014.11.001.
(108) Gowers, R.; Linke, M.; Barnoud, J.; Reddy, T.; Melo, M.; Seyler, S.; Domański, J.;
Dotson, D.; Buchoux, S.; Kenney, I.; et al. MDAnalysis: A Python Package for the
Rapid Analysis of Molecular Dynamics Simulations. Proc. 15th Python Sci. Conf.
2016, No. Scipy, 98–105. https://doi.org/10.25080/majora-629e541a-00e.
(109) Sokal, R. R. A Statistical Method for Evaluating Systematic Relationships; 1958;
114
Vol. 38.
(110) Liu, W.; Chun, E.; Thompson, A. a; Chubukov, P.; Xu, F.; Katritch, V.; Han, G. W.;
Roth, C. B.; Heitman, L. H.; IJzerman, A. P.; et al. Structural Basis for Allosteric
Regulation of GPCRs by Sodium Ions. Science 2012, 337 (6091), 232–236.
https://doi.org/10.1126/science.1219218.
(111) Varoquaux, G.; Buitinck, L.; Louppe, G.; Grisel, O.; Pedregosa, F.; Mueller, A.
Scikit-Learn. GetMobile Mob. Comput. Commun. 2015, 19 (1), 29–33.
https://doi.org/10.1145/2786984.2786995.
(112) Van Der Walt, S.; Colbert, S. C.; Varoquaux, G. The NumPy Array: A Structure for
Efficient Numerical Computation. Comput. Sci. Eng. 2011, 13 (2), 22–30.
https://doi.org/10.1109/MCSE.2011.37.
(113) McKinney, W. Data Structures for Statistical Computing in Python. PROC. 9th
PYTHON Sci. CONF 2010, 51–56. https://doi.org/10.3828/ajfs.41.3.62.
(114) Latorraca, N. R.; Venkatakrishnan, A. J.; Dror, R. O. GPCR Dynamics: Structures
in Motion. Chem. Rev. 2017, 117 (1), 139–155.
https://doi.org/10.1021/acs.chemrev.6b00177.
(115) Ramírez-anguita, J. M. GPCRmd Uncovers the Dynamics of the 3D-GPCRome.
2019.
(116) Staus, D. P.; Strachan, R. T.; Manglik, A.; Pani, B.; Kahsai, A. W.; Kim, T. H.;
Wingler, L. M.; Ahn, S.; Chatterjee, A.; Masoudi, A.; et al. Allosteric Nanobodies
Reveal the Dynamic Range and Diverse Mechanisms of G-Protein-Coupled
Receptor Activation. Nature 2016, 535 (7612), 448–452.
115
https://doi.org/10.1038/nature18636.
(117) Ring, A. M.; Manglik, A.; Kruse, A. C.; Enos, M. D.; Weis, W. I.; Garcia, K. C.;
Kobilka, B. K. Adrenaline-Activated Structure of β 2-Adrenoceptor Stabilized by an
Engineered Nanobody. Nature 2013, 502 (7472), 575–579.
https://doi.org/10.1038/nature12572.
Appendix
A1. List of publications resulted from this doctoral research work
(1) Claff, T.; Yu, J.; Blais, V.; Patel, N.; Martin, C.; Wu, L.; Han, G. W.; Holleran, B. J.; Poorten,
O. Van Der; White, K. L.; et al. Elucidating the Active -Opioid Receptor Crystal Structure with
Peptide and Small-Molecule Agonists. Sci. Adv. 2019, 5 (November), eaax9115.
https://doi.org/10.1126/sciadv.aax9115.
In this project, I analyzed the delta-opioid receptor structure in complex with two agonists,
performed molecular modeling and docking to understand the ligand-receptor interactions, the
role of water network at the ligand binding pocket, and the activation state related feature analysis.
(2) Gusach, A.; Luginina, A.; Marin, E.; Brouillette, R. L.; Besserer-offroy, É.; Longpré, J.;
Ishchenko, A.; Popov, P.; Patel, N.; Fujimoto, T.; et al. Mutations in Cysteinyl Leukotriene
Receptors. Nat. Commun. 2019, 1–9. https://doi.org/10.1038/s41467-019-13348-2.
116
I performed molecular modeling and MD simulations to analyze the conformational landscape of
CysLT2R structures, focusing at the intracellular regions near the Helix-8 and activation related
TM-VII regions.
(3) Luginina, A.; Gusach, A.; Marin, E.; Mishin, A.; Brouillette, R.; Popov, P.; Shiriaeva, A.;
Besserer-Offroy, É.; Longpré, J. M.; Lyapina, E.; et al. Structure-Based Mechanism of Cysteinyl
Leukotriene Receptor Inhibition by Antiasthmatic Drugs. Sci. Adv. 2019, 5 (10).
https://doi.org/10.1126/sciadv.aax2518.
In this project, my work focused on evaluating the structural features observed in the CysLT 1R
structures and evaluate their biological relevance using MD simulation techniques.
(4) Stauch, B.; Johansson, L. C.; McCorvy, J. D.; Patel, N.; Han, G. W.; Huang, X. P.; Gati,
C.; Batyuk, A.; Slocum, S. T.; Ishchenko, A.; et al. Structural Basis of Ligand Recognition at the
Human MT 1 Melatonin Receptor. Nature 2019, 569 (7755), 284–288.
https://doi.org/10.1038/s41586-019-1141-3.
My contribution in this project was to analyze the structural features of MT 1 structure
experimentally solved using XFEL technique, analyze the interactions of co-crystallized and high-
affinity ligands using molecular modeling and docking, suggest mutations to interrogate the lipidic
channel observed in the structures, simulate and analyze the ligand interactions and structural
feature stabilities using MD simulations, evaluate the relationship between the ligands targeting
Melatonin and 5-HT2C receptors using structural assessment and molecular docking.
(5) Johansson, L. C.; Stauch, B.; McCorvy, J. D.; Han, G. W.; Patel, N.; Huang, X. P.; Batyuk,
A.; Gati, C.; Slocum, S. T.; Li, C.; et al. XFEL Structures of the Human MT 2 Melatonin Receptor
Reveal the Basis of Subtype Selectivity. Nature 2019, 569 (7755), 289–292.
https://doi.org/10.1038/s41586-019-1144-0.
117
In this project, I analyzed the ligand-receptor interactions using molecular modeling, docking, and
MD simulations to gain insights into the structural basis of ligand selectivity at the MT receptors.
(6) Eddy, M. T.; Gao, Z.-G.; Mannes, P.; Patel, N.; Jacobson, K. A.; Katritch, V.; Stevens, R.
C.; Wüthrich, K. Extrinsic Tryptophans as NMR Probes of Allosteric Coupling in Membrane
Proteins: Application to the A
2A
Adenosine Receptor. J. Am. Chem. Soc. 2018, 140.
https://doi.org/10.1021/jacs.8b03805.
Using the molecular modeling techniques to evaluate the stability of tryptophan mutations at the
critical structural regions in Adenosine A 2A receptor related to activation mechanism, we
suggested potential positions in the TM domain for NMR based eval uations of the receptor
activation process and ligand profiling.
(7) White, K. L.; Eddy, M. T.; Gao, Z.-G.; Han, G. W.; Lian, T.; Deary, A.; Patel, N.; Jacobson,
K. A.; Katritch, V.; Stevens, R. C. Structural Connection between Activation Microswitch and
Allosteric Sodium Site in GPCR Signaling. Structure 2018.
https://doi.org/10.1016/j.str.2017.12.013.
My contribution in this project is to analyze the structural features at the sodium binding site in
class A GPCRs using molecular modeling, in this case, for Adenosine A 2A receptor, to evaluate
its impact on the structural stability of the receptor.
(8) Ishchenko, A.; Wacker, D.; Kapoor, M.; Zhang, A.; Han, G. W.; Basu, S.; Patel, N.;
Messerschmidt, M.; Weierstall, U.; Liu, W.; et al. Structural Insights into the Extracellular
Recognition of the Human Serotonin 2B Receptor by an Antibody. Proc. Natl. Acad. Sci. U. S. A.
2017, 114 (31). https://doi.org/10.1073/pnas.1700891114.
In this project, I analyzed the 5-HT2B receptor structural interactions with a receptor subtype-
selective extracellular antibody using molecular modeling to understand the features from protein-
protein interaction underpinning the selectivity of these interactions.
118
(9) Zhang, H.; Han, G. W.; Batyuk, A.; Ishchenko, A.; White, K. L.; Patel, N.; Sadybekov, A.;
Zamlynny, B.; Rudd, M. T.; Hollenstein, K.; et al. Structural Basis for Selectivity and Diversity in
Angiotensin II Receptors. Nature 2017, 544 (7650), 327. https://doi.org/10.1038/nature22035
In this paper (as discussed in case study 1 in chapter 3), I analyzed the structural and interaction
stability using molecular modeling, docking and MD simulations for Angiotensin-II receptor Type2
and also performed biased sampling MD simulations to study helix-VIII interactions in order to
gain insights into its structural and biological relevance.
(10) Zhang, H.; Unal, H.; Desnoyer, R.; Han, G. W.; Patel, N.; Katritch, V.; Karnik, S. S.;
Cherezov, V.; Stevens, R. C. Structural Basis for Ligand Recognition and Functional Selectivity
at Angiotensin Receptor. J. Biol. Chem. 2015, 290 (49), 29127–29139.
https://doi.org/10.1074/jbc.M115.689000.
The focus of my work in this paper is on understanding the interactions of widely studied and
marketed drugs (belonging to -sartan series of antihypertensives) targeting Angiotensin-II type1
receptor using molecular modeling and docking approaches.
(11) Miller, R. L.; Thompson, A. A.; Trapella, C.; Guerrini, R.; Malfacini, D.; Patel, N.; Han, G.
W.; Cherezov, V.; Caló, G.; Katritch, V.; et al. The Importance of Ligand-Receptor Conformational
Pairs in Stabilization: Spotlight on the N/OFQ G Protein-Coupled Receptor. Structure 2015.
https://doi.org/10.1016/j.str.2015.07.024.
In this project, I performed molecular modeling and docking to evaluate the Nociceptin N/OFQ
receptor interactions with its ligands and their impact on the stability and crystallizability of such
complexes.
Manuscript in Preparation :
(1) Patel, N.; Katritch, V. GAUGE: An Automated, Machine Learning Based Tool for Activation
Based Class A GPCR Structural Analysisclass A.
119
(2) Patel, N.; Huang, X.P.; Grandner, J. M.; Johansson, L.C.; Stauch, B.; McCorvy, J. D.;
Roth, B. L.; Cherezov, V; Katritch, V. Structure-Based Discovery of Potent & Selective Melatonin
Receptor Agonists. eLife (under Review)
120
A2. GAUGE Tool for activation related conformation analysis
Installation requirements:
A. Download the GAUGE scripts from https://www.github.com/nil16/GAUGE
B. GAUGE is written in Python v. 3.6, and it requires following dependencies
- MDAnalysis (https://www.mdanalysis.org)
- Sci-kit learn (https://scikit-learn.org)
- Ipyvolume (https://github.com/maartenbreddels/ipyvolume)
- Chart-studio Plotly (https://chart-studio.plot.ly)
- Scipy (https://www.scipy.org)
- NumPy (https://numpy.org)
- Pandas (https://pandas.pydata.org)
- Matplotlib (https://matplotlib.org)
- Seaborn (https://seaborn.pydata.org)
Representative structures (PDB codes) for training the SVM classification models:
A. Actives (12)
- 5G53, 4MQS, 4MQT, 3POG, 4LDL, 5C1M, 2X72, 3SN6, 6B73, 3PXO, 2YDO,
5XRA
B. Inactive (18)
- 3EML, 4EIY, 2VT4, 4DKL, 1U19, 4IB4, 2R4R, 5JQH, 4DJH, 4N6H, 5UO9,
3UON, 2HPY, 1F88, 2G87, 3REY, 3RFM, 2RH1
Abstract (if available)
Linked assets
University of Southern California Dissertations and Theses
Conceptually similar
PDF
Structural studies of the human leukotriene B4 receptor 1
PDF
Structure and function of the orphan G protein-coupled receptor 6
PDF
Exploring the structure of G protein-coupled receptors
PDF
Computational analysis of the spatial and temporal organization of the cellular environment
PDF
Mapping 3D genome structures: a data driven modeling method for integrated structural analysis
PDF
Structure-based computational analysis and prediction of TCR CDR3 loops in the TCR-peptide-MHC complex using solvation parameters and peptide molecular dynamics.
PDF
Machine learning of DNA shape and spatial geometry
PDF
Nanomaterials under extreme environments: a study of structural and dynamic properties using reactive molecular dynamics simulations
PDF
Controlling electronic properties of two-dimensional quantum materials: simulation at the nexus of the classical and quantum computing eras
PDF
High throughput computational framework for synthesis and accelerated discovery of dielectric polymer materials using polarizable reactive molecular dynamics and graph neural networks
PDF
Photoexcitation and nonradiative relaxation in molecular systems: methodology, optical properties, and dynamics
PDF
Theoretical studies of lipid bilayer electroporation using molecular dynamics simulations
PDF
Shock-induced poration, cholesterol flip-flop and small interfering RNA transfection in a phospholipid membrane: multimillion atom, microsecond molecular dynamics simulations
PDF
Molecular dynamics studies of protein aggregation in unbounded and confined media
PDF
Quantum molecular dynamics and machine learning for non-equilibrium processes in quantum materials
PDF
Molecular simulations of water and monovalent ion dynamics in the electroporation of phospholipid bilayers
PDF
Profiling transcription factor-DNA binding specificity
PDF
Multi-scale quantum dynamics and machine learning for next generation energy applications
PDF
Integrative approach of biological systems analysis on regulatory pathways, modules, protein essentialities, and disease gene classification
PDF
Structure and kinetics of the Orb2 functional amyloid
Asset Metadata
Creator
Patel, Nilkanth Nareshbhai
(author)
Core Title
Structure – dynamics – function analysis of class A GPCRs and exploration of chemical space using integrative computational approaches
School
College of Letters, Arts and Sciences
Degree
Doctor of Philosophy
Degree Program
Molecular Biology
Publication Date
01/10/2021
Defense Date
12/02/2019
Publisher
University of Southern California
(original),
University of Southern California. Libraries
(digital)
Tag
computational structural biology,drug discovery,GPCR,machine learning,molecular dynamics,molecular modeling,OAI-PMH Harvest,structure-based drug design
Language
English
Contributor
Electronically uploaded by the author
(provenance)
Advisor
Katritch, Vsevolod (
committee chair
), Cherezov, Vadim (
committee member
), Fokin, Valery (
committee member
), Nakano, Aiichiro (
committee member
), Stevens, Raymond (
committee member
)
Creator Email
nilkantp@usc.edu
Permanent Link (DOI)
https://doi.org/10.25549/usctheses-c89-260928
Unique identifier
UC11673818
Identifier
etd-PatelNilka-8114.pdf (filename),usctheses-c89-260928 (legacy record id)
Legacy Identifier
etd-PatelNilka-8114.pdf
Dmrecord
260928
Document Type
Dissertation
Rights
Patel, Nilkanth Nareshbhai
Type
texts
Source
University of Southern California
(contributing entity),
University of Southern California Dissertations and Theses
(collection)
Access Conditions
The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...
Repository Name
University of Southern California Digital Library
Repository Location
USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA
Tags
computational structural biology
drug discovery
GPCR
machine learning
molecular dynamics
molecular modeling
structure-based drug design