Prospective Evaluation of Free Energy Calculations for the Prioritization of Cathepsin L Inhibitors

Article abstract of DOI:10.1021/acs.jmedchem.6b01881

Improving the binding affinity of a chemical series by systematically probing one of its exit vectors is a medicinal chemistry activity that can benefit from molecular modeling input. Herein, we compare the effectiveness of four approaches in prioritizing building blocks with better potency: selection by a medicinal chemist, manual modeling, docking followed by manual filtering, and free energy calculations (FEP). Our study focused on identifying novel substituents for the apolar S2 pocket of cathepsin L and was conducted entirely in a prospective manner with synthesis and activity determination of 36 novel compounds. We found that FEP selected compounds with improved affinity for 8 out of 10 picks compared to 1 out of 10 for the other approaches. From this result and other additional analyses, we conclude that FEP can be a useful approach to guide this type of medicinal chemistry optimization once it has been validated for the system under consideration.

Full text of DOI:10.1021/acs.jmedchem.6b01881

Article
pubs.acs.org/jmc
Prospective Evaluation of Free Energy Calculations for the
Prioritization of Cathepsin L Inhibitors
Bernd Kuhn,†,# Michal Tichy,‡,# Lingle Wang,§,# Shaughnessy Robinson,§ Rainer E. Martin,†
́
Andreas Kuglstatter,† Jorg Benz,† Maude Giroud,‡ Tanja Schirmeister,∥ Robert Abel,*,§
̈
Franco
̧
is Diederich,*,‡ and Jerome Hert*,†
†Roche Pharmaceutical Research and Early Development (pRED), Roche Innovation Center Basel, F. Hoffmann-La Roche Ltd.,
́
̂
Grenzacherstrasse 124, 4070 Basel, Switzerland
‡Laboratorium fur Organische Chemie, ETH Zurich, Vladimir-Prelog-Weg 3, 8093 Zurich, Switzerland
̈
§Schrodinger, Inc., 120 West 45th Street, New York, New York 10036, United States
̈
∥Institut fur Pharmazie und Biochemie, Johannes Gutenberg-Universitat Mainz, Staudinger Weg 5, 55128 Mainz, Germany
̈
̈
S
* Supporting Information
ABSTRACT: Improving the binding affinity of a chemical series by systematically probing one of its exit vectors is a medicinal
chemistry activity that can benefit from molecular modeling input. Herein, we compare the effectiveness of four approaches in
prioritizing building blocks with better potency: selection by a medicinal chemist, manual modeling, docking followed by manual
filtering, and free energy calculations (FEP). Our study focused on identifying novel substituents for the apolar S2 pocket of
cathepsin L and was conducted entirely in a prospective manner with synthesis and activity determination of 36 novel
compounds. We found that FEP selected compounds with improved affinity for 8 out of 10 picks compared to 1 out of 10 for the
other approaches. From this result and other additional analyses, we conclude that FEP can be a useful approach to guide this
type of medicinal chemistry optimization once it has been validated for the system under consideration.
is to engage in parallel synthesis: A common setting involves a
scaffold with one or several defined exit vectors and a set of
chemical reactions with the goal of optimizing side chains. The
number of suitable reactants accessible (internally or
purchasable) can be very large (hundreds, thousands, or
more). The task of molecular modeling consists then in
prioritizing building blocks with respect to binding affinity in
order to limit the amount of synthesis and experimental testing
required.
A prerequisite for this exercise is the availability of an initial
ligand together with structural information, an experimental
cocrystal structure describing the binding mode to the protein.
We picked human cathepsin L (hCatL), a cysteine protease,
which can be inhibited by ligands with an activated nitrile group
forming a covalent thioimidate adduct with the catalytic Cys25.
Previous SAR and structural studies with aryl nitriles (Figure 1)
INTRODUCTION
■
Free energy calculation approaches, such as free energy
perturbation (FEP), have been around for a long time1−5 but
had only limited impact in the drug discovery process so far.
Likely reasons for their historically restrained use include
lengthy simulation times not practical in fast-paced project
environments combined with overstated accuracy levels based
on small test set retrospective analyses which did not translate
when employed prospectively in real-world systems. FEP has
now taken advantage of improved sampling algorithms6,7 and
force-field quality8 and is profiting from the increased
availability of low-cost parallel computing. Speed and accuracy
appear to have progressed significantly.7,9,10 This has in turn led
to recent accounts of successful industrial applications of FEP
in active drug discovery projects.11−13
Here, we investigate the application of FEP in a typical drug
discovery use case where the goal is to prioritize compounds for
synthesis. One way for therapeutic project teams to further
explore the structure−activity relationship (SAR) of a hit series
Received: December 28, 2016
© XXXX American Chemical Society
A
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry
Article
Figure 1. (a) Triazine nitrile reference compound 1.16 The substituents in the S1, S2, and S3 pockets of cathepsin L are highlighted with gray
dashed-line boxes. (b) Binding mode of compound 1 in hCatl (PDB code 4AXM). The dashed orange circle highlights the region of the S2 pocket.
(c) Fluoropyrimidine nitrile reference compound 2.15 (d) Fluoropyrimidine nitrile scaffold used in this study.
revealed that apart from the covalent bond the main protein−
ligand interactions are made in the S2 and S3 pockets.14−16
While the reactivity of the aryl nitrile warhead and substituents
for the S3 pocket have been explored already to a larger extent,
SAR knowledge in the S2 pocket is very limited for this scaffold.
As one of the exit vectors at the central nitrogen atom is
pointing into the S2 pocket (Figure 1) and is amenable to
diverse chemical modifications, we chose this as model system
for our investigation.
We evaluated the current state-of-the-art of FEP in
prioritizing novel substituents with improved target binding
for the S2 pocket of hCatL. The approach was first validated
using the S2 pocket SAR available from the literature. The
evaluation was then conducted by setting up a realistic drug
discovery scenario: Prioritization was carried out in a
prospective manner, and FEP predictions were delivered in a
constraint time frame that matches project requirements. The
performance of FEP is compared to that of established
computational approaches as well as to the selection of a
medicinal chemist which is used as a benchmark to the
theoretical methods. Given the nature of the hCatL system, we
believe this study also represents the first account of FEP
calculations applied to covalent inhibitors.
The S2 substituent is incorporated in the first of a three-step
synthetic route, calling for primary amines for the reductive
amination of piperonal (see Experimental Section). We
compiled a list of 3325 amines that would lead to
fluoropyrimidine compounds with molecular weight of ≤500
Da and that were available in sufficient quantity to carry out the
synthesis. The challenge consisted for each submitter to pick 10
building blocks that would result in fluoropyrimidine nitrile
derivatives with improved hCatL potency compared to
compound 2.
Each participant had access to the full list of 3325 amines, the
list of corresponding fluoropyrimidine compounds (see
Supporting Information Table S1), all available SAR on the
series (see Supporting Information Table S2), and 1 month of
time.14−16 Building blocks were selected by one of the
following four approaches (full details are available from the
Experimental Section):
(1) A medicinal chemist without knowledge of hCatL
biostructure information browsed through the list of
available building block taking previously known SAR
into account; this corresponds to our null hypothesis,
since no modeling approach was employed (hereafter
referred to as MC).
(2) Experts visually assessed the binding mode of com-
pounds minimized with Moloc17 (MAB force field) and
manually selected building blocks based on assumed
favorable interactions (hereafter MM).
(3) An expert visually assessed the binding mode of
compounds docked with Gold, refined using the MAB
force field,17 rescored with ScorpionScore,18 and further
filtered for highly strained conformations and unfavor-
able interactions (hereafter DOMF).
(4) As a prerequisite, the FEP approach was validated using
the S2 pocket SAR previously published in the literature
(see Supporting Information S3). The selection was
based on FEP binding free energy estimations for a set of
93 ligands which remained after processing the full
reactant set with Glide SP docking, MM-GB/SA scoring,
and by-eye diversity analysis (hereafter FEP).
RESULTS AND DISCUSSION
■
Drug Target, Reference Scaffold, and Building Block
Selection. The goal of this study is to compare different
approaches in their ability to select 10 reactants from a sizable
list of available building blocks for array chemistry, which when
incorporated into the ligand will yield tight binding with the
target protein. Due to the availability of a relevant protein−
ligand crystal structure, a reliable biochemical assay for
measuring binding affinities, and established chemical synthesis,
we chose the drug target hCatL in complex with an aryl nitrile
scaffold as model system (Figure 1). Only few substituents have
been reported for the S2 pocket, most of them on the triazine
nitrile series.16 The reference compound 1 has an inhibition
constant of 13 nM against hCatL, and the binding mode is
known from its cocrystal structure (Figure 1a,b). Since this
compound already exhibits a high affinity for the target, we
decided to base our exploration of the S2 pocket on reference
compound 2 (see Figure 1c) which has a reduced reactivity of
the warhead resulting in a lower hCatL affinity (Ki = 200
nM).15 The scaffold resulting from compound 2 used in this
experiment is shown in Figure 1d.
The four sets of selected building blocks are listed in Table 1.
Three building blocks (3, 13, 20) were selected by more than
one approach.
Ligand Binding Affinities, Physical Properties, and
Comparison of Selection Approaches. All corresponding
fluoropyrimidine nitriles were then synthesized except com-
pound 12 for which the synthesis was not successful (see
B
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry
Article
e
Table 1. List of the Building Blocks Selected by Each of the Four Approaches
a
b
Compound 12 could not be synthesized. Compound 25 is missing a methyl on the nitrogen compared to the intended molecule, still maintaining
c
the positive charge. Compound 29 is an unintended product but was approved in the DOMF approach before compounds were tested. It shows
high similarity to compound 28. The original submission for compound 38 had one more carbon atom in the linker. Three building blocks (3, 13,
20) were prioritized by more than one approach.
d
e
Experimental Section for full details), and their inhibitory
constants (Ki) were determined experimentally in a fluo-
rescence-based assay.14,19 Table 2 provides all details of
calculated and measured properties as well as hCatL Ki values,
ligand efficiencies20 (LE), lipophilic efficiencies21 (LipE) and
similarities to the nearest neighbor (see details in the Methods
Section) in the set of known fluoropyrimidine nitriles. Of the
36 compounds synthesized and tested, 9 achieved better
potency than the reference compound 2. Compound 3 was
picked by three approaches (MC, MM, and FEP), compound
22 was selected by the DOMF approach, and all other
compounds were selected by FEP.
Figure 2 shows scatter plots of the inhibitory constants,
ligand efficiencies, and lipophilic efficiencies grouped by
selection method. The FEP approach was able to pick building
blocks that corresponded to compounds with affinities of <200
nM for 8 out of its 10 choices. The other approaches, MC,
MM, and DOMF, were each only successful in one case.
Looking at a less stringent criterion, the ability to prioritize
active compounds (compounds with submicromolar activity, Ki
≤ 1000 nM, are classified as active, and 13 and 20 with Ki
C
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry
Article
Table 2. Calculated (Molecular Weight, AlogP) and Measured (log D and Kinetic Solubility) Properties, hCatL Inhibitory
Constant, Corresponding Ligand and Lipophilic Efficiencies (LE and LipE, Respectively), and Nearest Neighbor Tanimoto
Similarity (NNT) for the Set of Fluoropyrimidines 3−38
a
compd
selected by
molecular weight
AlogP
log D
kin. sol.
hCatL Ki (nM)
LE
LipE
3.90
NNT
3
MC, MM, FEP
MC
354.4
396.5
326.3
348.3
366.3
422.4
422.4
372.4
396.8
394.4
340.4
376.4
377.4
377.4
384.4
418.4
446.4
441.4
362.4
326.3
382.4
394.4
405.4
453.5
454.5
455.3
401.4
442.4
424.9
462.5
420.4
458.4
392.4
454.9
368.3
426.8
4.0
5.1
3.0
3.9
4.1
4.8
4.8
2.9
4.5
4.4
3.8
4.2
3.2
3.0
3.1
4.7
4.9
4.3
3.9
3.1
3.9
4.4
3.7
2.9
5.5
4.9
4.1
4.8
5.3
3.0
4.6
4.0
3.9
5.6
3.0
4.6
<0.2
<0.2
0.3
12
3020
>5100
1800
952
0.42
0.26
<0.24
0.30
0.30
0.25
<0.20
0.33
<0.21
0.31
0.33
0.33
0.26
0.31
0.27
<0.19
<0.18
0.25
<0.22
0.40
0.31
0.26
<0.195
0.27
<0.17
0.30
0.30
0.29
0.34
0.23
0.33
0.28
0.34
0.25
0.40
0.31
0.733
0.860
0.672
1.000
0.852
0.875
0.875
0.873
0.741
0.717
0.721
0.717
0.692
0.694
0.662
0.719
0.730
0.723
0.768
0.729
0.662
0.770
0.734
0.716
0.778
0.652
0.662
0.614
0.672
0.557
0.629
0.597
0.672
0.638
0.717
0.642
4
2.74
3.38
3.44
3.49
0.41
<1.31
1.86
5
MC
6
MC
<0.2
<0.2
<0.2
<0.2
0.4
7
MC
1.93
8 (cis)
8 (trans)
9
MC
3500
>5100
515
0.70
MC
<−0.46
3.41
MC
3.61
10
11
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
MC
<0.2
<0.2
<0.2
<0.2
4.2
>5100
279
<−0.25
2.15
MC
MM, DOMF
MM
1010
217
2.24
2.47
MM
3.39
3.15
3.46
3860
505
2.23
MM
<0.2
<0.2
<0.2
<0.2
<0.2
<0.2
<0.2
<0.2
<0.2
<0.2
0.4
3.25
MM
2790
>5100
>5100
1020
>5100
77
2.46
MM
<−0.39
<−0.58
1.70
MM
MM, DOMF
MM
<0.42
4.00
DOMF
DOMF
DOMF
DOMF
DOMF
DOMF
DOMF
DOMF
FEP
3.02
671
2.32
5100
>5100
411
0.87
<0.55
3.49
3.49
2.12
2.96
<0.2
<0.2
<0.2
<0.2
<0.2
18
>5100
358
<−1.16
1.50
304
2.44
123
2.08
FEP
27
2.25
FEP
1750
30
2.75
FEP
<0.2
2.89
FEP
91
3.01
FEP
<0.2
<0.2
0.5
77
3.19
FEP
2.32
3.16
1430
25
0.29
FEP
4.59
FEP
<0.2
167
2.19
a
Result is the average of three independent measurements.
values of 1010 nM and 1020 nM, respectively, are included),
MC has a success rate of 40%, MM of 50%, DOMF of 70%, and
FEP of 80%. In total, there were only three compounds 3, 22,
and 37 that had better ligand efficiencies than reference
compound 2 with FEP identifying two of those. This is perhaps
surprising as the ligand efficiency of 2 was not that high to
begin with (0.35), even though improving ligand efficiency was
not one of the initially stated objectives. Lipophilic efficiency
was better than compound 2 for 10 of the new fluoropyr-
imidines. MC, MM, and DOMF each identified two of those,
and FEP identified six of them.
From these results, it is clear that FEP outperformed all other
approaches in prioritizing compounds for this particular system.
The scatterplots in Figure 2 show that this is most pronounced
when inhibitory constants are compared. For LE and LipE, the
difference from the other methods is smaller, reflecting the fact
that submissions using FEP generally involved larger and more
lipophilic substituents. The S2 pocket is quite hydrophobic, and
it is possible that the trend of identifying more lipophilic R-
groups may simply be a function of the specifics of that pocket
rather than a feature of the optimization method. Additional
studies in more polar binding sites are required to address this
point.
All fluoropyrimidine compounds obtained from the complete
set of 3325 building blocks have an average similarity to
previously reported inhibitors of 0.643 (see Experimental
Section for details on the similarity calculation). The average
similarity of the molecules of the MC submission is significantly
higher, 0.820, suggesting that the selected molecules are closely
related to the known inhibitors and more so than could be
expected from a representative selection. The average similarity
of the MM and DOMF molecules is also higher than that of the
whole set with 0.716 and 0.715, respectively. In contrast, the
average similarity of molecules selected by FEP is 0.647, in line
with the average of the whole pool of candidates.
The nine newly synthesized hCatL inhibitors with improved
binding have topologies that fall into two classes, hereafter
termed A and B. Molecules of class A contain a small lipophilic
motif adjacent to a methylene linker (3, 22, 37). This topology
was previously unknown and unreported though ligands 3, 22,
D
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry
Article
Figure 2. Scatter plots displaying (a) the inhibitory constant Ki, (b) ligand efficiency,20 and (c) lipophilic efficiency21 of the fluoropyrimidine nitrile
compounds 3−38 grouped by selection method (MC = pink, MM = orange, DOMF = blue, and FEP = green). The dashed lines indicate the Ki,
ligand efficiency, and lipophilic efficiency of reference compound 2 (see Figure 1b, Ki = 200 nM, LE = 0.352, and LipE = 2.55). The second dashed
line on scatter plot (a) shows the arbitrary threshold of 1 μM beyond which compounds are considered inactive.
Figure 3. (a) Cocrystal structure of hCatL with compound 35 (PDB code 5MQY, resolution 1.13 Å). Short intermolecular contacts with distances of
<4.0 Å between the ligand fragment (cyan) in the S2 pocket and protein residues (green) are shown as dashed blue lines. The covalent bond
between the sulfur atom of Cys25 and the nitrile carbon atom of 35 is shown in magenta. Only water molecules with distances less than 5.0 Å from
the ligand are shown. (b) Superimposition of the cocrystal structure of 35 in hCatL (ligand, cyan; protein, green) with a representative dominant
pose from the FEP trajectory (ligand, magenta; protein, salmon).
and 37 are close neighbors of known hCatL inhibitors
(Supporting Information Table S2) with Tanimoto similarities
of 0.733, 0.729, and 0.717, respectively. The relative difficulty to
identify this class is further put into context by the fact that it
was prioritized by all four approaches. Molecules in class B are
bigger and share a terminal aryl ring system which is connected
to the main scaffold by a flexible propyl or oxyethyl linker (30,
31, 33−35, 38). This class of compounds is also novel but was
only prioritized by the FEP approach; 30, 31, 33−35, 38 have
Tanimoto similarities to known inhibitors of 0.614, 0.672,
0.629, 0.597, 0.672, and 0.642, respectively.
Analysis of the Cocrystal Structure of Fluoropyrimi-
dine 35 in Complex with hCatL. The crystallographic
binding mode of compound 35 in hCatL was determined to
better understand the improved binding of these compounds
(Figure 3a). As expected, the fluoropyrimidine nitrile warhead
is engaged in a covalent bond with the sulfur atom of Cys25
and the piperonyl group sits in the S3 pocket engaging in π−π
stacking interactions with the Gly67-Gly68 dipeptide at its base.
The newly introduced ethoxybenzene substituent forms
dispersive contacts with several hydrophobic side chains
(Leu69, Ala135, Ala214) and the apolar surface of the
backbone amide carbonyl of Asp162. In particular, the phenyl
ring sits in a narrow cleft between residues Leu69 and Ala135.
Several short contacts with distances from 3.4−3.7 Å suggest
that this interaction region in the S2 pocket is important for the
improved binding. It is also noteworthy that no obvious ligand
strain is visible in the bound conformation. The Car−Car−O−C
dihedral angle is close to planar (τ = 13°) and the O−C−C−N
torsional angle adopts a gauche-like conformation (τ = 88°),
which is in line with the preferred angle ranges derived from
small molecule crystal structures.22 Representative dominant
E
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry
Article
Figure 4. (a) Representative dominant pose from the FEP trajectory for compound 3 (ligand, magenta; protein, salmon). (b) Top: Cocrystal
structure of hCatL with compound 1 (PDB code 4AXM). Short intramolecular contacts with distances of <2.8 Å in the ligand (cyan) are highlighted
by red, dashed lines. The angle between the two planes p1 and p2 is 13°. Bottom: Small molecule crystal structures of representative molecules
(green). The angles between p1 and p2 are 56° (8, trans), 89° (CSD code SULDEK), 56° (CSD code YAYKEQ), respectively.
Figure 5. Novel substituents identified for the S2 pocket of hCatL. All fragments are shown for which the following two conditions for the full
molecule are fulfilled: (a) Ki < 0.6 μM, (b) similarity to known fluoropyrimidine nitrile of <0.8. The fragments are colored by the selection method
(MC = red, MM = orange, DOMF = blue, and FEP = green).
poses can be extracted from the molecular dynamics (MD)
snapshots of the FEP trajectories. This was achieved by
inspecting by eye 20 equally spaced frames of the molecular
dynamics and selecting a consensus frame capturing the most
common features. Figure 3b reveals a good agreement between
the dominant FEP pose for compound 35 and its experimental
binding mode (the pose was selected prior to the availability of
the crystallographic data). While the protein−ligand inter-
actions in the S2 pocket are very well reproduced, the main
difference occurs in the orientation of the piperonyl ring system
in the S3 pocket. The second, flipped orientation is also part of
the FEP trajectory but to a smaller extent.
Ligand Conformations and Interactions of Most
Active Molecules. Overall, the largest gain in binding affinity
was achieved by a cyclopentylmethyl substituent in the S2
pocket (3 Ki = 12 nM), which is 17-fold more active than the
reference 2 (Ki = 200 nM). A representative dominant pose of
the FEP trajectory of 3 is shown in Figure 4a. In contrast to the
more extended topologies of class B compounds, the
cyclopentyl group does not reach into the cleft between
Leu69 and Ala135 but seems to optimally fill the front part of
the S2 pocket. Five short dispersive interactions can be
detected between the cyclopentyl group and protein atoms in
this region (Met70, Ala135, backbone of Asp162). Apart from
these favorable interactions, we suspect that some ligand strain
in the bound conformation of reference 2 contributes to the
large activity increase of compound 3. This hypothesis is
motivated by two very short intramolecular distances (d < 2.8
Å) observed in the hCatL cocrystal structure of 1 between
atoms of the cyclohexyl ring and the scaffold (Figure 4b). The
angle between the planes formed by the cyclohexyl substituent
(p1) and the aminotriazine scaffold (p2) is close to planar (α =
13°), which is considerably smaller than those observed for
similar topologies in the unbound conformation. For example,
the small molecule crystal structure of compound 8 (trans) has
α = 56° while other molecules in the CSD have even wider
interplanar angles.
The three compounds (3, 22, 37) that improve LE relative to
the reference all belong to class A occupying only the front part
of the S2 pocket. This suggests that the selected class B
molecules, which extend to the aromatic cleft and beyond into
the back region of the S2 pocket, do not pick up as many
favorable interactions as the best class A molecules and/or
suffer from enthalpic or entropic penalties. The most active
F
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry
Article
class B molecules feature a relatively flexible propyl or ethoxy
linker, which only partly occupies the front part of the S2
pocket. We suggest that in a next design step, the linker
flexibility should be reduced, ideally in combination with
improved protein interactions in the linker region. If successful,
this would very likely increase the LE of class B molecules
above the value of the reference molecule.
linker, for example, by replacement of sp3 by sp2 atoms or
introduction of a substituent, can indeed reduce the propensity
for hydrophobic collapse and its associated energy penalty.
However, this depends strongly on the topology of the linker
and further analyses on the specific linker system, as described
here, would be needed to estimate the actual effect.
Novel hCatL Inhibitor Topologies and Hydrophobic
Collapse. From our library selection exercise, several new
molecular topologies to fill the S2 pocket of hCatL emerged
(Figure 5). While aryl ring systems connected to the scaffold by
a propyl or ethoxy linker are generally more potent compared
to the ethyl linker analogues, the latter still show decent
submicromolar hCatL activity and could provide an alternative
starting point for medicinal chemistry optimization. Another
interesting motif with a relatively high LipE of 3.5 for the full
molecule is a thiadiazole attached to a piperidyl substituent.
After the release of the Ki values for the newly synthesized
molecules, we revisited the results from the DOMF procedure
which consisted of docking with subsequent postprocessing
including minimization, rescoring, and filtering steps (see
details in the Experimental Section). Since docking scores are
notoriously uncorrelated with experimental binding affin-
ities,23−25 it is common practice to consider a rather large
number of top-ranked docking solutions and to prune down to
the final selection based on visual inspection of the suggested
binding modes. Interestingly, we found that compound 33,
which is very potent in the enzymatic assay (Ki = 30 nM), was
ranked at position 1 after postprocessing of the GOLD docking
solution with Moloc17 as well as after postprocessing with the
Scorpion scoring function.18 Although the potential favorable
interactions of the terminal phenyl ring in the aromatic cleft of
the S2 pocket were recognized in the docking pose, the
compound was not selected because of potential hydrophobic
collapse. This was triggered by earlier data mining studies in the
CSD, which showed that phenyl rings connected by an
unsubstituted propyl linker assume a collapsed conformation in
27% of all crystal structures with this motif (see details in the
Methods Section). As the ligand needs to bind to hCatL in an
extended conformation, clearly some energy penalty is involved
in disrupting the collapsed conformation.
We then probed whether the FEP simulations of the
unbound ligands in water show any signs of hydrophobic
collapse (see details in the Methods Section). Taking the same
geometric definition for collapse as in the CSD search, we
found that for compounds 31 and 33, which both have propyl
linkers, the ratio of FEP snapshots with collapse to the total
number of FEP snapshots is 27% and 64%, respectively. For the
inhibitor 35, which has an oxyethyl linker, the ratio from the
FEP trajectory analysis is 14% while the CSD statistics for two
phenyl rings connected by an unsubstituted oxyethyl linker is
9%. While the comparison between CSD statistics and FEP
simulations is not trivial (they are carried out at very different
temperatures and experimental conditions), it is encouraging to
see that the FEP results are in the right ballpark and that the
simulations are able to reproduce the higher propensity of
unsubstituted propyl linkers to engage in hydrophobic collapse
compared to the analogous oxyethyl linkers. The combined
CSD and FEP analyses consistently suggest that an energy
penalty has to be paid to break the collapsed conformations
occurring in the unbound state but that this is more than
compensated by favorable interactions with the protein, mainly
in the aromatic cleft of the S2 pocket. Rigidification of the
CONCLUSION
■
The purpose of our study was to evaluate the ability of FEP and
other methods to prioritize potential hCatL reversible covalent
inhibitors. This experiment was run in a fully prospective
manner with synthesis and Ki determination of novel
compounds and with time constraints relevant for therapeutic
projects. Arguably the most important aspect of our FEP
evaluation is the direct comparison with benchmark methods.
In our case, we chose other computational selection methods
well established in the industry as well as the selection of a
medicinal chemist as references. We believe that similar types of
benchmark comparisons are required to objectively assess the
prediction power of FEP and to learn about its deficiencies.
In our study, FEP successfully picked 8 out of 10 compounds
that were more potent than the reference, outperforming the
different baselines, a medicinal chemist, and two typical
modeling approaches. Interestingly, FEP was also able to
identify one additional novel topology for S2 pocket moieties
by correctly anticipating the extent to which compounds would
be subject to hydrophobic collapse. That this effect was
correctly handled is a clear sign of the soundness and usefulness
of this approach. It also resulted in the most compelling
knowledge gain from this study for that target: Additional
designs based on the propyl- and oxyether-linked aryls would
be an obvious follow-up.
From these results, we conclude that FEP is an attractive
approach to prioritize compounds for synthesis, though our
observation is based on this one single system. It should also be
noted that the S2 binding pocket of hCatL has a strong
lipophilic character and that additional studies are required to
further substantiate if the type of results obtained here is also
observed for more polar active sites.
One key point is that FEP scoring for this system was
validated with retrospective literature data for this target before
it was used prospectively. Without retrospective validation
about the relevance and quality of its predictions, the outcome
of the prospective FEP scoring is much less certain. With good
evidence of its suitability, FEP appears well positioned to add
significant value to industrial drug discovery projects. We also
hope that our study design will further stimulate other
prospective evaluations of FEP in different binding site
environments and for alternative drug discovery use cases.
METHODS SECTION
■
Selection by a Medicinal Chemist (MC). Compounds
were selected primarily on the basis of the existing SAR
(Supporting Information Table 2 and refs 14−16), which
already showed low double digit nM potencies for R =
cyclohexyl and phenyl. Thus, phenyl 6 and p-fluorophenyl 7 as
well as the structurally related p-chlorobenzyl 10 and p-
fluorophenethyl 11 motifs were highlighted for synthesis. In
order to probe whether a five-membered heterocyclic motif
such as an oxazole is an appropriate replacement for the phenyl
group, the structurally more distinct cyclopropyloxazole 12 was
included. In addition, cyclohexyl side chains containing small
G
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry
Article
para substituents such as isopropyl- and CF3-derivatives (4 and
8) as well as the tetrahydrothiopyran 9 were selected. To probe
whether less spacious groups compared to cyclohexyl are
tolerated, the cyclopentylmethyl 3 and methylcylopropyl 5
were picked. It should be noted that physicochemical properties
such as solubility and lipophilicity were not guiding influences
in this first round of compound selection. The whole process
took our in-house expert less than 4 h.
and 10 molecules, which were thought to have multiple
favorable interactions with the protein while assuming a strain-
free binding conformation, were finally selected.
Computation times for the individual steps were ∼13 h for
docking and ∼2.5 h for postprocessing using a single CPU on a
multicore Linux workstation (Intel Xeon E5-2650, 2.2 GHz).
Manual filtering of the solutions required approximately 0.5
days.
Selection Based on FEP Scoring of Covalent Inhibitors
(FEP). Theoretical Framing of FEP Scoring of Covalent
Ligands. Covalent ligands bind to the protein through a two-
step process. In the first step, the ligand binds to the protein
active site due to electrostatic and hydrophobic interactions,
forming a noncovalently bound protein−ligand prereaction
complex. During the second step, a chemical reaction occurs
between the warhead of the ligand and a protein residue in the
binding pocket resulting in the formation of a covalent bond.
Two distinct free energies are relevant to covalent ligand
binding. The first is the binding free energy to form the
prereaction complex, which is defined as the free energy
difference between the protein−ligand prereaction complex and
the separately solvated ligand and protein, denoted by ΔGn1 in
Figure 6. The second is the free energy for covalent bond
Selection Based on Manual Modeling (MM). The MM
selection was performed using the program Moloc and the all-
atom MAB force field.17 For each optimization, the protein
coordinates and the covalent thioimidate were kept fixed. Next,
the new ligands with the various S2 substituents were energy-
minimized. While not all compounds of the library comprising
3325 amines could be modeled manually, an initial preselection
of the ligands was performed based on their chemical structure
and on the knowledge of S2 pocket SARs of the previously
reported triazine nitriles.14,16 In addition, conformationally
flexible linkers to the amino group were limited to ethane-1,2-
diyl and substituents with longer chains were not considered. A
first set of ligands, 3, 13, and 17, were selected to fill the S2
pocket with mainly aliphatic groups. A second set, 14−16 and
20, features aromatic S2 substituents to undergo C−H···π
interactions with Leu69. The substitution of phenyl (in 14) by
pyridyl (in 15 and 16) reduces the predicted log P (clogP) of
the ligand from 4.19 down to 3.04. Compounds 17−19 feature
substituents which penetrate deeper inside the S2 pocket,
compared to the previously reported ligands,14,16 in order to
harvest additional interactions with the protein. Biaryl ether 20
(clogP = 4.54), with appropriate torsional angle,22 reaches
deepest into the S2 pocket enabling the pyridine ring to
intercalate between the side chains of Leu69 and Ala214.
The MM selection took a total of approximately 2 days: 1
day to manually prioritize compounds from the list of available
building blocks and 1 day to finalize the selection by modeling
compounds in Moloc.
Figure 6. Thermodynamic cycle of covalent inhibitor binding.
Selection Based on Docking Followed by Manual
Filtering (DOMF). One approach relied on small molecule
docking into the hCatL crystal structure with subsequent
filtering of the top-ranked solutions. We used GOLD26,27 with
the ChemPLP scoring function to generate binding poses in the
protein structure 4AXM. To ensure proper orientation of the
fluoropyrimidine nitrile warhead around the reactive Cys25, we
employed a scaffold constraint with a weight of 5000 for the
constant substructure of Figure 1d with 3D coordinates
generated from the very similar ligand of complex structure
4AXM.16 The nitrile was modeled with sp hybridization, and
the side chain atoms of Cys25 were removed to avoid large
steric clashes. An additional hydrophobic pharmacophore
constraint (sphere radius of 2.0 Å, score contribution of 10
per atom in region) was placed in the S2 pocket at the center of
the cyclohexyl ring of compound 2 to make sure that this apolar
pocket, which is important for hCatL binding, is adequately
sampled. Additional postprocessing was done by minimizing
the docked solutions with Mol3d using the MAB force field,17
and subsequent energy estimation with the ScorpionScore18
energy function. The 100 top-ranked molecules from each of
the three energy evaluations (ChemPLP, MAB, ScorpionScore)
were combined, and highly constrained ligands with at least one
torsional angle clearly outside the dihedral angle distribution of
CSD small molecule crystal structures (red flag in program
TorsionAnalyzer28,29) were removed. The docking poses of the
remaining compounds were manually inspected using Moloc,17
formation, which is defined as the free energy difference
between the covalently bonded protein−ligand complex and
protein−ligand prereaction complex, denoted by ΔGc1 in
Figure 6. The overall binding free energy, which is the free
energy difference between the covalently bonded protein−
ligand complex and the state when the protein and ligand are
separated in solution, is equal to the summation of the above
two free energies.
In alchemical free energy calculations, the free energy
difference between two states is calculated by gradually
transforming the Hamiltonian from the first state to the second
state. Through alchemical FEP simulations, we can calculate the
free energy to mutate from ligand 1 to ligand 2 (in this context
free energy refers to the Gibbs free energy) in bulk solution
(ΔGs in Figure 6), in the protein binding pocket (ΔGn in
Figure 6), as well as in the covalently bonded protein−ligand
complexes (ΔGc in Figure 6). Through the thermodynamic
cycles shown in Figure 6, the relative binding free energy
difference between two ligands for prereaction complex
formation can be calculated by performing FEP simulations
corresponding to the first two vertical legs of Figure 6, as is
done is in the typical alchemical free energy calculations for
noncovalent ligand binding. Similarly, by performance of FEP
simulations corresponding to the first and the third vertical legs
of Figure 6 (i.e., the difference of ΔGs and ΔGc), the relative
overall binding free energy difference between the two
reversible covalent ligands can also be directly obtained.
H
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry
Article
An example perturbation path for the calculation of the
relative binding free energy between a pair of covalent ligands
binding to hCatL is shown in Figure 7. In the first leg, the free
FEP/REST method reported in a prior publications was used
to enhance the sampling of protein−ligand binding pocket,7,30
and the Bennett acceptance ratio method was used to calculate
the free energy difference between neighboring λ windows.
Docking and Compound Selection for FEP Analysis. The
PDB structure 4AXM was prepared with the Protein
Preparation Wizard program assuming a pH of 7 and used as
the starting structure for all docking analyses. The cocrystallized
ligand was modified to coincide with its prereactive form, and
the catalytic residue Cys25 was mutated to alanine to
accommodate the additional steric bulk of the reactive nitrile.
All 3325 ligand design ideas were docked in the active site of
the modified 4AXM structure using Glide SP core constrained
docking to the modified cocrystallized ligand.31 Of these 3325
virtual ligands, 2258 did not return a valid docking pose, due to
severe steric clashes and other distortions while respecting the
core constraints. The 1067 compounds that were found to dock
successfully were carried forward to MM-GB/SA scoring.32 The
Z-score data fusion method was used to combine the MM-GB/
SA and Glide SP scores for these ligands into a single
composite score.32 A human expert was asked to select a
diverse set of ∼100 molecules for FEP+ scoring. In consultation
with the human expert, 93 molecules were ultimately scored,
and the top scoring of these were recommended for synthesis.
A table detailing the ranks of the synthesized molecules in the
original 1067 Glide SP, MM-GB/SA, and Z-score rank-ordered
lists may be found in Table 3.
Figure 7. Example perturbation pathway for the calculation of relative
binding free energy between a pair of covalent ligands binding to
hCatL. Note, in the ΔGs, that a fictitious dipeptide reference is used to
ensure the interactions assigned to the dummy atoms are strictly
identical in the complex and solvent legs of the FEP simulations.
energy to mutate from the covalently bonded protein−ligand 1
complex to protein−ligand 2 complex is calculated, and in the
second leg, the free energy to mutate from ligand 1 to ligand 2
in bulk solution is calculated. The difference between the two
free energies gives the relative binding free energy of the two
ligands.
When comparing the relative overall binding free energy
calculated in this way to experimental data, we are assuming
that the concentration of the protein−ligand prereaction
complex does not appreciably contribute to the experimentally
measured affinity. In practice, the protein−ligand prereaction
complex might also contribute to the measured affinity;
however, if the equilibrium concentration of the protein−
ligand prereaction complex is much smaller than that of the
covalently bonded complex, the contribution should be
negligible.
It should also be noted that the covalent bond between the
warhead of the ligand and the corresponding protein residue is
modeled by a harmonic potential in classical-mechanics-based
force fields. Therefore, the free energy difference between the
ligands due to the changing reactivity of the different warheads,
which is mainly due to quantum effects, cannot be accurately
modeled in this protocol. Thus, the above method should only
be accurate for the calculation of relative binding free energies
between covalent ligands with the same warhead and when the
mutation is relatively far away from the warhead. This is the
case in the study presented here.
Details of the FEP Simulation. The FEP+ product
implemented in Schrodinger software suite was used for all
the simulations.7 The PDB structure 4AXM was prepared with
the Protein Preparation Wizard program assuming a pH of 7
and used as the starting structure for all simulations. The
OPLS3 force field was used for the protein and ligands, and the
SPC water model was used for solvent.8 The systems were
relaxed by a series of short MD simulations through the
standard relaxation protocol implemented in FEP+ reported in
prior publications.7,30 A total number of 12 λ windows were
used for both the complex and solvent simulations, and the
production simulations lasted 5 ns for each λ window. The
Table 3. Ranks of the Molecules Selected for Synthesis on
the Basis of the FEP+ Scoring in the Glide SP, MM-GB/SA,
a
and Z-Score Rank Ordered Lists
compd
Glide SP rank
MM-GB/SA rank
Z-score rank
3
189
99
62
17
31
8
501
691
303
37
267
265
96
37
35
33
31
30
36
34
32
8
177
488
284
133
94
35
55
28
44
412
56
38
187
a
Molecule 38 differed from the originally submitted molecule by
omission of a CH2 linker carbon, and its rank in the triage process is
thus unavailable.
Computation Time. The docking and MM-GB/SA steps,
including the human time investment, took approximately 2
days using standard computing resources. Determining the
binding free energies of the 93 selected molecules took 7360
GPU hours on K80 equivalent GPUs or approximately 5 days
of exclusive access to a 64 GPU cluster. Human time required
to ensure atom mappings and interaction mappings were
correct for all FEP calculations reported herein was significant.
However, the laborious parts of this step have been fully
automated over the course of completing this work, and only
minimal human time would be required to run similar covalent
inhibitor R-group optimizations in the future.
Similarity Calculation. The similarity between two
molecules was determined by calculating the Tanimoto
value33 between their ECFP4 descriptors.34 The set of known
hCatL fluoropyrimidine nitrile inhibitors was compiled by
grouping the 13 different substituents previously reported (see
I
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry
Article
a
Scheme 1. Synthesis of Fluoropyrimidine Nitriles 3−38
a
DMSO = dimethylsulfoxide, DMF = N,N-dimethylformamide, DABCO = 1,4-diazabicyclo[2.2.2]octane, TEA = triethylamine.
Supporting Information Table S2) and replacing their core by
the fluoropyrimidine nitrile used in this study. The resulting set
of 13 compounds is virtual but correctly defines the known
chemical space (this is why compound 6 has a Tanimoto
similarity of 1.0 in Table 2). The similarity between each
molecule prioritized and the closest known fluoropyrimidine
nitrile was calculated by determining the similarity between that
molecule and all eight reference molecules and keeping the
maximum value. The similarity between sets of molecules
(complete list and MC, MM, DOMF, and FEP selections) was
computed by calculating the maximum similarity between each
molecule in the set, in turn, and the eight reference compounds
and determining the arithmetic mean of those Tanimoto scores.
Investigation of Hydrophobic Collapse. Experimental
propensities to form conformations of hydrophobic collapse
were derived by data mining of the CSD, version 5.33, using the
ConQuest 1.14 program. The following general search flags
were set: R ≤ 0.10, “3D coordinates determined”, “not
disordered”, “no ions”, “no errors”, “not polymeric”, and
“only organic”. The two molecular topologies investigated for
this study were two phenyl rings A and B connected by an
unsubstituted n-propyl and oxyethyl linker, respectively. The
propensity for a given topology to form hydrophobic collapse
was calculated from the ratio of CSD entries showing a
collapsed conformation to the total number of entries. A CSD
crystal structure conformation was considered collapsed if at
least three distances d between any atom of phenyl ring A and
any atom of phenyl ring B fulfilled d ≤ 4.0 Å. We verified this
definition in model systems of parallel-displaced and edge-to-
face arrangements, and this seemed to be more general
compared to definitions using ring centroids.
EXPERIMENTAL SECTION
■
Synthesis of the Fluoropyrimidine Derivatives. The series of
37 5-fluoropyrimidinenitriles 3−38 (except 12) was prepared by a
three-step protocol based on reductive amination, nucleophilic
aromatic substitution, and palladium-catalyzed cyanation (see Scheme
1). The first two steps relied on procedures previously reported for the
synthesis of triazine nitriles.16 Reductive amination of piperonal and
corresponding primary amines gave all the desired secondary amines
41−75 in good to excellent yields (Supporting Information Tables
S4−S6). Nucleophilic substitution of 2,4-dichloro-5-fluoropyrimidine
with secondary amines in the presence of Hunig base in i-PrOH/
EtOH 2:1 furnished 2-chloropyrimidin-4-amines 76−109 with yields
in broad range 23−97% (Supporting Information Table S5).
̈
Final nitriles were synthesized by palladium-catalyzed cyanation
using Buchwald third generation Pd-precatalyst t-BuXPhos and
Zn(CN)2 as a cyanide source in a mixture of THF and water 1:5 at
50 °C.35 In the case of ligands containing a second chlorine atom in
the molecule (10, 31, 36, 38), only 0.5 equiv of Zn(CN)2 was used to
avoid dicyanation. DMF was employed as a solvent for the synthesis of
ligands 4 and, 9, which are insoluble in the THF/water mixture.
Palladium-catalyzed cyanation failed in the cases of ligands 28 and,
34. These two ligands were instead prepared by nucleophilic
substitution with KCN in the presence of DABCO in a DMSO/
water mixture. Reactions had to be performed at lower temperature
(60 °C) to avoid decomposition and for longer reaction times (3 days)
than in the published procedure15 to get ligands 28 and, 34 in low
yields (14% and 43%, respectively). Detailed synthesis, yields,
cyanation optimization, and characterization of all compounds are
available in the Supporting Information.
Chemical Synthesis. All reagents were purchased from
commercial suppliers and used without further purification. Primary
amines for the synthesis were provided by F. Hoffmann-La Roche Ltd.
and were used as received. Solvents for the synthesis were purchased
in HPLC quality, and solvents for extractions and chromatography
were of technical grade and were distilled prior to use. Thin layer
chromatography (TLC) was conducted on aluminum sheets coated
with SiO2 60 F254 (Merck), visualization with UV lamp (254 nm) or by
ninhydrin staining (1.5 g in 100 mL of n-butanol and 3 mL of glacial
acetic acid). MPLC was performed on an ISCO Combiflash Rf 200
system on RediSep Rf normal-phase silica flash columns (particle size
35−70 μm, 230−400 mesh). HPLC analysis was performed on a
Merck Hitachi LaChrom HPLC system (D-7000 interface, L-7100
pump, L-7200 autosampler, and L-7400 UV detector). Separations
under acidic conditions (HPLC method A) were performed using
FEP simulations of the solvent leg were analyzed for
hydrophobic collapse using the same geometric criteria as in
the CSD searches. A total of 22 snapshots were taken into
account for each ligand investigated. These snapshots were
taken from 20 equally spaced frames of the 5 ns solvent leg
production stage FEP simulation and included data from λ0 and
λ11 windows where a single ligand species is simulated in each
window.
J
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry
Article
0.1% vol formic acid in eluent H2O and MeCN on a Merck
Superspher 100 RP-18e (100 Å, 4 μm), 250 mm × 4 mm column.
Separations under basic conditions (HPLC method B) were
performed using a 10 mM K2HPO4 buffer (pH 8) in eluent H2O
and MeCN/H2O (7:3) on a SGE ProteCol C18H (120 Å, 5 μm) 250
mm × 4.6 mm column. A flow rate of 1 mL/min was used, and the UV
detector was set to 254 nm (optical path length of 10 mm). Melting
equation36 from the IC50 determined in an enzymatic assay following
the same procedure as previously described.15
Cocrystallization of hCatL in Complex with Ligand 35 (PDB
code 5MQY). Protein at a concentration of 6.7 μM in 100 mM
sodium acetate, pH 5.5, 5 mM DTT, 5 mM EDTA, 0.02% NaN3
(buffer A) was incubated with 35 in a 12-fold molar excess overnight at
4 °C under argon. Subsequently the protein−ligand solution was
diluted by addition of 100 mL of buffer A and further incubated for 2 h
at 21 °C. Precipitated material was removed by centrifugation. Prior to
crystallization experiments, the protein−ligand mixture was concen-
trated to 30.4 mg/mL and centrifuged at 20 000g. The crystallization
droplets were set up at 21 °C by mixing 0.1 μL of protein solution
with 0.17 μL of reservoir solution and 0.03 μL of seed solution of CatL
crystals in sitting drop vapor diffusion experiments. Crystals were
obtained within 1 day out of 25% PEG 3350, 0.1 M Bis-Tris, pH 6.5.
points (mp) were measured on a Buchi M-560 melting point apparatus
̈
in open capillaries and are uncorrected. Infrared spectra (IR) were
recorded on a PerkinElmer Spectrum Two FT-IR device fitted with an
ATR unit (4000−600 cm−1). Absorption bands are given in
wavenumbers (cm−1), and signal intensities are specified with s
(strong), m (medium), and w (weak). NMR spectra were measured
on a Varian Gemini 300 MHz spectrometer (300 MHz for 1H and 75
MHz for 13C) or a Bruker DRX 400 MHz spectrometer (400 MHz for
1H, 100 MHz for 13C, and 376 MHz for 19F) in CDCl3 (referenced to
ASSOCIATED CONTENT
the residual solvent signal, CDCl3, δ H = 7.26, δ C = 77.16 ppm) at 25
°C. Chemical shifts are given in ppm (δ scale), coupling constants (J)
in Hz. Complete assignment of all NMR signals was performed using a
combination of H,H-COSY, H,C-HSQC, and H,C-HMBC experi-
ments. High-resolution electrospray ionization mass spectrometry
(HR-ESI-MS) was conducted with a Bruker Daltonics maXis ESI/
Nanospray Qq TOF instrument. All mass spectra were acquired by the
MS service of the Laboratory of Organic Chemistry at ETH Zurich.
Elemental analyses were performed by the Mikrolabor of the
Laboratory of Organic Chemistry at ETH Zurich with a LECO
CHNS/932 instrument. Purity of all final compounds (>95%) was
determined by elemental analyses or by analytical HPLC. The
nomenclature was obtained with ACD/Name software (ACD/Labs).
General Procedure A: Reductive Amination. A solution of the
primary amine (1.0 equiv) and piperonal (39) (0.9 equiv) in CH2Cl2
(3.5 mL/1 mmol of amine) over molecular sieves (4 Å) (if starting
material was a hydrochloride salt, triethylamine (1.0 equiv) was added)
was stirred at 25 °C for 10−18 h, treated with NaBH(OAc)3 (2.5
equiv), and stirred at 25 °C for another 10−18 h. The mixture was
diluted with CH2Cl2 and sat. NaHCO3 solution and filtered. The
organic layer was washed with sat. NaHCO3 solution (3×), dried over
Na2SO4, filtered, and evaporated. The crude material was purified by
MPLC.
■
S
* Supporting Information
. The Supporting Information is available free of charge on the
ACS Publications website at DOI: 10.1021/acs.jmed-
chem.6b01881.
Molecular strings and some data (CSV)
List of all virtual molecules to prioritize, previsously
known SAR for the hCatL triazine series, retrospective
validation of FEP+ approach, chemical diagrams of CSD
structures SULDEK and YAYKEQ, compound character-
ization, yields, cyanation optimization, small molecule X-
ray crystallography (PDF)
Accession Codes
Atomic coordinates and experimental data for the cocrystal
structure of 35 (PDB code 5MQY) in complex with hCatL will
be released upon article publication.
AUTHOR INFORMATION
■
Corresponding Authors
*R.A.: phone, +1 212 548 2352; e-mail, robert.abel@
schrodinger.com.
General Procedure B: Nucleophilic Aromatic Substitution.
The secondary amine (1.0 equiv) and 2,4-dichloro-5-fluoropyrimidine
(40) (1.0 equiv) were dissolved in i-PrOH or i-PrOH/EtOH (2:1)
(5.0 mL/1 mmol of amine), i-Pr2NEt (1.2 equiv) was added, and the
mixture was heated at 85 °C for 18−24 h. Solvents were evaporated,
and the crude material was purified by MPLC.
*F.D.: phone, +41 44 63 22992; e-mail, diederich@org.chem.
ethz.ch.
*J.H.: phone, +41 61 68 89773; e-mail, jerome.hert@roche.
com.
General Procedure C: Palladium-Catalyzed Cyanation in
THF/H2O .35 A Biotage microvawe reaction vial (2−5 mL) was
charged with a 2-chloro-5-fluoro-4-pyrimidinamine (1.0 equiv),
Zn(CN)2 (0.66 equiv or 0.5 equiv if the pyrimidinamine contained
an additional chlorine atom), and t-BuXPhos G3 Pd-precatalyst (0.05
equiv). The vial was capped with septum, evacuated, and refilled with
argon. THF/degassed H2O 1:5 (3 mL/0.2 mmol of pyrimidinamine)
was added, and the mixture was vigorously stirred (1200 rpm) at 60
°C for 12−18 h. Efficient stirring is crucial for complete conversion.
The mixture was diluted with sat. NaHCO3 solution (2 mL) and
EtOAc (2 mL) and stirred for 5 min. The layers were separated, and
the aqueous layer was extracted with EtOAc (3 × 2 mL). The
combined organic layers were dried over Na2SO4, filtered, and
evaporated. The crude material was purified by MPLC.
ORCID
Michal Tichy: 0000-0002-8460-7680
Lingle Wang: 0000-0002-8170-6798
Rainer E. Martin: 0000-0001-7895-497X
́
Jerome Hert: 0000-0003-3062-8029
́
̂
Author Contributions
#B.K., M.T., and L.W. contributed equally. The manuscript was
written through contributions of all authors. All authors have
given approval to the final version of the manuscript.
Notes
The authors declare no competing financial interest.
General Procedure D: Palladium-Catalyzed Cyanation in
DMF. A Biotage microwave reaction vial (2−5 mL) was charged with a
2-chloro-5-fluoro-4-pyrimidinamine (1.0 equiv), Zn(CN)2 (0.66
equiv), and t-BuXPhos G3 Pd-precatalyst (0.05 equiv). The vial was
capped with septum, evacuated, and refilled with argon. DMF (4 mL/
0.5 mmol of pyrimidinamine) was added, and the mixture was stirred
at 60 °C for 1−12 h. The solvent was coevaporated with toluene, and
the crude product was purified by MPLC.
ACKNOWLEDGMENTS
■
We thank Dr. Bruno Bernet for help with the compound
characterizations, Dr. Nils Trapp for small molecule crystallo-
graphy, and Dr. Julian Fuchs for help in deriving CSD statistics
for the investigation of hydrophobic collapse. Work at ETH
Zurich was supported by a grant from the ETH Research
Council (Grant ETH-01 13-2) and by an F. Hoffmann-La
Roche−ETH Zurich research agreement.
Determination of hCatL Activity. The inhibitory constants (Ki)
of the ligands against hCatL were calculated using the Cheng−Prusoff
K
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry
Article
(14) Ehmke, V.; Heindl, C.; Rottmann, M.; Freymond, C.;
Schweizer, W. B.; Brun, R.; Stich, A.; Schirmeister, T.; Diederich, F.
Potent and Selective Inhibition of Cysteine Proteases from
Plasmodium Falciparum and Trypanosoma Brucei. ChemMedChem
2011, 6, 273−278.
ABBREVIATIONS USED
■
CSD, Cambridge Structural Database; DOMF, docking and
manual filtering; ECFP, extended-connectivity fingerprints;
FEP, free energy perturbation; hCatL, human cathepsin L;
HR-ESI-MS, high-resolution electrospray ionization mass
spectrometry; LipE, lipophilic efficiency; MC, medicinal
chemist; MM, manual modeling; MM-GB/SA, molecular
mechanics generalized Born/surface area; MPLC, medium
pressure liquid chromatography
(15) Ehmke, V.; Quinsaat, J. E. Q.; Rivera-Fuentes, P.; Heindl, C.;
Freymond, C.; Rottmann, M.; Brun, R.; Schirmeister, T.; Diederich, F.
Tuning and Predicting Biological Affinity: Aryl Nitriles as Cysteine
Protease Inhibitors. Org. Biomol. Chem. 2012, 10, 5764−5768.
(16) Ehmke, V.; Winkler, E.; Banner, D. W.; Haap, W.; Schweizer, W.
B.; Rottmann, M.; Kaiser, M.; Freymond, C.; Schirmeister, T.;
Diederich, F. Optimization of Triazine Nitriles as Rhodesain
Inhibitors: Structure-Activity Relationships, Bioisosteric Imidazopyr-
idine Nitriles, and X-Ray Crystal Structure Analysis with Human
Cathepsin L. ChemMedChem 2013, 8, 967−975.
REFERENCES
■
(1) Postma, J. P. M.; Berendsen, H. J. C.; Haak, J. R.
Thermodynamics of Cavity Formation in Water. A Molecular
Dynamics Study. Faraday Symp. Chem. Soc. 1982, 17, 55−67.
(2) Tembre, B. L.; Mc Cammon, J. A. Ligand-Receptor Interactions.
Comput. Chem. 1984, 8, 281−283.
(3) Jorgensen, W. L.; Ravimohan, C. Monte Carlo Simulation of
Differences in Free Energies of Hydration. J. Chem. Phys. 1985, 83,
3050−3054.
(4) Kollman, P. Free Energy Calculations: Applications to Chemical
and Biochemical Phenomena. Chem. Rev. 1993, 93, 2395−2417.
(5) Kollman, P. A.; Massova, I.; Reyes, C.; Kuhn, B.; Huo, S.; Chong,
L.; Lee, M.; Lee, T.; Duan, Y.; Wang, W.; Donini, O.; Cieplak, P.;
Srinivasan, J.; Case, D. A.; Cheatham, T. E. III. Calculating Structures
and Free Energies of Complex Molecules: Combining Molecular
Mechanics and Continuum Models. Acc. Chem. Res. 2000, 33, 889−
897.
(17) Gerber, P. R.; Muller, K. MAB, a Generally Applicable
̈
Molecular Force Field for Structure Modelling in Medicinal
Chemistry. J. Comput.-Aided Mol. Des. 1995, 9, 251−268.
(18) Kuhn, B.; Fuchs, J. E.; Reutlinger, M.; Stahl, M.; Taylor, N. R.
Rationalizing Tight Ligand Binding through Cooperative Interaction
Networks. J. Chem. Inf. Model. 2011, 51, 3180−3198.
(19) Vicik, R.; Busemann, M.; Baumann, K.; Schirmeister, T.
Inhibitors of Cysteine Proteases. Curr. Top. Med. Chem. 2006, 6, 331−
353.
(20) Hopkins, A. L.; Groom, C. R.; Alex, A. Ligand Efficiency: A
Useful Metric for Lead Selection. Drug Discovery Today 2004, 9, 430−
431.
(21) Leeson, P. D.; Springthorpe, B. The Influence of Drug-like
Concepts on Decision-Making in Medicinal Chemistry. Nat. Rev. Drug
Discovery 2007, 6, 881−890.
(22) Brameld, K. A.; Kuhn, B.; Reuter, D. C.; Stahl, M. Small
Molecule Conformational Preferences Derived from Crystal Structure
Data. A Medicinal Chemistry Focused Analysis. J. Chem. Inf. Model.
2008, 48, 1−24.
(6) Gallicchio, E.; Levy, R. M. Advances in All Atom Sampling
Methods for Modeling Protein−Ligand Binding Affinities. Curr. Opin.
Struct. Biol. 2011, 21, 161−166.
(7) Wang, L.; Wu, Y.; Deng, Y.; Kim, B.; Pierce, L.; Krilov, G.;
Lupyan, D.; Robinson, S.; Dahlgren, M. K.; Greenwood, J.; Romero,
D. L.; Masse, C.; Knight, J. L.; Steinbrecher, T.; Beuming, T.; Damm,
W.; Harder, E.; Sherman, W.; Brewer, M.; Wester, R.; Murcko, M.;
Frye, L.; Farid, R.; Lin, T.; Mobley, D. L.; Jorgensen, W. L.; Berne, B.
J.; Friesner, R. A.; Abel, R. Accurate and Reliable Prediction of Relative
Ligand Binding Potency in Prospective Drug Discovery by Way of a
Modern Free-Energy Calculation Protocol and Force Field. J. Am.
Chem. Soc. 2015, 137, 2695−2703.
(8) Harder, E.; Damm, W.; Maple, J.; Wu, C.; Reboul, M.; Xiang, J.
Y.; Wang, L.; Lupyan, D.; Dahlgren, M. K.; Knight, J. L.; Kaus, J. W.;
Cerutti, D. S.; Krilov, G.; Jorgensen, W. L.; Abel, R.; Friesner, R. A.
PLS3: A Force Field Providing Broad Coverage of Drug-like Small
Molecules and Proteins. J. Chem. Theory Comput. 2016, 12, 281−296.
(9) Homeyer, N.; Stoll, F.; Hillisch, A.; Gohlke, H. Binding Free
Energy Calculations for Lead Optimization: Assessment of Their
Accuracy in an Industrial Drug Design Context. J. Chem. Theory
Comput. 2014, 10, 3331−3344.
(23) Perola, E.; Walters, W. P.; Charifson, P. S. A Detailed
Comparison of Current Docking and Scoring Methods on Systems of
Pharmaceutical Relevance. Proteins: Struct., Funct., Genet. 2004, 56,
235−249.
(24) Plewczynski, D.; Łazniewski, M.; Augustyniak, R.; Ginalski, K.
́
Can We Trust Docking Results? Evaluation of Seven Commonly Used
Programs on PDBbind Database. J. Comput. Chem. 2011, 32, 742−
755.
(25) Gathiaka, S.; Liu, S.; Chiu, M.; Yang, H.; Stuckey, J. A.; Kang, Y.
N.; Delproposto, J.; Kubish, G.; Dunbar, J. B., Jr; Carlson, H. A.;
Burley, S. K.; Walters, W. P.; Amaro, R. E.; Feher, V. A.; Gilson, M. K.
D3R Grand Challenge 2015: Evaluation of Protein−Ligand Pose and
Affinity Predictions. J. Comput.-Aided Mol. Des. 2016, 30, 651−668.
(26) Jones, G.; Willett, P.; Glen, R. C. Molecular Recognition of
Receptor Sites Using a Genetic Algorithm with a Description of
Desolvation. J. Mol. Biol. 1995, 245, 43−53.
(10) Steinbrecher, T. B.; Dahlgren, M.; Cappel, D.; Lin, T.; Wang, L.;
Krilov, G.; Abel, R.; Friesner, R.; Sherman, W. Accurate Binding Free
Energy Predictions in Fragment Optimization. J. Chem. Inf. Model.
2015, 55, 2411−2420.
(27) Jones, G.; Willett, P.; Glen, R. C.; Leach, A. R.; Taylor, R.
Development and Validation of a Genetic Algorithm for Flexible
Docking. J. Mol. Biol. 1997, 267, 727−748.
(28) Scharfer, C.; Schulz-Gasch, T.; Ehrlich, H.-C.; Guba, W.; Rarey,
̈
(11) Ciordia, M.; Per
́
ez-Benito, L.; Delgado, F.; Trabanco, A. A.;
M.; Stahl, M. Torsion Angle Preferences in Druglike Chemical Space:
A Comprehensive Guide. J. Med. Chem. 2013, 56, 2016−2028.
(29) Guba, W.; Meyder, A.; Rarey, M.; Hert, J. Torsion Library
Reloaded: A New Version of Expert-Derived SMARTS Rules for
Assessing Conformations of Small Molecules. J. Chem. Inf. Model.
2016, 56, 1−5.
(30) Wang, L.; Deng, Y.; Knight, J. L.; Wu, Y.; Kim, B.; Sherman, W.;
Shelley, J. C.; Lin, T.; Abel, R. Modeling Local Structural
Rearrangements Using FEP/REST: Application to Relative Binding
Affinity Predictions of CDK2 Inhibitors. J. Chem. Theory Comput.
2013, 9, 1282−1293.
(31) Friesner, R. A.; Banks, J. L.; Murphy, R. B.; Halgren, T. A.;
Klicic, J. J.; Mainz, D. T.; Repasky, M. P.; Knoll, E. H.; Shelley, M.;
Perry, J. K.; Shaw, D. E.; Francis, P.; Shenkin, P. S. Glide: A New
Tresadern, G. Application of Free Energy Perturbation for the Design
of BACE1 Inhibitors. J. Chem. Inf. Model. 2016, 56, 1856−1871.
(12) Lovering, F.; Aevazelis, C.; Chang, J.; Dehnhardt, C.; Fitz, L.;
Han, S.; Janz, K.; Lee, J.; Kaila, N.; McDonald, J.; Moore, W.; Moretto,
A.; Papaioannou, N.; Richard, D.; Ryan, M. S.; Wan, Z.-K.;
Thorarensen, A. Imidazotriazines: Spleen Tyrosine Kinase (Syk)
Inhibitors Identified by Free-Energy Perturbation (FEP). ChemMed-
Chem 2016, 11, 217−233.
(13) Abel, R.; Mondal, S.; Masse, C.; Greenwood, J.; Harriman, G.;
Ashwell, M. A.; Bhat, S.; Wester, R.; Frye, L.; Kapeller, R.; Friesner, R.
A. Accelerating Drug Discovery through Tight Integration of Expert
Molecular Design and Predictive Scoring. Curr. Opin. Struct. Biol.
2017, 43, 38−44.
L
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX

