SAR-Guided Scoring Function and Mutational Validation Reveal the Binding Mode of CGS-8216 at the α1+/γ2- Benzodiazepine Site

Article abstract of DOI:10.1021/acs.jcim.8b00199

The structural resolution of a bound ligand-receptor complex is a key asset to efficiently drive lead optimization in drug design. However, structural resolution of many drug targets still remains a challenging endeavor. In the absence of structural knowl

Full text of DOI:10.1021/acs.jcim.8b00199

Subscriber access provided by - Access paid by the | UCSB Libraries
Pharmaceutical Modeling
SAR-Guided Scoring Function and Mutational Validation Reveal
Binding Mode of CGS-8216 at the #1+/#2- Benzodiazepine Site
David C. B. Siebert, Marcus Wieder, Lydia Schlener, Petra Scholze, Stefan Boresch, Thierry
Langer, Michael Schnürch, Marko D. Mihovilovic, Lars Richter, Margot Ernst, and Gerhard F. Ecker
J. Chem. Inf. Model., Just Accepted Manuscript • DOI: 10.1021/acs.jcim.8b00199 • Publication Date (Web): 20 Jul 2018
Downloaded from http://pubs.acs.org on July 21, 2018
Just Accepted
“Just Accepted” manuscripts have been peer-reviewed and accepted for publication. They are posted
online prior to technical editing, formatting for publication and author proofing. The American Chemical
Society provides “Just Accepted” as a service to the research community to expedite the dissemination
of scientific material as soon as possible after acceptance. “Just Accepted” manuscripts appear in
full in PDF format accompanied by an HTML abstract. “Just Accepted” manuscripts have been fully
peer reviewed, but should not be considered the official version of record. They are citable by the
Digital Object Identifier (DOI®). “Just Accepted” is an optional service offered to authors. Therefore,
the “Just Accepted” Web site may not include all articles that will be published in the journal. After
a manuscript is technically edited and formatted, it will be removed from the “Just Accepted” Web
site and published as an ASAP article. Note that technical editing may introduce minor changes
to the manuscript text and/or graphics which could affect content, and all legal disclaimers and
ethical guidelines that apply to the journal pertain. ACS cannot be held responsible for errors or
consequences arising from the use of information contained in these “Just Accepted” manuscripts.
is published by the American Chemical Society. 1155 Sixteenth Street N.W.,
Washington, DC 20036
Published by American Chemical Society. Copyright © American Chemical Society.
However, no copyright claim is made to original U.S. Government works, or works
produced by employees of any Commonwealth realm Crown government in the course
of their duties.

Page 1 of 46
Journal of Chemical Information and Modeling
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
SARꢀGuided Scoring Function and Mutational
Validation Reveal Binding Mode of CGSꢀ8216 at the
α1+/γ2− Benzodiazepine Site
David C. B. Siebert1, ‡, Marcus Wieder2,3, ‡, Lydia Schlener2, Petra Scholze4, Stefan Boresch3, Thierry
Langer2, Michael Schnürch1, Marko D. Mihovilovic1, Lars Richter2,*, Margot Ernst5, Gerhard F. Ecker2
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
1Institute of Applied Synthetic Chemistry, TU Wien, Getreidemarkt 9/163, 1060 Vienna, Austria
2Department of Pharmaceutical Chemistry, University of Vienna, Althanstrasse 14, 1090 Vienna,
Austria
3Faculty of Chemistry, Department of Computational Biological Chemistry, University of Vienna,
Währingerstrasse 17, 1090, Vienna, Austria
4Department of Pathobiology of the Nervous System, Center for Brain Research, Medical University of
Vienna, Spitalgasse 4, 1090 Vienna, Austria
5Department of Molecular Neurosciences, Center for Brain Research, Medical University of Vienna,
Spitalgasse 4, 1090 Vienna, Austria
ACS Paragon Plus Environment
1

Journal of Chemical Information and Modeling
Page 2 of 46
ABSTRACT: The structural resolution of a bound ligandꢀreceptor complex is a key asset to efficiently
drive lead optimization in drug design. However, structural resolution of many drug targets still remains
a challenging endeavor. In the absence of structural knowledge, scientists resort on structureꢀactivity
relationship (SAR) to promote compound development. In this study we incorporated ligandꢀbased
knowledge to formulate a docking scoring function which evaluates binding poses for their agreement
with known SAR. We showcased this protocol by identifying the binding mode of the
pyrazoloquinolinone (PQ) CGSꢀ8216 at the benzodiazepine binding site of the GABAA receptor.
Further evaluation of the final pose by molecular dynamics and free energy simulations revealed a close
proximity between PQ’s pending phenyl ring and γ2D56 congruent with the low potency of carboxyꢀ
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
phenyl analogues. Ultimately, we introduced the γ2D56A mutation and in fact observed a 10ꢀfold
potency increase in the carboxyꢀphenyl analogue, providing experimental evidence in favour of our
binding hypothesis.
ACS Paragon Plus Environment
2

Page 3 of 46
Journal of Chemical Information and Modeling
INTRODUCTION
1
2
3
4
5
6
7
8
9
The structural resolution of ligandꢀreceptor complexes dramatically accelerates the elaborate
multiparameter optimization of a lead structure in drug discovery. The gained insight into binding mode
provides a rational basis to drive potency and selectivity improvements1. Moreover, medicinal chemists
can use the binding orientation of the ligand to optimize ADME parameters, such as solubility and
metabolic stability, without disrupting essential ligandꢀreceptor interactions1,2. In addition, this
information can assist to identify molecular determinants leading, for instance, to agonist and antagonist
behavior of ligands 3.
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
However, structural elucidation of a bound complex itself is a timeꢀconsuming endeavor and failed for
many important drug targets4. In the absence of protein structure information scientists can use different
methods to investigate lead compounds. Two prominent approaches are (1) homology modeling and
structureꢀbased design5, as well as (2) structureꢀactivity relationships6 (SAR) derived from biological
evaluation of compound libraries.
Among other factors, the success of homology modeling and subsequent structureꢀbased design
depends on the sequence similarity between target and template proteins and the resolution of the
template structure7. With decreasing structural certainty of the target protein, the chance to correctly
identify the correct binding mode declines8. In general, docking algorithms have been shown to be able
to explore the conformational space sufficiently well to generate the correct binding orientation,
however scoring functions often fail to correctly rank them9 . In this context the proper algorithmic
representation of complex processes like entropy, solvation and protein flexibility remains a big
challenge in pose scoring10,11. However, more elaborate methods like molecular dynamics (MD)
simulations and free energy simulations (FES) describe these processes more accurately and should
perform better in such specific cases12,13. Nevertheless, over the last years the Drug Design Data
Resource (D3R) challenges14 for proteinꢀligand pose and affinity predictions have shown impressively
that there is no one method 'to rule them all'. In the latest D3R Grand Challenge the two best performing
methods for pose prediction were molecular docking and MD simulations (according to the D3R Grand
ACS Paragon Plus Environment
3

Journal of Chemical Information and Modeling
Page 4 of 46
Challenge 2)15. Thus, a combined methodology implementing the best of both worlds might be a
promising approach.
1
2
3
4
5
6
7
8
9
In the absence of a reliable structural model, scientists focused on the identification of structureꢀ
activity relationships derived from biological profiling of scaffold analogues to drive lead optimization6.
In this context the GABAA receptors represent prominent examples of clinically important targets for
drugs like benzodiazepines (BZ), barbiturates, neuroactive steroids, anesthetics and anticonvulsants16. In
mammals, GABAA receptors are a heterogeneous group of pentameric receptors assembled from a pool
of nineteen possible subunits. In the central nervous system the majority of receptors is composed of
two α, two β and one γ subunits (Figure 1A)17. BZs exert their anxiolytic, muscleꢀrelaxant, sedativeꢀ
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
hypnotic and anticonvulsant effects by binding to an allosteric site located at the extracellular α+/γ−
subunit interface (Figure 1B, C)18. Generally, a high number of different chemotypes was developed
that targets various high affinity BZꢀbinding sites in the search for selective agents to be used, e.g. for
anxiolysis19. Exhaustive SARꢀanalyses were performed for several compound classes which led to an
high amount of pharmacological activity data20–23
.
Among the promising chemotypes, the PQs displayed high potential as nonꢀsedative anxiolytics and
as benzodiazepine antagonists respectively23. From this development two valuable tool compounds,
CGSꢀ8216 (Figure 1D) and its chloroꢀderivative CGSꢀ9896, emerged for in vitro and ex vivo work24.
Pharmacokinetic issues and inconclusive in vivo profiling of partial agonism and antagonism across
different GABAA receptor subtypes precluded further developments of these compounds into clinics25.
Recent work indicated that modulatory effects of many PQs in most αβγ receptors are predominantly
exerted via an alternative allosteric binding site at the “homologues” α+/β− interface26 while effects at
the high affinity BZꢀsites in most subtypes are “antagonistic”, i.e. flumazenilꢀlike27,28. Thus, revisiting
the PQs as valuable class of benzodiazepine antagonists due to their subꢀnanomolar affinity towards the
α+/γ− interfaces appears to be promising.
ACS Paragon Plus Environment
4

Page 5 of 46
Journal of Chemical Information and Modeling
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
Figure 1: (A) Side view and (B) top view of an α1β3γ2 GABAA receptor homology model based on the
human β3 homopentameric crystal structure 4COF without the intracellular domain. The agonist
binding sites (GABA) at the β3+/α1− interfaces and the benzodiazepine binding site (BZ site) at the
α1+/γ2− interface are indicated by arrows. α subunits are depicted in yellow, β subunits in red and the γ
subunit in blue. (C) α1+/γ2− extracellular interface shaped by loops AꢀC of the α+ side and loops DꢀG
of the γ2− side. (D) 2D representation of CGSꢀ8216.
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
CGSꢀ8216 (2ꢀphenylꢀ2Hꢀpyrazolo[4,3ꢀc]quinolinꢀ3(5H)ꢀone) has been widely used as an essential
tool compound for the benzodiazepine binding site. In the past 30 years, scientists have synthesized and
tested PQ analogues as BZꢀsite ligands leading to a highly differentiated SAR of this scaffold23,29
.
Interestingly, the PQs show both, characteristics of ‘continuous’ as well as ‘discontinuous’ SAR30,
depending on the respective scaffold substitution sites. Continuous SAR is characterized by a flat
activity hypersuperface ꢀ here even large structural changes on a molecular scaffold only slightly
influence the biological response. Discontinuous SAR, in contrast, is described by a hilly activity
hypersuperface, small structural changes can lead to drastic changes in potency.
ACS Paragon Plus Environment
5

Journal of Chemical Information and Modeling
Page 6 of 46
The described scenario, lowꢀreliability in protein structural information on the one hand and SAR
knowledge on the other hand is common for many scaffolds binding to membraneꢀbound drug targets.
While the combination of ligandꢀbased approaches, such as, 3Dꢀshapeꢀmatching, pharmacophore
1
2
3
4
5
6
7
8
mapping and 2DꢀQSAR with structureꢀbased approaches are widely applied in virtual screening31–33
,
9
they are less common to assist binding mode prediction. Varela et al34 utilized 3Dꢀsimilarity and
machine learning to augment binding mode prediction. In another study35, the authors coupled
electronic reactivity calculation, multipleꢀinstance learning and molecular docking to predict CYPꢀ
mediated sites of metabolism. However, most ligandꢀ, structureꢀbased combinations that rely on SAR
data from a congeneric series, analyze only a small number of poses, mostly a single, bestꢀscore pose36–
38. Hence, a potential pose with high SAR congruence but a slightly higher energetic scoring value
might be missed by such a procedure.
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Given the current challenges of molecular docking into homology models10 and the simultaneous
availability of a highly differentiated SAR for PQ23, we developed in this study a protocol which
integrates the ligand based knowledge in the process of docking pose evaluation. For this purpose, we
formulated a SAR scoring function which enabled the systematic assessment of a large CGSꢀ8216
docking pose library for their SAR congruency. The favorite binding hypothesis, characterized by a high
SARScore, was reconfirmed by sophisticated MD analysis composed of relative free energy calculations39
and stability analysis. Ultimately, we provide experimental evidence for our binding hypothesis by
synthesis and testing of PQ analogues in radioligand displacement assays. In line with our
computationally derived binding hypothesis, the PQꢀanalogue carrying a carboxyꢀphenyl moiety
showed a 10ꢀfold potency increase in the α1β3γ2ꢀD56A mutant compared to the α1β3γ2 receptor wild
type.
ACS Paragon Plus Environment
6

Page 7 of 46
Journal of Chemical Information and Modeling
1
2
3
METHODS
4
5
6
7
8
Homology model of the α1β3γ2 GABAA receptor. In order to determine the binding mode of PQs at
the extracellular α1+/γ2− interface (benzodiazepine binding site) a homology model based on 4COF40
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
was created consisting of two α1, two β3 and one γ2 subunits using Modeller41 as described
previously42. The transmembrane domains of the heteropentamer were truncated to obtain a model
consisting of the extracellular domains only in order to accelerate MD simulations.
Molecular Docking with GOLD. Molecular docking was performed using GOLD v1.6.243. The
putative binding site was defined by a cutꢀoff distance of 11.5 Å around the residue S204 of the Cꢀloop
of the α1 subunit. Further, we selected ten amino acids with flexible side chains (γ2Y58, γ2F77,
γ2T142, α1H101, α1Y159, α1V202, α1S204, α1S205, α1T206 and α1Y209) and set a soft potential to
increase to backbone flexibility of the Cꢀloop (α1S204, α1S205, α1T206 and α1G207). The ligand
CGSꢀ8216 was minimized using the MMFF94 force field44 within MOE2016.0845 prior to docking. To
ensure convergence of the sampling, 100 genetic algorithm runs were performed and the GoldScore46
scoring function was used to retain the 100 best scored poses. Finally, each ligandꢀreceptor complex of
the final pose pool was minimized using the LigX energy minimize tool within MOE to relieve potential
clashes caused by the soft potential protocol in the docking procedure.
Post-docking derivatization and SAR scoring. The postꢀdocking derivatization was conducted by a
customized SVL script using the MOEꢀSVL builtꢀin function sm_Build. This function enables the
addition and deletion of atoms on a molecular scaffold. On the basis of the 3D coordinates of each CGSꢀ
8216 docking pose (p1 ꢀ p100), an array of eight (2ꢀ9) analogue placements with sm_Build was created.
This step expanded the dataset from initial 100 CGSꢀ8216 docking poses to 800 additional postꢀdocking
derived analogue placements.
Steric clash energies (clashbb and clashsc) were calculated for the derivatized placements pi>>cpd 2-4
and pi>>cpd 6-9 using the MOE built in SVL function pro_CheckVDWContacts. The function returns
ACS Paragon Plus Environment
7

Journal of Chemical Information and Modeling
Page 8 of 46
clash scores defined as repulsion interaction energy (kcal/mol) for the respective analogue placement.
The Hꢀbond interaction strength (hbond) between the N5 nitrogen in CGSꢀ8216 docking poses and the
binding site were calculated using the SVL builtꢀin function prolig_Calculate.
1
2
3
4
5
6
7
8
9
Interatomic distance (d) calculation. The distances between atoms of the carboxylate moiety of
pi>>cpd 5 placements and the carboxylate group atoms of the closest aspartate or glutamate residue
were calculated using the SVL builtꢀin function aDist. From these distances the shortest interatomic
distance between the carboxylate group of the ligand and the receptor was kept.
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
The SAR scoring function is composed of four subfunctions: InAct.clash, N5.eval, pcarboxy.dist and
Act.clash. Each of these take the postꢀdocking engineered analogue geometries, pi>>cpd 2-9 as its input
(Figure S1) to evaluate the SAR agreement of a docking pose. To assess the SAR of steric hindrance
(R6, R6ꢀ7), the InAct.clash evaluates clash energies of pi>>cpd 2-3 placements (Figure S1B). As the
protocol does not include a minimization step after the postꢀdocking derivatization, clashes with rigid
structures such as backbone and Cβ atoms (clashbb) and sidechain atoms (clashsc) were weighted
distinctly. The N5.eval function evaluates both, potential steric hindrance or a disrupted hydrogen bond
interaction with quinoline nitrogen (Figure S1C). The PQꢀSAR suggests a strong electrostatic repulsion
between the 4’ꢀcarboxylate moiety and the pocket environment. The pcarboxy.eval function evaluates if
the 4’ꢀcarboxylate moiety of pi>>cpd 5 is placed in the vicinity (< 4 Å) of an aspartate or glutamate
carboxylate moiety (Figure S1D) by simple distance calculation. The steric tolerance at R7ꢀ8, R8, R8ꢀ9
and R4’ (Figure 2) was assessed via the Act.clash subfunction, utilizing the clash energies of pi>>cpd 6-
9 placements (Figure S1E).
Multidimensional Scaling (MDS) of the pose space. The SVLꢀexchange script MOL_RMSD.SVL47
was used to compute the RMSD matrix of the 100 CGSꢀ8216 docking poses. The matrix served as input
for MDS to visualize the geometric similarity between poses. The MDS was conducted within the R
environment48 utilizing the mds function of the smacof R package49. The first two dimensions were used
to visualize the pose space.
ACS Paragon Plus Environment
8