Journal of Medicinal Chemistry
Article
Approach for Rapid, Accurate Docking and Scoring. 1. Method and
Assessment of Docking Accuracy. J. Med. Chem. 2004, 47, 1739−1749.
(32) Knight, J. L.; Krilov, G.; Borrelli, K. W.; Williams, J.; Gunn, J. R.;
Clowes, A.; Cheng, L.; Friesner, R. A.; Abel, R. Leveraging Data
Fusion Strategies in Multireceptor Lead Optimization MM/GBSA
End-Point Methods. J. Chem. Theory Comput. 2014, 10, 3207−3220.
(33) Willett, P.; Barnard, J. M.; Downs, G. M. Chemical Similarity
Searching. J. Chem. Inf. Comput. Sci. 1998, 38, 983−996.
(34) Rogers, D.; Hahn, M. Extended-Connectivity Fingerprints. J.
Chem. Inf. Model. 2010, 50, 742−754.
(35) Cohen, D. T.; Buchwald, S. L. Mild Palladium-Catalyzed
Cyanation of (Hetero)aryl Halides and Triflates in Aqueous Media.
Org. Lett. 2015, 17, 202−205.
(36) Cheng, Y.-C.; Prusoff, W. H. Relationship between the
Inhibition Constant (KI) and the Concentration of Inhibitor Which
Causes 50 per Cent Inhibition (I50) of an Enzymatic Reaction.
Biochem. Pharmacol. 1973, 22, 3099−3108.
M
DOI: 10.1021/acs.jmedchem.6b01881
J. Med. Chem. XXXX, XXX, XXX−XXX