Page 9 of 46
Journal of Chemical Information and Modeling
Molecular dynamics simulations and postprocessing. The CHARMMꢀGUI web interface50 was used
1
2
3
to solvate the homology model of the intracellular GABAꢀligand complex and to set up the simulations.
4
5
All MD simulations were carried out with CHARMM51, utilizing the CHARMM/openMM coupling52,53
.
6
7
8
9
Parameters and molecular topologies for the ligands were generated based on the CGenFF force field54
obtained from ParamChem (https://cgenff.paramchem.org/). The proteinꢀligand complexes were
solvated in cubic boxes of TIP3P water55. Ions were added to compensate the net charge of the proteinꢀ
ligand complex and to set the ion concentration to 0.15 M KCI. Electrostatic interactions were
computed using the particleꢀmeshꢀEwald method56. SHAKE was used to keep all bonds involving
hydrogen atoms rigid.
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
After initial equilibration for 500 ps with a 1 fs time step, each system was simulated for 50 ns using
Langevin dynamics at 303.15 K; the pressure was kept around 1 atm by a Monte Carlo barostat. The
time step of the production calculation was 2 fs; coordinates were saved every 10 ps, resulting in 5,000
coordinate sets per simulation. For every system ten MD simulations started from the same initial
coordinates but with different, randomly seeded velocities were performed. The stability of the
simulations was monitored by computing root mean square deviations for the protein and ligand, using
the MDAnalysis package57, as well as visual inspection of the trajectories. The RMSDs were calculated
as follows: all coordinates saved during the MD were fitted against the starting structure based on the
coordinates of the Cαꢀatoms of the protein. Using the starting structure as reference, we computed for
these reoriented coordinates the RMSD of the Cαꢀatoms for the protein and the RMSD of the heavy
atoms of the ligand. For further calculation the average of the RMSD values of the ligand for the
different coordinate sets of a simulation is used (subsequently called LꢀRMSD in accordance with Liu et
al.58). The protein and ligand RMSD plots of the simulations (RMSD vs simulation time) are available
in the SI section. Two tests were performed using the LꢀRMSD: first, the average of the LꢀRMSD
values of the 10 independent simulations was calculated, i.e. one LꢀRMSD value was generated that
represents the 10 simulations and indicates the stability of the ligand for the binding ode. This value is
indicated in the manuscript by using angular brackets and typed as <LꢀRMSD>. Furthermore, the
ACS Paragon Plus Environment
9

Journal of Chemical Information and Modeling
Page 10 of 46
number of MD simulations that have a LꢀRMSD value below a threshold value (in this manuscript we
used 2.0 in accordance with Liu et al.58) were also evaluated. This test is a further indicator for the
stability of the ligand in the binding pocket ꢀ if the LꢀRMSD value of the majority of the simulations is
below the threshold a single simulation with a high LꢀRMSD value can be set in context.
Additionally, the types of proteinꢀligand interactions (ionic interactions, hydrophobic interactions,
aromatic and hydrogen bond interactions) were analysed for the MD simulations of the different poses
using an interaction matrix. This matrix was generated by obtaining a structureꢀbased pharmacophore
model for every saved coordinate set of the MD simulation for a specific binding site and subsequently
analysing the frequency of the individual features. The pharmacophore feature definitions are in
accordance with the definitions used by LigandScout59.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
The maximum possible frequency value corresponds to the length of the simulation, i.e., for 50 ns the
maximum number a specific feature can appear is 5,000. Since the ten MD simulations for each pose
were accumulated in a single interaction map, the maximum number an interaction could be formed was
50,000. The columns of the interaction matrix indicate all amino acid residues that are involved in an
interaction at some point during the MD simulation, the rows designate different parts of the ligand and
the color in the matrix indicate how often a specific amino acid was involved in a specific
pharmacophore feature. In this way it is possible to analyze the number of interaction partners and also
their frequency.
The naming for the columns was in accordance with the following conventions: the letter before the
colon indicate the type of interaction. H for hydrophobic, HBA and HBD for hydrogen bond
acceptor/donor, PI for positive ionizable interaction and AR for aromatic interaction. This was followed
by the one letter residue code and by the residue number of the amino acid. The naming for the rows are
in accordance to the labels using in Figure 1. Interactions that appear in fewer than 300 coordinate sets
(out of 50,000) were discarded. The interaction map was generated using the python package
matplotlib60 and the chemoinformatic toolkit CDPkit (https://github.com/aglanger/cdpkit).
ACS Paragon Plus Environment
10

Page 11 of 46
Journal of Chemical Information and Modeling
Free Energy Calculations. Relative free energy differences61,62 between 5 and 10 were computed
employing the usual thermodynamic cycle63 shown in Figure S2; i.e., we computed the alchemical free
energy difference between 5 and 10 in the binding pocket and in solution, respectively. The calculations
were carried out for each of the two binding modes considered. To compute the free energy differences
of interest, the use of intermediate states in addition to the physical end states of 5 and 10 was required
(Figure S3). The technical details and simulation protocols as well as individual results (Table S1) are
summarized in the SI section.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Benchmark molecular docking protocols. In total, six different docking protocols were utilized to
dock CGSꢀ8216 into the α1β3γ2 GABAA homology model. From each protocol the best scored pose
was kept for comparison studies. The best performing GoldScore46 and ChemScore46 poses were filtered
from the initial GOLD docking run. In MOE2016.08 the structure preparation protocol was utilized to
initialize the system. Then the InducedFit protocol within MOEdock was used to dock CGSꢀ8216 using
the default parameters. The best scored poses according to the London dG64 and the GBVI/WSA dG65
scoring function were kept. Within the Schrödinger Suite 2015.03 the Protein Preparation Wizard was
used to prepare the protein for Schrödinger’s Induced Fit docking protocol66 utilizing default parameters
including the Glide SP scoring function.
Molecular docking was also performed using Autodock Vina67 .For the molecular docking run
standard parameters were used, protein flexibility was regarded using flexible side chains (γ2Y58,
γ2F77, γ2T142, α1H101, α1Y159, α1V202, α1S204, α1S205, α1T206 and α1Y209) and the
exhaustiveness level was set to 16.
GABAA receptor subunits and mutated subunits. cDNA’s of rat GABAA receptor subunits α1, β3
and γ2S were cloned as described68. Mutant γ2SꢀD56A was constructed using the Q5 SiteꢀDirected
Mutagenesis Kit (New England Biolabs) following manufacturer’s instruction. We used the wildꢀtype
rat γ2ꢀpCI vector as template and the primers ATTCATACAGcTATGTACGTGAAC and
ACS Paragon Plus Environment
11

Journal of Chemical Information and Modeling
Page 12 of 46
TAATGTTGGTTTCACTCC resulting in a substitution of amino acid D56 (GAT) to A (GCT). The
mutated subunit was confirmed by sequencing.
1
2
3
4
5
6
7
8
9
Culturing of human embryonic kidney 293 cells. Human embryonic kidney (HEK) 293 cells
(American Type Culture Collection ATCC® CRLꢀ1574TM) were maintained in Dulbecco’s modified
Eagle medium (DMEM, high glucose, GlutaMAXTM supplement, Gibco 61965ꢀ059, ThermoFisher,
Waltham, Massachusetts, USA) supplemented with 10 % fetal calf serum (SigmaꢀAldrich F7524, St.
Louis, Missouria, USA), 100 U/ml PenicillinꢀStreptomycin (Gibco 15140ꢀ122, ThermoFisher,
Waltham, Massachusetts, USA) and MEM (NonꢀEssential Amino Acids Gibco 11140ꢀ035,
ThermoFisher, Waltham, Massachusetts, USA) on 10cm Cell culture dishes (Cell+, Sarstedt,
Nürnbrecht, Germany) at 37 °C and 5 % CO2. HEK293 cells were transfected with cDNAs encoding rat
GABAAꢀreceptor subunits subcloned into pCI expression vectors. The ratio of plasmids used for
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
69
transfection with the calcium phosphate precipitation method were 3 µg α1 : 3 µg β3 and 15 µg γ2
per 10 cm dish. Medium was changed 4ꢀ6 hours after transfection. Cells were harvested 72 days after
transfection by scraping into phosphate buffered saline. After centrifugation (10 min, 3 000 g, 4 °C)
cells were resuspended in TC50 (50 mM TrisꢀCitrate pH=7.1), homogenized with an ULTRAꢀ
TURRAX® (IKA, Staufen, Germany) and centrifuged (10 min, 3 000 g). Membranes were frozen at ꢀ20
°C until use.
Radioligand binding assay. [3H]Flunitrazepam (specific activity 83 Ci/mmmol) was purchased from
Perkin Elmer NEN (New England Nuclear) (Waltham, Massachusetts, USA). Diazepam (7ꢀchloroꢀ1,3ꢀ
dihydroꢀ1ꢀmethylꢀ5ꢀphenylꢀ2Hꢀ1,4, benzodiazepineꢀ2ꢀone) from Nycomed (Opfikon, Switzerland).
Standard chemicals came from SigmaꢀAldrich (St. Louis, Missouri, USA). Frozen membranes were
thawed, resuspended in TC50 and incubated for 90 min at 4 °C in a total of 500 µL of a solution
containing 50 mM Tris/citrate buffer, pH=7.1, 150 mM NaCl and 2 nM [3H]Flunitrazepam in the
absence of presence of either 5 µM diazepam (to determine unspecific binding) or various
concentrations of receptor ligands (dissolved in DMSO, final DMSOꢀconcentration 0.5 %). Membranes
were filtered through Whatman GF/B filters and the filters were rinsed twice with 4 mL of iceꢀcold 50
ACS Paragon Plus Environment
12

Page 13 of 46
Journal of Chemical Information and Modeling
mM Tris/citrate buffer. Filters were transferred to scintillation vials and subjected to scintillation
counting after the addition of 3 mL Ultima gold liquid scintillation cocktail. Nonspecific binding
determined in the presence of 5 µM Diazepam was subtracted from total [3H]Flunitrazepam binding to
result in specific binding. In order to determine the equilibrium binding constant KD of
[3H]Flunitrazepam for the various receptorꢀsubtypes, membranes were incubated with various
concentrations of [3H]Flunitrazepam in the absence or presence of 5 µM Diazepam. Saturation binding
experiments were analyzed using the equation Y=Bmax*X/(KD+X). Nonlinear regression analysis of
the displacement curves used the equation: log(inhibitor) vs. response ꢀ variable slope with Top=100 %
and Bottom=0 % Y=100/(1+10^((logIC50ꢀx)*Hillslope)). Both analyses were performed using
GraphPad Prism version 5.0a for Mac OS X, GraphPad Software, La Jolla California USA, www.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
graphpad.com. Drug concentrations resulting in half maximal inhibition of specific [3H]Flunitrazepam
69,70
binding (IC50) were converted to Ki values by using the ChengꢀPrusoff relationship
Ki=
IC50/(1+(S/KD)) with S being the concentration of the radioligand (2 nM) (KD values (MW ± SEM, n =
3ꢀ4): α1β3γ2SꢀWT = 7.2 ± 0.2 nM, α1β3γ2SꢀD56A = 10.3 ± 1.7 nM).
Synthesis. Commercially available reagents were used without further purification. Reactions were
monitored by thin layer chromatography with silica gel 60 F254 plates (E. Merck, Darmstadt, Germany).
1H and 13C NMR spectra were recorded on Bruker Avance Ultrashield 400 (1H: 400 MHz, 13C: 101
MHz) or Bruker Avance IIIHD 600 spectrometer equipped with a Prodigy BBO cryo probe (1H: 600
13
MHz, C: 151MHz). Chemical shifts are reported in parts per million (ppm) and were calibrated using
DMSOꢀd6 as internal standard. Multiplicities are denoted by s (singlet), br d (broad doublet), d
(doublet), dd (doublet of doublet) and m (multiplet). Melting points were determined with a Büchi
Melting Point Bꢀ545 apparatus. HRꢀMS was measured on an Agilent 6230 LC TOFMS mass
spectrometer equipped with an Aglient Dual AJS ESIꢀSource. The chromatographic purities of the
tested compounds were determined by HPLC analyses on a Nexera system from Shimadzu equipped
with 2 binary LCꢀ30AD pumps plus degassers, a CTOꢀ20 column oven and a SPDꢀM30A PDA
detector. Analytical separation was made using a KINETEX C18 column with 1.7 ꢁm pore size and 2.1
ACS Paragon Plus Environment
13

Journal of Chemical Information and Modeling
Page 14 of 46
x 50 mm length. As solvents HPLC grade water and acetonitrile plus 0.1% formic acid were used
(method: gradient of 5% à 95% acetonitrile in water, 12 min, 40°C, 0.5 mL flow, injection volume 3
µL). The purity of all final compounds was >95% (see SI for HPLC analysis).
1
2
3
4
5
6
7
8
9
Diethyl 2-(((4-methoxyphenyl)amino)methylene)malonate. 4ꢀMethoxyaniline (5 g, 36.9 mmol, 1 eq.)
and diethyl(ethoxymethylene)malonate (7.46 mL, 36.9 mmol, 1 eq.) were dissolved in toluene (50 mL)
and the reaction mixture was heated to reflux. After 22 h the solvent was removed under reduced
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
pressure and the residue was purified by column chromatography (gradient of 10%ꢀ30% EtOAc in LP)
1
to give the desired product as yellow crystals (7.90 g, 26.9 mmol, 73%). H NMR (400 MHz, CDCl
3
) δ
1.31 (t, J = 7.1 Hz, 3H), 1.37 (t, J = 7.1 Hz, 3H), 3.80 (s, 3H), 4.23 (q, J = 7.1 Hz, 2H), 4.29 (q, J = 7.1
Hz, 2H), 6.86 – 6.94 (m, 2H), 7.05 – 7.10 (m, 2H), 8.43 (d, J = 13.9 Hz, 1H), 10.98 (br d, J = 13.8 Hz,
1H). 13C NMR (101 MHz, CDCl
166.0, 169.4. TLC (PE/EtOAc = 3/1): R
Ethyl 6-methoxy-4-oxo-1,4-dihydroquinoline-3-carboxylate.
3
) δ 14.5, 14.6, 55.7, 60.1, 60.4, 92.6, 115.1, 119.0, 132.9, 152.8, 157.3,
f
= 0.63. Mp: 32ꢀ34 °C.
Diethyl
2ꢀ(((4ꢀ
methoxyphenyl)amino)methylene)malonate (2.5 g, 8.52 mmol) was dispersed in diphenylether (15 mL),
flushed with argon for 5 min and heated to 235 °C for 1 h. The reaction mixture was poured into LP, the
formed precipitate was collected by filtration and washed with LP/EtOAc (1/1, 3 x 40 mL) to yield the
desired product as brown powder (1.14 g, 4.62 mmol, 54%). 1H NMR (600 MHz, DMSOꢀd
) δ 1.27 (t, J
6
= 7.1 Hz, 3H), 3.84 (s, 3H), 4.20 (q, J = 7.1 Hz, 2H), 7.34 (dd, J = 8.9, 3.0 Hz, 1H), 7.56 (d, J = 3.0 Hz,
13
1H), 7.58 (d, J = 9.0 Hz, 1H), 8.49 (d, J = 6.7 Hz, 1H), 12.30 (br d, J = 6.7 Hz, 1H). C NMR (151
MHz, DMSOꢀd
TLC (2% MeOH in CH
Ethyl 4-chloro-6-methoxyquinoline-3-carboxylate. Ethyl 6ꢀmethoxyꢀ4ꢀoxoꢀ1,4ꢀdihydroquinolineꢀ3ꢀ
carboxylate (400 mg, 1.62 mmol) was dispersed in POCl (3 mL) and heated to reflux. After 2 h the
reaction mixture was poured onto ice, neutralized with satd. NaHCO , extracted with CH Cl (3 x 20
mL), washed with brine (1 x 20 mL), dried over Na SO , filtered and evaporated. The residue was
6
) δ 14.4, 55.5, 59.5, 105.5, 108.7, 120.6, 122.2, 128.5, 133.4, 143.7, 156.6, 165.0, 172.9.
2
Cl ): R = 0.33. Mp: 265ꢀ267 °C.
2
f
3
3
2
2
2
4
purified by FC (gradient of 5%ꢀ15% EtOAc in LP to give the desired product as colorless crystals (284
ACS Paragon Plus Environment
14

Page 15 of 46
Journal of Chemical Information and Modeling
mg, 1.08 mmol, 66%). 1H NMR (400 MHz, CDCl
) δ 1.30 – 1.56 (m, 3H), 3.92 (s, 3H), 4.43 (q, J = 7.1
3
1
2
3
4
5
Hz, 2H), 7.40 (dd, J = 9.2, 2.5 Hz, 1H), 7.53 (d, J = 2.5 Hz, 1H), 7.96 (d, J = 9.2 Hz, 1H), 8.98 (s, 1H).
13C NMR (101 MHz, CDCl
3
) δ 14.4, 55.9, 62.2, 103.1, 123.7, 125.1, 127.6, 131.7, 141.4, 145.9, 147.6,
6
7
8
159.5, 164.9. TLC (2% MeOH in CH
2
Cl
2
): R
f
= 0.84. Mp: 88ꢀ90 °C.
9
4-(8-Methoxy-3-oxo-3,5-dihydro-2H-pyrazolo[4,3-c]quinolin-2-yl)benzonitrile.
4ꢀ(8ꢀMethoxyꢀ3ꢀ
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
oxoꢀ3,5ꢀdihydroꢀ2Hꢀpyrazolo[4,3ꢀc]quinolinꢀ2ꢀyl)benzonitrile was synthesized according to literature
1
23,27,71 in 57% yield (yellow solid, 202 mg, 0.64 mmol). H NMR (400 MHz, DMSOꢀd6) δ 3.93 (s, 3H),
7.33 (dd, J = 9.1, 2.9 Hz, 1H), 7.60 (d, J = 2.9 Hz, 1H), 7.69 (d, J = 9.1 Hz, 1H), 7.87 – 7.95 (m, 2H),
8.45 – 8.52 (m, 2H), 8.73 (d, J = 6.4 Hz, 1H), 12.97 (br d, J = 6.2 Hz, 1H. 13C NMR (101 MHz, DMSOꢀ
d6) δ 55.7, 102.7, 104.6, 105.3, 118.1 (2C), 119.2, 119.8, 120.2, 121.4, 130.0, 133.2 (2C), 138.6, 143.5,
+
+
144.3, 157.7, 162.5. HRꢀMS: calculated [C18H13N4O2 ]: 317.1033; found [C18H13N4O2 ]: 317.1038
(diff.: −1.55 ppm). TLC (10% MeOH in CH2Cl2): Rf = 0.50. M.p.: decomposes > 300 °C.
4-(8-Methoxy-3-oxo-3,5-dihydro-2H-pyrazolo[4,3-c]quinolin-2-yl)benzoic acid. 4ꢀ(8ꢀMethoxyꢀ3ꢀ
oxoꢀ3,5ꢀdihydroꢀ2Hꢀpyrazolo[4,3ꢀc]quinolinꢀ2ꢀyl)benzonitrile (40 mg, 0.13 mmol) and NaOH (35 mg,
0.89 mmol) were dissolved in 3 mL EtOH/H2O (v/v) and the reaction mixture was heated to reflux.
After 18 h the mixture was acidified wit 2 M HCl and the precipitate was collected by filtration, washed
with water (3 mL), LP (15 mL), EtOAc (20 mL) and dried under reduced pressure to give 4ꢀ(8ꢀ
methoxyꢀ3ꢀoxoꢀ3,5ꢀdihydroꢀ2Hꢀpyrazolo[4,3ꢀc]quinolinꢀ2ꢀyl) benzoic acid as yellow solid (28 mg,
0.084 mmol, 66%). 1H NMR (400 MHz, DMSOꢀd6) δ 3.95 (s, 3H), 7.33 (dd, J = 9.1, 2.9 Hz, 1H), 7.62
(d, J = 2.9 Hz, 1H), 7.70 (d, J = 9.1 Hz, 1H), 8.04 (d, J = 8.5 Hz, 2H), 8.42 (d, J = 8.5 Hz, 2H), 8.73 (d,
J = 6.6 Hz, 1H), 12.79 (br s, 1H), 12.92 (br d, J = 6.6 Hz, 1H). 13C NMR (151 MHz, DMSOꢀd6) δ 55.7,
102.7, 104.9, 117.6 (2C), 119.9, 120.0, 121.4, 125.5, 129.9, 130.3 (2C), 138.3, 143.6, 143.9, 157.7,
+
+
162.3, 167.0. HRꢀMS: calculated [C18H14N3O4 ]: 336.0979; found [C18H14N3O4 ]: 336.0984 (diff.:
−1.46 ppm). TLC (20% MeOH in CH2Cl2): Rf = 0.57. M.p.: decomposes > 300 °C. HPLC: tR = 4.64 min
(>99.9% purity).
ACS Paragon Plus Environment
15

Journal of Chemical Information and Modeling
Page 16 of 46
2-(4-Aminophenyl)-8-methoxy-2,5-dihydro-3H-pyrazolo[4,3-c]quinolin-3-one. 2ꢀ(4ꢀAminophenyl)ꢀ
1
2
3
8ꢀmethoxyꢀ2,5ꢀdihydroꢀ3Hꢀpyrazolo[4,3ꢀc]quinolinꢀ3ꢀone was synthesized according to literature72 in
4
5
1
17% yield (yellow solid, 16 mg, 0.052 mmol). H NMR (600 MHz, DMSOꢀd6) δ 3.90 (s, 3H), 4.99 (br
6
7
8
9
s, 2H), 6.58–6.64 (m, 2H), 7.21 (dd, J = 9.0, 2.9 Hz, 1 H), 7.52 (d, J = 2.9 Hz, 1H), 7.63 (d, J = 9.0 Hz,
1H), 7.78–7.83 (m, 2H), 8.54 (s, 1H). 13C NMR (151 MHz, DMSOꢀd6) δ 55.6, 102.3, 105.4, 113.6 (2C),
119.2, 120.2, 120.9 (2C), 121.5, 129.9, 130.0, 137.6, 142.0, 145.6, 157.3, 160.6. HRꢀMS: calculated
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
+
+
[C17H15N4O2 ]: 307.1190; found [C17H15N4O2 ]: 307.1196 (diff.: −2.21 ppm). TLC (5% MeOH in
CH2Cl2): Rf = 0.25. M.p.: decomposes > 300 °C. HPLC: tR = 3.82 min (96.6% purity).
ACS Paragon Plus Environment
16

Page 17 of 46
Journal of Chemical Information and Modeling
RESULTS
1
2
3
4
5
6
7
8
9
Molecular docking of CGS-8216. As a first step in our workflow to identify the PQ binding mode
we docked the highly potent CGSꢀ821624 (cpd 1, Figure 1D) into the BZꢀsite of an α1β3γ2 GABAA
receptor homology model40 , using flexible side chains as well as soft potentials on the tip of the flexible
loop C (Figure 1C) (see Methods). The 100 best scored and minimized CGSꢀ8216 docking poses (p1 ꢀ
p100) were used for further evaluation.
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Structure-Activity Relationship of PQ (PQ-SAR). Based on the differentiated pool of SAR data for
the PQs we deduced a model which shows in the substitution positions R5, R6 and R4’ a discontinuous
SAR with activity cliffs while at positions R8 and R9 a continuous SAR is observed (Figure 2 top). In
detail, it has been shown that the introduction of substituents at position R6 leads to a dramatic loss of
potency as the substituent’s vanꢀderꢀWaals volume increases (Figure S4). For instance, the potency of
the R6 = CF3 analogue29 (2, Figure 2 bottom) drops by nearly 4 orders of magnitude compared to the
unsubstituted scaffold (1). Another considerable loss was described in the R5 = CH3 analogue 4. The
disruption of either a hydrogen bond or the introduction of a steric clash by the additional methyl group
might account to the drop of more than 2 log units in potency23. Interestingly, analogue 5 (R8 = OCF3,
R4’ = COOH) showed another dramatic potency cliff with IC50 = 3200 nM, whereas its amino analogue
10 (R8 = OCF3, R4’ = NH2) remains highly potent with an IC50= 0.80 nM. The drop might derive from
strong electrostatic repulsion between the carboxylate group in position R4’ and the surrounding binding
pocket.
In contrast to the rugged, discontinuous SAR landscape found at positions R5, R6 and R4’, at positions
R8 and R9 large structural changes have minimal influences on compound potency (Figure 2 bottom).
For example, the R8 = tert-butyl, R8,9 = benzofused and R7,8 = methylenedioxo analogues23,73 (6-8)
remain highly potent despite their large and rigid substituents.
ACS Paragon Plus Environment
17

Journal of Chemical Information and Modeling
Page 18 of 46
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
cpd
1
ref
R5
H
R6
H
R7
H
R8
H
R9
H
H
H
H
H
H
R4’
potency (nM)
Savini et al. 23
Savini et al. 29
Savini et al. 23
Savini et al. 23
Savini et al. 23
He at al. 73
H
IC50
IC50
IC50
IC50
IC50
=
=
>
=
=
0.4
1860
1000
1600
3200
0.1
H
CF3
H
H
H
2
H
6,7ꢀbenzofused
H
H
H
3
CH3
H
H
H
H
H
H
H
H
H
H
H
H
H
4
OCF3
tBu
COOH
H
5
H
Ki =
6
Savini et al.23
Anathan et al. 74
He at al. 73
H
8,9ꢀbenzofused
H
IC50
IC50
=
7.2
7
H
7,8ꢀOCH2O
H
H
H
H
=
3.0
8
H
H
H
H
C≡CCH3
NH2
Ki =
IC50
2.0
9
Savini et al. 29
H
OCF3
=
0.8
10
Figure 2. Top: 2D representation of PQꢀSARs. The letters A, B, C and D refer to the different rings in
the scaffold. The different positions and numbering of the different residues are shown by R5ꢀR9 and R4’
respectively. Positions which tolerate steric bulk (R8, R9 and R4’) are indicated by gray spheres, whereas
positions which do not tolerate steric bulk (R6) are indicated by red spheres. Electrostatic repulsion is
indicated by a red outline and prevented Hꢀbond formation by a light blue outline. Bottom: table of the
chemical structures and α1β3γ2 GABAA binding affinities of PQ cpds 1ꢀ10.
SAR guided pose selection protocol. The correct CGSꢀ8216 binding mode should align with the PQꢀ
SAR (Figure 2 top). Thus, we selected a set of 4 weak and 4 strong PQ binders (2-9, Figure 2 bottom)
that reflect the essentials of the PQꢀSARs. The correct binding mode should provide a rationale for the
ACS Paragon Plus Environment
18

Page 19 of 46
Journal of Chemical Information and Modeling
observed potency differences within these eight analogues. In detail, the binding orientation should
align with the concepts of i) steric hindrance at position R6 (2 and 3), ii) steric hindrance or loss of
hydrogen bonding in position R5 (4), iii) high affinity despite bulky substituents at position R8, R9 and
R4’ (6ꢀ9) and iv) negative electrostatic repulsion at position 4’ (9 and 10).
1
2
3
4
5
6
7
8
9
Considering the current difficulties of docking into homology models, we extended the conventional
docking by developing a protocol which evaluates docking poses for their agreement with PQꢀSAR, in
an automated manner. In the first, preparatory step, the post-docking derivatization tool, symbolized by
“>>“ in the following, takes the 3D coordinates of each retrieved CGSꢀ8216 (1) docking pose (p1ꢀp100)
and adds substituents to the PQ scaffold leading to postꢀdocking derivatized proteinꢀligand complexes
of analogues 2ꢀ9. Figure 3 displays the output of the derivatization tool for 2, 4 and 9 on docking pose
p88, resulting in the engineered complexes p88>>cpd 2, p88>>cpd 4 and p88>>cpd 9.
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
In the second, evaluative step, we defined a SAR scoring function which assesses the derivatized
analogue placements of a particular pose (pi>>cpd 2-9) for agreement with the PQꢀSAR (Figure S1).
The input of the SARScoring function takes calculated clash energies, Hꢀbond interaction strength
energies and distance calculations of analogues 2ꢀ9 as its input. For each analogue that aligns with the
PQꢀSAR, the function adds one point to the SARScore, resulting in a maximum score of eight points for a
given pose. For a detailed description of this method we refer to the Methods section.
ACS Paragon Plus Environment
19

Journal of Chemical Information and Modeling
Page 20 of 46
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Figure 3. Postꢀdocking derivatization and PQ analogue assessment. Original CGSꢀ8216 docking pose
p88 and its postꢀdocking derivatized analogue placements p88>>cpd 2, p88>>cpd 4, p88>>cpd 5 and
p88>>cpd 9. Clash energy (clashsc, clashbb), interatomic distances (d) and hydrogen bond strength
interactions (hbond) were calculated and serve as input for the SAR scoring function (see Methods and
Figure S1). For each analogue that aligns with the PQꢀSAR (Figure 2 top) one point is added to the
SARScore of a particular pose, leading to a maximum SARScore of eight. yellow, α1 subunit; blue, γ2
subunit; black spaceꢀfilling atoms, carbon atoms introduced in the postꢀdocking derivatization step;
clashbb, ligand clashes with backbone and Cβ atoms; clashsc, ligand clashes with sidechain atoms; d,
shortest distance between pi>>cpd 5 carboxylate group atoms and an ASP/GLU carboxylate group
atoms; hbond, Hꢀbond interaction strength between quinolone nitrogen and binding pocket (see Methods).
SAR scoring function identifies two candidate binding modes. To obtain an overview of the
geometric heterogeneity of the 100 CGSꢀ8216 docking poses, we applied multidimensional scaling
(MDS) methodology on the basis of the ligand RMSD distance matrix (see Methods). Docking poses
which are in vicinity to each other in the MDS plot share a similar binding mode (BM) while dissimilar
poses are distant to each other (Figure 4A). The MDS plot showed a KruksalꢀStress of 0.1975 which is
sufficient for visualization purpose. SAR assessment of 100 poses resulted in 14 poses with a SARScore
greater than 6 (Table S2). These poses can be grouped into two geometrically different binding modes,
BM I (Figure 4B) and BM II (Figure 4C). Twelve of the 14 best SAR scored poses belonged to BM I,
ACS Paragon Plus Environment
20

Page 21 of 46
Journal of Chemical Information and Modeling
containing two poses with the maximum SARScore of 8. The remaining two poses belonged to BM II,
with best score of 6 (Table S2, Figure S5).
1
2
3
4
5
6
7
The main geometric difference between BM I and BM II is the orientation of the CGSꢀ8216’s D ring.
In BM I, the carbonyl oxygen and the D ring point both to the γ2 subunit, while in BM II, they point
8
9
towards the α1 subunit. With respect to the interaction profiles of BM I and BM II, we found similar
interactions of their A and B rings. In both binding orientations the quinoline nitrogen (N5) is engaged
in Hꢀbond interactions with the backbone of α1Y159, while γ2F77 and α1Y209 show hydrophobic
and/or π−π interactions with the quinoline ring. In contrast, the interactions and position of the D ring
differ significantly. In BM I, the Dꢀring showed strong interactions with γ2Y58 and is pointing towards
the Fꢀloop (Figure 4B). In BM II it mainly interacts with α1H101 and is located close to loop B of the α
subunit (Figure 4C).
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Figure 4. Identification of two major binding modes BM I and BM II. (A) Visualization of the CGSꢀ
8216 pose space utilizing MDS (see Methods, Figure S5). The color of the dots reflects the agreement
of these poses with the SAR scoring from blue (low agreement) to red (high agreement). In the left
center of the MDS plot, a cluster of highly scored poses, containing two poses with a maximum
SARScore of 8, was found. These poses were defined as BM I. A second group of poses, in the lower
right corner of the MDS plot, showed a moderate score of 6. All other poses had a SARScore below 6
ACS Paragon Plus Environment
21

Journal of Chemical Information and Modeling
Page 22 of 46
points. A 2D representation of CGSꢀ8216 can be found in the upper left corner, with its rings labeled.
(B) Top and side view of the representative BM I pose p88, with the side chain rotamers as seen in this
particular pose. (C) Top and side view of the representative BM II p18, with the side chain rotamers as
seen in this particular pose. It should be noted that side chain rotamers cannot be predicted reliably from
homology models at this level of sequence similarity. yellow, α1 subunit; blue, γ2 subunit
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Analysis of BM I and BM II in the light of PQ-SAR. Two BM I poses (p88 and p64) displayed the
maximum SARScore of 8 and perfectly fit the PQꢀSAR. For characterizing BM I we selected p88 (Figure
5A) and its derivatized analogue placements, p88>>cpds 2-9 (Figure 5B). The low potency of 2ꢀ5 is
well reflected in the p88. The derivatized placements p88>>cpds 2-3 displayed severe backbone clashes
with α1Y159 and α1S158, congruent with the dramatic potency drop seen in R6 substituted analogues
(Figure S4). The strong Hꢀbond interaction between the quinolone nitrogen and the carbonyl oxygen of
α1Y159 in p88 aligns with the low potency of the Nꢀmethylated analogue 4. The carboxylate group of
the Dꢀring in p88>>cpd 5 was found to be close to γ2D56 and points towards loop F, on which γ2E189
and γ2D192 are positioned (Figure S6). This orientation is in good agreement with analogue’s 5 low
potency. Finally, postꢀengineered placements for the highly active, bulky R8 = tertbuyl, R8,9
=
benzofused, R7,8 = OCH2O and R4’ = C≡CCH3 (6-9) analogues (p88>>cpds 6-9) showed minor to no
clashes. Substitution sites R8 and R9, though still in the pocket already point towards the pocket entry
and bulk water. CGSꢀ8216’s ring D was placed below the tip of loop C and is partly solvent exposed.
Thus, analogue 9 with its rigid and bulky ꢀC≡CCH3 substituent at R4’ fits well into the pocket.
ACS Paragon Plus Environment
22

Page 23 of 46
Journal of Chemical Information and Modeling
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Figure 5: BM I in the light of the PQꢀSAR model. (A) Homology model of the high affinity binding
site at the extracellular α1+/γ2− interface showing p88 in BM I orientation (α1 subunit = yellow,
γ2 subunit = blue, pocket surface = gold grid). Focused perspective of the ligand orientation in the
pocket surface grid. (B) Table of weak and strong PQ binders in BM I orientation. Agreement of the
derivatized p88>>cpd 2-9 with the PQꢀSAR (Figure 2) is displayed by green hooks. Unfavorable
interactions such as and putative electrostatic repulsion or steric clashes are displayed by red crosses.
black spaceꢀfilling atoms, carbon atoms introduced in the postꢀdocking derivatization step (Figure 3);
clashbb, ligand clashes with backbone and Cβ atoms; clashsc, ligand clashes with sidechain atoms; d,
ACS Paragon Plus Environment
23

Journal of Chemical Information and Modeling
Page 24 of 46
shortest distance between pi>>cpd 5 carboxylate group and an ASP/GLU carboxylate group; hbond, Hꢀ
1
2
3
4
bond interaction strength between quinolone nitrogen and binding pocket (see Methods).
5
6
7
8
9
The best BM II pose displayed a SARScore of 6 and hence could not fully explain the PQꢀSAR. For
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
characterizing BM II we selected p18 (Figure S7) and its derivatized analogue placements, p18>>cpds 2-
9. In p18>>cpds 2-3 we found severe clashes with the side chain atoms of γ2M130 and γ2T142.
Likewise, to BM I, the quinoline nitrogen showed hydrogen bond interaction with the backbone
carbonyl oxygen of α1Y159. Hence, the potency loss seen in p18>>cpd 4 is plausible as well.
In contrast to BM I, BM II fails to provide a rationale for the dramatic loss in affinity seen in the
carboxy analogue 5. Among the closest residues to p18>>cpd 5 carboxyl group were the weak and
strong basic amino acids, α1H101 and α1K155 (Figure 4C). The vicinity of these basic amino acids to
ring D in p18>>cpd 5 is in contradiction to the assumption of strong negative electrostatic repulsion
between the analogue 5’s carboxylate moiety and the pocket environment. Finally, the four active
analogues 6-9 are well represented in BM II by p1>>cpds 6-9 (Figure S7).
Stability and affinity evaluation of 5 and 10 in BM I and BM II. BM I is in agreement with the PQꢀ
SAR, while BM II fails to provide a rational for the poor activity of 5. To investigate this further
additional analyses were conducted. We used the postꢀdocking derivatized complexes of analogues 5
and 10 for the best scored BM I (p64 and p88) as well as the BM II (p18 and p75) poses (Table S2, Figure
S5) to investigate the stability of the ligands using MD simulations and employing simple RMSD
criteria76. Additionally, the difference in relative binding free energy for 5 and 10 in BM I and BM II
was calculated.
The stability of the selected poses was evaluated using an adopted validation scheme described by Liu
et al.76. Ten independent MD simulations (same coordinates but different initial velocities) were
conducted for the derivatized BM I (p64>>cpd 5, p64>>cpd 10, p88>>cpd 5, p88>>cpd 10) and BM II
ACS Paragon Plus Environment
24

Page 25 of 46
Journal of Chemical Information and Modeling
(p18>>cpd 5, p18>>cpd 10, p75>>cpd 5, p75>>cpd 10) placements. To distinguish between stable and
unstable poses we considered the average ligandꢀRMSD of the ten MD simulations (<LꢀRMSD>) and
the number of MD simulations with a LꢀRMSD below 2 Å (Table 1 and Table S3). In addition to the
1
2
3
4
5
6
7
8
geometric analysis the free energy differences of each pi>>cpd 5 ꢀ pi>>cpd 10 pair were calculated
9
(Figure S2ꢀS3 and Table S1).
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Table 1. MDꢀbased stability analysis and interaction profile of 5 and 10 in BM I and BM II,
respectively (Figure S8ꢀ S9).
cpd 5 (R4’ = COOH)
BM I BM II
cpd 10 (R4’ = NH2)
BM I BM II
p88
2.3
p64
2.6
p18
1.7
p75
1.9
p88
1.5
p64
1.7
p18
2.4
p75
2.6
<LꢀRMSDa>
< 2.0 Åb >
7/10
2/10
8/10
7/10
10/10
8/10
3/10
2/10
α1F99, α1H101,
α1V202, α1Y209,
γ2F77
α1F99, α1H101,
α1V202, α1Y209,
γ2F77
α1Y159,
α1Y209, γ2F77
ABꢀringc
γ2F77, α1Y209
N5ꢀnitrogend
Cꢀringd
ꢀ
ꢀ
ꢀ
α1Y159
ꢀ
ꢀ
α1S204, γ2T142
α1S204, γ2T142
α1V202, γ2Y58,
α1H101
α1V202, γ2Y58,
α1H101, α1F99
Dꢀringc
γ2Y58
γ2Y58
γ2D56
Position R4’:
COOH/NH2
α1K155, γ2R194,
γ2R197
γ2K184, γ2R194
ꢀ
e
a average LꢀRMSD of 10 independent MD simulations.
b
Number of MD simulations for which the LꢀRMSD was below 2.0 Å
c Hydrophobic and aromatic interactions.
d Hydrogen bond interactions.
e Ionic interactions.
First, the stability and relative free energy difference for cpd 5 and 10 in BM I (p64 and p88) was
evaluated. For cpd 10 the <LꢀRMSD> values for the 10 MD simulations are 1.5 Å in p88 and 1.7 Å in p64
suggesting an overall stable complex for BM I, which is additionally reflected by the number of
simulations with an LꢀRMSD below 2 Å (18 out of 20 for both poses). For cpd 5 we obtained
ACS Paragon Plus Environment
25

Journal of Chemical Information and Modeling
Page 26 of 46
ambiguous results even though the <LꢀRMSD> of 2.3 Å in p88 and 2.6 Å in p64 are above the stability
criteria and thus suggesting a low stability of 5. However, considering the simulations in p88 we
observed only three simulations with an LꢀRMSD higher than 2 Å (LꢀRMSD between 4.3 Å and 5.6 Å)
putting this quantitative analysis in doubt. The relative binding free energy difference between 5 and 10
in BM I was calculated as ꢀ2.5 kcal/mol (standard deviation of 1.3 kcal/mol) reflecting the higher
affinity of 10. The proximity between the carboxylic group of 5 and γ2D56 seemingly leads to a
destabilization of the binding pose in BM I.
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Interestingly, in BM II a contrary stability profile for 5 and 10 (p18 and p75) was observed. For 5 <Lꢀ
RMSD> values of 1.7 Å in p18 and 1.9 Å in p75 were obtained and the majority of the individual MD
simulations were below 2.0 Å indicating a high stability of 5 in BM II. In contrast, 10 displayed <Lꢀ
RMSD> values of 2.4 Å in p18 and 2.6 Å in p4. In addition, most individual simulations (15 out of 20
simulations) showed LꢀRMSD values higher than 2.0 Å suggesting a low stability of 10 in BM II.
Moreover, the calculated difference in binding free energy between 5 and 10 is 0.02 kcal/mol (standard
deviation of 1.88 kcal/mol). These results provide further evidence against BM II as it failed to
rationalize the difference in potency between 5 and 10.
SAR guided docking vs. conventional molecular docking protocols – a comparison. We were
interested whether conventional docking protocols would favor a similar binding mode (BM I) as our
new SARꢀdriven protocol. For this purpose, we conducted molecular docking using six different
docking protocols. We performed flexible docking of CGSꢀ8216 into our α1β3γ2 GABAA homology
model utilizing MOE 2013.0845, GOLD 5.243, Schrödinger's Induced Fit66 and AutoDock Vina
1.1.22067. The poses within each run were ranked by the software’s internal scoring function
encompassing the GBVI/WSA dG65, London dG64, ChemScore46, GoldScore46, Glide SP77 and the
AutoDock Vina scoring function67. We kept the best scored pose per scoring function and compared
these orientations against p88 (BM I) from the SAR guided protocol. While AutoDock Vina (RMSD ~
1.5 Å) and Induced Fit protocol (RMDS ~ 2.8 Å) were able to identify similar binding modes as our
ACS Paragon Plus Environment
26

Page 27 of 46
Journal of Chemical Information and Modeling
SAR guided protocol, the other four protocols showed completely different orientations with RMSD
distances greater than 6 Å (Table S4).
1
2
3
4
5
6
7
8
Prospective γ2D56A mutant strengthens BM I binding hypothesis. Based on the conclusive results
of our in silico analysis we aimed to strengthen our binding hypothesis (BM I) by a prospective
mutation. In BM I, CGSꢀ8216’s pending phenyl ring is in the vicinity to three acidic residues, γ2D56,
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
γ2E189 and γ2D192 (Figure S6). We selected the residue closest to CGSꢀ8216 in BM I, γ2D56 and
converted it into an alanine residue (γ2D56A). This point mutation should lead to a loss of the putative
electrostatic repulsion feature of our previous analysis (Figure 2). To address this missing feature we
synthesized two ligands, namely compound 18 and 19, which differ from 5 and 10 only in the position
R8 (for 18 and 19 R8 = OCH3; for 5 and 10 R8 = OCF3) while possessing the same electrostatic
properties in position R4’ (Scheme 1). Converting the 4ꢀmethoxyaniline with diethylethoxymethylene
malonate via the GouldꢀJacobs78 reaction yielded the quinoline 14. Treatment with POCl3 led to 15
which was reacted further with 4ꢀnitrophenylhydrazine and 4ꢀcyanophenylhydrazine respectively to
obtain the precursors 16 and 17. The reduction of 16 yielded the desired amine 18 and basic hydrolysis
of 17 led to the corresponding carboxylic acid derivative 19 (Scheme 1).
ACS Paragon Plus Environment
27

Journal of Chemical Information and Modeling
Page 28 of 46
1
2
3
4
5
6
7
Scheme 1. Synthetic route to compound 18 and 19.
8
9
First, we tested our new PQ ligands in α1β3γ2 wild type receptors to verify that 19 with the
carboxylic acid moiety in position R4’ is a weak binder compared to 18 with the amino group in the
same position. As expected 18 (Ki = 0.21 ± 0.03 nM) displayed a three order of magnitude higher
potency compared to 19 (Ki = 161 ± 27 nM) (Figure 6, grey lines). Thus, we were able to investigate our
new compounds in the α1β3γ2ꢀD56A mutant construct. Remarkably, we observed for 19 an increase of
potency by a factor of higher than 10 (Ki = 15 ± 2 nM) due to the introduction of our single point
mutation whereas 18 showed only a very weak potency shift by a factor of 2 (Ki = 0.094 ± 0.005 nM)
(Figure 6, dotted lines). Hence, these results support our binding hypothesis that PQs favor BM I.
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
Figure 6: (A) [3H]Flunitrazepam displacement assay of 18 in α1β3γ2-WT receptors (ꢀ, Ki = 0.21 ± 0.03 nM) and
α1β3γ2S-D56A mutants (●, Ki = 0.094 ± 0.005 nM) (mean ± SEM, n = 3-4). (B) [3H]Flunitrazepam displacement
assay of 19 in α1β3γ2S-WT receptors (ꢀ, Ki = 161 ± 27 nM) and α1β3γ2S-D56A mutants (●, Ki = 15 ± 2 nM)
(mean ± SEM, n = 3-4). Drug concentrations resulting in half maximal inhibition of specific [3H]Flunitrazepam
binding (IC50) were converted to Ki values by using the Cheng-Prusoff relationship (Figure S10).
ACS Paragon Plus Environment
28

Page 29 of 46
Journal of Chemical Information and Modeling
DISCUSSION
1
2
3
4
5
6
7
8
9
The experimental elucidation of drugꢀtarget complexes is laborꢀintensive and even impossible for
many membraneꢀbound drug targets. Here, computational approaches such as homology modeling
combined with molecular docking can generate reasonable binding hypotheses, but frequently docking
scoring functions struggle to filter them from a pool of alternative orientations9.
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
To increase the credibility of docking pose selection for lowꢀreliability target structures, we
established in this study a fully automatized routine applicable to molecules for which a distinct SAR is
available. Specifically, we defined a docking scoring function that solely evaluates docking poses by
their degree of SAR agreement. To evaluate binding orientations, we derived analogue placements from
four weak (2-5) and four strong (6-9) binders (Figure 2 bottom) on the basis of the coordinates of each
CGSꢀ8216 docking pose (Figure 3). Subsequently, the SAR scoring function utilized clash analysis,
distance calculation and Hꢀbond interaction strength calculation to assess the congruence between the
analogue placement and the PQꢀSAR (see Methods and Figure S1).
Our SARꢀdriven docking protocol led to one favorable binding mode, BM I (Figure 4B), fully
compatible with the PQꢀSAR (Figure 5B), reflected by a maximum SARScore of eight. A second,
moderately performing binding mode, BM II (Figure 4C), failed to provide a structural evidence for the
assumed electrostatic repulsion between carboxy analogue (5) and the pocket environment (Figure S6).
However, the SAR scoring function estimates the electrostatic repulsion assumption by simple distance
calculations (see Methods).
To evaluate the electrostatic influence more sophisticatedly, we performed MD based stability
analysis for the two best SAR scored BM I and BM II poses and relative binding free energy
calculations for the carboxy analogue (5) and the amino analogue (10) in each binding mode. In BM I
the results of the relative binding free energy calculations as well as the stability analysis indicate that
10 is stable in contrast to 5 – thus BM I can give a rationale for the reported potency discrepancy
between these analogues while BM II failed. Here, 10 is more stable than 5 and there is no difference in
the relative binding free energy between 5 and 10.
ACS Paragon Plus Environment
29