Molecular Dynamics Simulations and Kinetic Measurements to Estimate and Predict Protein-Ligand Residence Times

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

Ligand-target residence time is emerging as a key drug discovery parameter because it can reliably predict drug efficacy in vivo. Experimental approaches to binding and unbinding kinetics are nowadays available, but we still lack reliable computational to

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

Subscriber access provided by La Trobe University Library
Article
Molecular dynamics simulations and kinetic measurements
to estimate and predict protein-ligand residence times
luca mollica, Isabelle Theret, Mathias Antoine, Francoise Perron-Sierra, Yves
Charton, Jean-Marie Fourquez, Michel Wierzbicki, Jean A. Boutin, Gilles Ferry,
Sergio Decherchi, Giovanni Bottegoni, Pierre Ducrot, and Andrea Cavalli
J. Med. Chem., Just Accepted Manuscript • DOI: 10.1021/acs.jmedchem.6b00632 • Publication Date (Web): 08 Jul 2016
Downloaded from http://pubs.acs.org on July 11, 2016
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 free 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 accessible to all
readers and 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.
Journal of Medicinal Chemistry 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 35
Journal of Medicinal Chemistry
1
2
3
4
5
6
Molecular dynamics simulations and kinetic measurements to estimate and
predict protein-ligand residence times
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
Luca Mollica1, Isabelle Theret2, Mathias Antoine2, Françoise Perron-Sierra2, Yves Charton2, Jean-
Marie Fourquez2, Michel Wierzbicki2, Jean A. Boutin2, Gilles Ferry2, Sergio Decherchi1,3, Giovanni
Bottegoni1,3, Pierre Ducrot*2, Andrea Cavalli*1,4
1CompuNet, Istituto Italiano di Tecnologia, via Morego, 30, 16163 Genova, Italy
2Institut de Recherches Servier, 125 Chemin de Ronde, 78290 Croissy-sur-Seine, France
3BiKi Technologies S.r.l. via XX Settembre 33/10 16121 Genova, Italy
4Dept. of Pharmacy and Biotechnology, University of Bologna, via Belmeloro 6, 40126 Bologna,
Italy
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 2 of 35
1
2
3
Abstract
4
5
6
7
8
Ligandꢀtarget residence time is emerging as a key drug discovery parameter because it can reliably
predict drug efficacy in vivo. Experimental approaches to binding and unbinding kinetics are
nowadays available, but we still lack reliable computational tools for predicting kinetics and
residence time. Most attempts have been based on bruteꢀforce molecular dynamics (MD) simulations,
which are CPUꢀdemanding and not yet particularly accurate. We recently reported a new scaledꢀMDꢀ
based protocol, which showed potential for residence time prediction in drug discovery. Here, we
further challenged our procedure’s predictive ability by applying our methodology to a series of
glucokinase inhibitors that could be useful for treating type 2 diabetes mellitus. We combined scaled
MD with experimental kinetics measurements and Xꢀray crystallography, promptly checking the
protocol’s reliability by directly comparing computational predictions and experimental measures.
The good agreement highlights the potential of our scaledꢀMDꢀbased approach as an innovative
method for computationally estimating and predicting drug residence times.
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
ACS Paragon Plus Environment

Page 3 of 35
Journal of Medicinal Chemistry
1
2
3
Introduction
4
5
6
7
8
Binding and unbinding kinetics are emerging as key factors for predicting drug efficacy in living
organisms.1,2 There are several experimental techniques (e.g. SPR, stoppedꢀflow CD, fluorescence
spectroscopy) for studying (un)binding kinetics, but currently no efficient computational approaches
for predicting absolute kinetic parameters. The few attempts reported in the literature have mainly
been based on bruteꢀforce molecular dynamics (MD) simulations,3 which are highly demanding in
terms of time and computational power.4ꢀ6 These methods are also unsuitable for industrial use,
where dozens of compounds must be prioritized in the hitꢀtoꢀlead and the lead optimization phases.3
Importantly, (un)binding rates cannot be directly computed for pharmacologically relevant systems,
since the residence time (tr) of molecules can be in the order of seconds, minutes, or even hours.
Thus, smarter algorithms and effective practical solutions are needed to tackle the problem of kinetic
rate estimation.
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
Researchers recently reported a detailed computational study of proteinꢀligand dissociation based on
extensive metadynamics simulations.7 This work demonstrated the possibility of studying the
mechanisms governing unbinding events, and of disclosing their pathways, rates, and rateꢀlimiting
steps. Although this approach has provided remarkable atomistic information, its practical
effectiveness is limited by the substantial computational resources required (i.e. many weeks on a
highly parallel computational infrastructure) to evaluate each single binding and unbinding kinetic
constant pair (i.e. kon and koff values). Moreover, the setup of metadynamics simulations is not trivial
and requires a strong expertise by the user. The high error in the koff estimation also prevents an easy
application of this approach within the drug discovery pipeline, especially referring to an automated
process. A possible alternative is to combine the kon obtained from unbiased simulations with the
binding free energy estimated using free energy methods. Although promising, this alternative
method is still too inaccurate and computationally demanding for any screening purpose.5
We recently reported a novel computational method that addresses the challenge of predicting
unbinding kinetics, particularly in the hitꢀtoꢀlead and lead optimization phases of the drug discovery
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 4 of 35
1
2
3
4
5
6
7
8
9
process.8 Rather than trying to predict the absolute offꢀrate value, koff ꢁ=ꢁtr−1, on individual complexes,
we proposed a procedure to identify the correct koffꢀbased ranking relationship among congeneric
compounds that bind to a given target, using (possibly) limited computational resources. Our solution
is rooted in enhancing the transition probability between different free energy minima during MD
simulations by smoothing the potential energy surface (PES) of the system using a scaling factor λ
(which ranges from 0 to 1) applied to all the potential energies computed during MD simulations.
This procedure flattens the PES associate to biomolecular events and samples a greater amount of
conformational space in simulation times that are much shorter than standard MD simulations.
The underlying rationale is that simulating a proteinꢀligand complex under scaled potential energy
conditions facilitates the rupture of the key physical interactions that confer stability to the complex,
leading to unbinding in much shorter simulation times. These transitions are considered in a
statistical framework that combines multiple replicas, a regressive predictive model, and statistical
analysis to establish the confidence of the predictions. The averaged residence times can and have
been successfully compared to the experimental ones obtained by inverting the koff values. The details
of the method, the analytical treatment of the residence times, and their comparison with
experimental data are reported in the Supporting Information 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
The proposed method is easy to implement and sufficiently accurate in reproducing experimental
data.8 However, several questions remain regarding its general applicability, especially for series of
chemically unrelated compounds. To clarify these aspects and to assess the validity of our
methodology in a predictive and/or prospective manner, we here focused on glucokinase (GK, GK1,
or hexokinase IV), a system of primary interest for both national healthcare systems and the
pharmaceutical industry. A member of the hexokinase (HK) family, GK9 catalyzes the
phosphorylation of glucose to glucose 6ꢀphosphate, the rateꢀlimiting step in glycolysis. The
pharmacological relevance of GK1 is related to its involvement in type 2 diabetes mellitus (T2DM), a
rapidly expanding public health issue affecting over 150 million people worldwide.10,11 GK plays a
critical role in glucose homeostasis, serving as a glucose sensor for glucoseꢀdependent insulin
ACS Paragon Plus Environment

Page 5 of 35
Journal of Medicinal Chemistry
1
2
3
4
5
6
7
8
9
secretion in pancreatic βꢀcells and regulating glucose uptake and glycogen synthesis in the liver.12
GK1 displays several features which are unique among the hexokinase family members: i) as a
monomeric enzyme, GK achieves positive substrate cooperativity solely through equilibration
between multiple conformations,13,14 ii) GK has a unique biochemical kinetics that accounts for its
role as a glucose sensor, iii) GK has a low affinity for glucose as a substrate with a Ka of ~7 mM,
which is within the physiological glucose range, and iv) the enzyme also has no direct feedback
inhibition from its product glucose 6ꢀphosphate. Together, these kinetic properties enable GK to be
highly responsive to fluctuations in blood glucose levels and to ensure that the glucose metabolic flux
is closely tied to glucose concentration. While several classes of diabetic therapies are available for
clinical use, there is still a significant need for new therapies with improved efficacy. A promising
area of current research involves the use of small molecule allosteric activators of GK (GKAs) to
lower blood glucose and normalize insulin secretion. Therapeutically, it seems that activating GK1 in
the liver and pancreas would be an effective strategy for lowering blood glucose by upregulating
hepatic glucose utilization, downregulating hepatic glucose output, and normalizing glucoseꢀ
stimulated insulin secretion.15,16
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 present work, we investigate the unbinding of a novel series of promising GK modulators at
the atomic level, using computational simulations coupled to experimental kinetic measurements and
Xꢀray crystallography. We were able to correctly rank most of the ligands according to their
residence time and to depict some structureꢀkinetics relationships (SKRs). In parallel, we challenge
our approach, highlighting its major advantages and limitations. Finally, we provide some
perspectives on the methodology in terms of its potential to become the computational method of
choice for predicting drug residence times.
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 6 of 35
1
2
3
Results and Discussion
4
5
6
7
8
The main hypothesis under which the original method8 is able to correctly predict residence time
distributions is related to chemical similarity, namely the presence throughout a series of only minor
chemical modifications (e.g. small functional group replacement relative to the parent compound).
The overall protocol comprises the following steps: i) an initial model for each proteinꢀligand
complex is built. Ideally, one proteinꢀligand complex has been experimentally determined by Xꢀray
crystallography, and chemical replacement is used to generate threeꢀdimensional models of
congeneric compounds bound at the target binding pocket. ii) Multiple replicas of scaled molecular
dynamics (MD) of the partially restrained system are performed and stopped when the ligand is
unbound. iii) The average simulated and experimental unbinding times are normalized with respect to
one randomly selected ligand (in this case the ligand for which the experimental structure is
available) in order to convert the experimental residence times ratio to a scaled ratio via an
exponential relation. iv) A bootstrap analysis on the simulated unbinding times per target is carried
out to assess the statistical significance of the observations and to establish whether an increased
number of replicas is needed.
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 present study, two different issues were addressed: i) determination of computational residence
times for all the available ligands for GK1, their ranking and comparison with the experimental data,
and ii) identification of the physicochemical determinants of residence time variations within the
series for the ligands and protein (i.e. depiction of the SKRs).
Measurement of GKA binding kinetics. Following our initial work,17 we exploited the intrinsic
fluorescence properties of 1 (RO4389620)17 (Figure 1) to use it as a fluorescent reporter in a
competitive displacement kinetic experiment to rapidly assess the residence time (1/koff) of several
GKAs, in the presence of glucose to lock GK in its closed conformation. The observed dissociation
transients were monophasic, suggesting a singleꢀstep dissociation. The GKAs’ association kinetics
was also measured, using the quenching of the GK intrinsic fluorescence as readout. As previously
ACS Paragon Plus Environment

Page 7 of 35
Journal of Medicinal Chemistry
1
2
3
4
5
6
7
8
observed17 with 2 (TAFMT) and 3 (LY121260) (Figure 1), association transients were mainly
monophasic. The equilibrium KD was calculated from the onꢀ and offꢀrate constants. Data are
summarized in Table 1.
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
Residence time prediction. Initially, we focused on molecules for which crystal structures in
complex with GK1 were available, namely compounds 1 and 2 (Figure 1). The presence of a wellꢀ
defined crystal structure is crucial for determining the physicochemical properties of binding; hence,
it is also a fundamental prerequisite for applying our method. Moreover, we decided to include in this
first round of calculations ligand 3 (Figure 1) that shares a scaffold similar to 1 and size with the
aforementioned molecules. We therefore investigated these three systems first, to challenge the
method with similar ligands and to understand the level of scaling needed to investigate whole series
of ligands.
Ligands 1 and 2 were coꢀcrystallized18 with GK1 (1, resolution 1.7 Å, PDB ID 4NO7; 2, resolution
2.3 Å, PDB ID 3F9M) without strongly altering the protein structure, which maintains the same
folding features of previously resolved GK1 structures (Figure 2a). These ligands bind the same
region of the protein, namely the pocket comprised between the C terminus, the 206 – 218 helix, and
the 65 – 69 loop. Despite its length, the 65 – 69 loop is wellꢀresolved in the experimental Xꢀray
structure and its geometrical features are clearly visible (Figure 2b, 2c). In the apoꢀ form of GK1 (e.g.
3IDH, resolution equal to 2.14 Å), the heavy atoms of this loop systematically show a higher B factor
relative to the holoꢀ forms. This clearly indicates that the presence of the ligand reduces the
conformational flexibility of this region. We assumed that this feature is retained also for the
complex of GK1 with ligand 3; hence we placed it in the binding site by superimposing its structure
with the scaffold of ligands 1 and 2 (Figure 2d).
We simulated the three systems using two different scaling factors, λ = 0.4 and λ = 0.5. A more
pronounced potential scaling (λ = 0.4) pushed the system in a largely diffusive regime, whereas a
gentler potential scaling (λ = 0.5) required longer simulations to overcome barriers and sample
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 8 of 35
1
2
3
4
5
6
7
8
9
different minima. In this sense, the two different levels of scaling can highlight different details of the
ligandꢀprotein interaction during release. By varying the level of potential scaling, we observed the
passage from a regime where minima were still present and the energy landscape was slightly more
realistic (λ = 0.5) to a dominantly diffusive behavior (λ = 0.4) where the differences in minima are
flattened. As shown in Table 1 and in Figure 3, the residence times in both cases have the same trend,
i.e. tr (1) > tr (3) > tr (2), in very good agreement with the experimental data. A lower value of λ (i.e.
0.4) has the obvious advantage of reducing the computational cost, translating to a roughly 10x
speedup compared to a λ value of 0.5. However, when considering residence times in close proximity
(as here), a larger λ value increases separation and returns a higher confidence in the ranking.
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
We then sought to address whether new and structurally different ligands could be adapted starting
from the coordinates of a different holo structure, and how this strategy could lead to correct
predictions or inaccuracies. We selected two ligands, 4 (GKA50) and 5 (Sꢀ44520) (Figure 1), which
share a common structural core but differ in the size and nature of their substituents. We built the
complexes for compounds 4 and 5 by chemical replacement (see the Supporting Information S2a/b),
since experimental crystal structures for these molecules in complex with GK1 were not available.
Their starting conformations were obtained by adapting the structure of GK1 in complex with 1, i.e.
by matching the orientation of the new ligand with the experimentally available pose according to
their common scaffold. In principle, starting conformations could also be generated by molecular
docking. Here, we adopted the aforementioned procedure because a conserved orientation within the
binding site is consistently displayed by quite large series of GK1 inhibitors, for which the threeꢀ
dimensional structure is available. We were thus confident that our molecules should have been
oriented accordingly. We performed a series of unbinding simulations of 5 with two λ factors and
compared the computational predictions (tr ~ 14 ns for λ = 0.4, tr ~ 45 ns for λ = 0.5) with the
experimental results obtained via fluorescence stoppedꢀflow ligand displacement kinetics (see the
Materials and Methods section for details). The correlation was strictly dependent on the λ factor
ACS Paragon Plus Environment

Page 9 of 35
Journal of Medicinal Chemistry
1
2
3
4
5
6
7
8
(Figure 3). In detail, simulations with λ = 0.5 showed a residence time that, as confirmed
experimentally, was between those of 1 and 3, whereas λ = 0.4 led to a rather high value for ligand 5
with respect to the other compounds.
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
This case shows that the effect of a larger molecular size (the radius of gyration of ligand 5 is 5.3 Å,
whereas the one for ligands 1, 2 and 3 is ~ 4.3 Å) becomes dominant if the scaling factor is too low
(i.e., the potential is too heavily scaled) and the related energy landscape goes flat. This suggests that,
if the potential energy scaling is too pronounced, the proposed method could provide residence times
that are longer than those expected according to the experimental results. When treating molecules
that are partially diverse in terms of structure from those present in the crystal structure, an
explorative approach could rely on a double run with two scaling factors and check consensus and
possibly discrepancy between them. An even higher and unexpected estimate of residence time was
found for 4 with λ = 0.4, for which the residence time value was almost equal to that of 5. However,
whereas computational predictions for 5 were in line with the experimentally determined residence
time when λ = 0.5, 4 could not be accurately ranked with respect to the other members of the series
even when the higher scaling factor was applied. This behavior could partially be related to the size
of the molecule: since it is even larger than 5, 4 (its radius of gyration is 5.8 Å) displays a less
diffusive behavior regardless of the λ scaling and was hence assigned a longer residence time even
when the system potential was more gently scaled (λ = 0.5). As recently reported,8 a considerable
change in size within a series of congeneric ligands of the same protein can greatly alter the
solvation/desolvation enthalpy of the unbinding process, greatly affecting the ability of the method to
predict the relative residence time for a much larger molecule. Together with a notable deviation
from the linearity of distributed data in similar cases (see the recent example of Grp788), this results
in a lower predictability of the molecule’s residence time. Moreover, despite its similarity with ligand
5, ligand 4 displays more flexibility and degrees of freedom, suggesting that these also could be
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 10 of 35
1
2
3
4
5
important physiochemical features that are able to influence the kinetic properties of a ligand during
its release from the binding site.
6
7
8
9
Finally, we focused on two racemic mixtures of enantiomers 6 (Sꢀ49164) and 7 (Sꢀ49513) (6aꢀ6b and
7aꢀ7b in Figure 1; see also the Supporting Information S2c/d), which required simulations on each
enantiomer to computationally predict residence times. In particular, we carried out 40 runs for each
pair of enantiomers, estimating the residence time averaging all of them. These compounds displayed
an overall structure that is partially shared by 3 (Figure 2d). As in the previous case, we adapted the
new structures on an existing crystal structure (GK1 in complex with 3), superimposing the amide
moiety and the heterocyclic ring, and keeping their orientation in the binding pocket. The relative
prioritization of the ligands is retained in the computational model, averaging the residence times
obtained by simulating both the enantiomers (Figure 3, Table 1). Again, agreement with the
experimental data was fairly good also for these compounds.
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
We challenged our novel scaledꢀMDꢀbased protocol in two different scenarios: i) ranking unbinding
kinetics for a series of strictly congeneric compounds, and ii) properly ranking and predicting kinetics
for structurally unrelated molecules. We further showed that the method can reliably rank ligands of
the same series bearing minor chemical modifications, in agreement with our recent report.8 When
the chemical similarity is lower, because of differences in size and in the chemical groups contacting
the target, the method provides predictions that are less accurate but still reliable (see below). With
the sole exception of 4, we can observe rather linear correlations for both the λ values, with the noise
increasing, as expected, for λ = 0.4 relative to λ = 0.5 (correlation coefficients for residence times
normalized as recently described8 (details are reported in the section S1 of the Supporting
Information) are 0.72 and 0.91 for λ = 0.4 and λ = 0.5, respectively; see Figure 3b). At the lower λ
value, the more diffusive nature of simulations could be responsible for larger errors, pointing to the
need for more accurate and CPUꢀintensive simulations to achieve more reliable predictions of
unbinding kinetics. This is perfectly in line with expectations of the approach, which could be
ACS Paragon Plus Environment

Page 11 of 35
Journal of Medicinal Chemistry
1
2
3
4
5
6
7
8
adopted using lower λ values for returning fast results in screening campaigns (especially when the
ligands display a high degree of chemical similarity), whereas higher λ values could be adopted for
accurate but more computationallyꢀdemanding predictions.
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
Structure-kinetics relationships (SKRs). Despite the loss of detail that is associated with scaled
MD simulations, it is still possible to identify the chemical features responsible for longer residence
times, thus establishing structureꢀkinetics relationships (SKRs). This could be of valuable help in
designing novel GK1 inhibitors with an improved kinetics profile.
Ligands 1, 2, and 3 share some similarities in the chemical structure (i.e. an amide moiety connected
at the carbonyl side to a 5ꢀ or 6ꢀatom aromatic ring). In particular, the amide group provides stability
to the ligandꢀprotein complex donating a hydrogen bond to the backbone of R63 (Fig. 4a). However,
due to several differences in flexibility and steric hindrance, these compounds overall establish quite
heterogeneous patterns of interaction (Figure 1cꢀe). In this respect, both 1 and 3, together with an
increased number of rotatable bonds, display a common polar sulphonyl group, which is missing in 2.
This group is accommodated in a region of the protein where it can contact multiple polar groups
(Figure 4a), creating a network of interactions that stabilizes the protein–ligand complex and that is
eventually conducive to longer computed residence times of 1 and 3 with respect to 2. This is further
confirmed by compounds 4 and 5, which, bearing electron donor oxygen atoms in the same region of
the sulphonylic groups of ligands 1 and 3, were associated with longer computed residence times
across the series. The larger size of ligand 4 and 5 might also contribute to extending computed
residence time even further, increasing the enthalpic contribution to the overall unbinding activation
energy.
In addition, we noticed that the shape of the ligands plays a role in influencing the residence times.
Interestingly, the ligands that display a more pronounced linear shape (2, 6, and 7) exhibit a
significantly shorter residence time than those with a Tꢀshaped geometry (1, 3, 4, and 5), with an
average tr ~ 27 ns for the former molecules and ~ 90 ns for the latter ones. The GK activators binding
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 12 of 35
1
2
3
4
5
6
7
8
9
site has a rather linear shape, and hence it is easily accessible to ligands that possess a linear
geometry, whereas molecules bearing other shapes than linear require an induced fit binding
mechanism. In particular, Tyr125 side chain should move apart to accommodate Tꢀshaped ligands, a
geometrical rearrangement that can likely help accounting for the longer residence times experienced
by Tꢀshaped molecules relative to the linear ones. Therefore, the present computational approach can
help discern among rather different scaffolds, and hence prioritize, for subsequent chemical
synthesis, those most promising from a residence time and lead efficacy standpoints.
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
Finally, we turned our attention to aliphatic rings of 1 and 3 (cyclopentanone in 1, cyclohexane in 3),
which were docked between T65 and Y215, and that pointed toward the solvent. These rings are
located on the other side of the loop T64 – V69 (see Figure 1) with respect to the sulphonyl
functional group in ligands 1 and 3. For ligand 2, the same region of the protein is occupied by an Nꢀ
methylꢀimidazole ring (which has much greater affinity for water molecules than the aforementioned
cycles) and by moieties rich in methyl groups in ligands 4 and 5. In this respect, a more hydrophobic
element of the ligand in this region seems to provide more stabilization to the complex and to be
another key determinant for the unbinding kinetics.
To understand in more detail which structural elements of the protein are relevant for the unbinding
process, we analyzed all the trajectories, considering their statistical distribution in terms of contacts
between ligands and single residues. Considering both scaling potentials adopted in the present study
(i.e. λ = 0.4 or λ = 0.5), we identified a set of residues that were transiently involved in all the
unbinding processes observed (Figure 5 and Figure S3). We considered a residue relevant for
unbinding when a contact with the ligand occurred in more than 66% of the replicas, for both sets of
simulations. We could thus monitor at a glance the unbinding events providing a residue level
description of the dissociation paths. As expected, simulations with λ = 0.4 had a more pronounced
diffusive character and thus provided more transient and aspecific proteinꢀligand contacts.
Conversely, simulations with λ = 0.5 contributed more significantly to the identification of specific
ACS Paragon Plus Environment

Page 13 of 35
Journal of Medicinal Chemistry
1
2
3
4
5
6
7
8
9
contacts along the unbinding paths. The contact maps reported in Figure 6 clearly exemplify this dual
behavior (all the other maps are available in Figures S4 and S5 of the Supporting Information). Our
analysis allowed the identification of two groups of amino acids: one encompassing residues that
were trivially clustered around the binding site (residues 56ꢀ66, 89ꢀ95, 154ꢀ156, 196ꢀ199, 203, 206ꢀ
218, 231ꢀ234, 444ꢀ454), and the other formed by a group of amino acids (residues 241ꢀ248, 393ꢀ394)
located roughly 20 to 25 Å away from the binding site.
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
Diabetes is heavily characterized by genetic mutations that correspond to a huge variety of
phenotypes. Interestingly, the most relevant pathological mutations identified in diabetes and glucose
homeostasis19 (Figure 5b) correspond to residues that, according to our analysis, appear to be relevant
in controlling the unbinding process. GKAs can only be accommodated in their binding site once the
glucose is bound to the protein via an allosteric modulation. At the same time, the protein’s
enzymatic efficiency is significantly reduced for many GK1 mutants19 that are related to diabetes and
similar metabolic diseases. Hence, good GKAs could be designed by accounting for the extreme
variability of the binding site and the portion of the surrounding region that is relevant for unbinding.
For example, mutants V91L, W99R, and Y214C (Figure 5b) severely unbalance the biochemistry of
glucose in patients affected by diabetes. Hence, their effect must be related to the change in local
dynamics and to steric hindrance in regulating the action of any potential GKA as well as their
dissociation.
Conclusions
Residence time is emerging as a key parameter for analyzing and optimizing lead candidates with
improved in vivo efficacy, within modern drug discovery programs that take it into account during
the hitꢀtoꢀlead and lead optimization phases. Here, we have shown that residence time predictions are
possible using scaledꢀMDꢀbased approaches, and that the computational protocol can be fine tuned
by comparing computational data and experimental kinetics measurements. Different scaling factors
can be used and, depending on the ligand set, the best tradeoff between speed and accuracy can be
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 14 of 35
1
2
3
4
5
6
7
8
9
sought. To properly identify scaling factors, we compared computational predictions with
experimental data, with λ = 0.5 emerging as more suitable for good correlations (relative to λ = 0.4).
This was mainly because the series of ligands included compounds that were not strictly analogous.
However, when investigating series of structurally analogous molecules, a lower λ value can be used,
allowing faster simulations.8 Building and expanding on the previously reported SMDꢀbased
method,8 the present work suggests, a practical approach to set the ideal scaling factor to adopt on the
basis of the needs of the users and the chemical nature of the ligands. We subsequently focused on
SKRs, an accurate description of the structural determinants responsible for a good residence time
can be of paramount importance when designing new compounds with improved in vivo profiles.
SKRs can be depicted from a ligandꢀbased standpoint, from a structureꢀbased analysis, or from both,
provided that structural information for the ligand and target is available. In this context, by means of
a statistical analysis of proteinꢀligand contacts during the unbinding events, we have been able to
identify those residues that affect the residence time and overall kinetic and release profile of a
compound. In addition, we observed that the ligand shape could be responsible for an induced fit
effect, which could explain the different residence times of differently shaped compounds. In
conclusion, the present computational methodology is very promising and practical. However, further
case studies are needed to highlight the potential and limitations of scaled MD simulations for drug
binding kinetics studies, with the goal of making this approach the method of choice for estimating
and predicting residence times.
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
Experimental section
Chemical synthesis
The synthesis of ligands 1ꢀ4 was already reported.17ꢀ18 The ones of ligands 6, 7 and 5 is reported
below, with reference (numbers in italics) to the schemes of the synthetic routes in the Supporting
Information section S6. All the ligands have been purified using HPLC and characterized using NMR
spectroscopy and MS, assessing for all the molecules purity more than 95%.
ACS Paragon Plus Environment

Page 15 of 35
Journal of Medicinal Chemistry
1
2
3
4
5
6
7
8
6-[[(1R,1aR,6aR)-1a-methyl-6,6a-dihydro-1H-cyclopropa[a]indene-1-carbonyl]amino]pyridine-
3-carboxylic acid (6). Cyclopropanation reaction of 3ꢀmethylindene 8 with ethyl diazoacetate was
performed in the presence of copper sulfate to afford a mixture of diastereoisomers 9 and 10 (70/30).
After separation by chromatography on silica gel, the ester function of diastereoisomer 9 was
hydrolyzed to the corresponding acid 11. Compound 11 was then transformed in acid chloride with
oxalyl chloride, and condensed with 2ꢀpyridin acid to afford the desired amide 12 as a racemic
mixture (compound 6 in the main text), with 98.7% HPLC purity.
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
1H NMR (400 MHz, 300K, DMSOꢀd6): δ 13.10 ( m, OH), 10.86 (s, NH), 8.78 (s, 1H), 8.22 (m, 2H),
7.37 (d, 1H), 7.19 (m, 3H), 3.24 (dd, 1H), 2.98 (d, 1H), 2.38 (dd, 1H), 1.66 (d, 1H), 1.62 (s, 3H);
HRMS (ESI/FIA) for C18H16N2O3: [M+H]+(calculated) = 309.123918, [M+H]+(found) = 309.1231.
6-[[(1R,1aR,6aR)-6a-methyl-1a,6-dihydro-1H-cyclopropa[a]indene-1-carbonyl]amino]pyridine-
3-carboxylic acid (7). Cyclopropanation reaction of 2ꢀmethylindene 13 with ethyl diazoacetate was
performed in the presence of rhodium acetate to afford a mixture of diastereoisomers 14 and 15
(63/37). After separation by chromatography on silica gel, the ester function of diastereoisomer 7 was
hydrolyzed to obtain acid 16. The acid function was then transformed in acid chloride with oxalyl
chloride, and condensed with 2ꢀaminopyridine methyl ester to afford, after saponification, the desired
amide 17 (compound 7a/b in the main text), as a racemic mixture, in 99% purity (HPLC).
1H NMR (400 MHz, 300K, DMSOꢀd6): δ 13.11 (m, OH), 10.93 (s, NH), 8.79 (s, 1H), 8.23 (m, 2H),
7.33 (m, 1H), 7.20 (m, 1H), 7.11 (m, 2H), 3.13 (dd, 2H), 2.84 (d, 1H), 1.70 (d, 1H), 1.48 (s, 3H);
HRMS (ESI/FIA) for C18H16N2O3: [M+H]+(calculated) = 309.123918, [M+H]+(found) = 309.1255.
(E)-3-[6-[(3,5-diisopropoxybenzoyl)amino]-3-pyridyl]prop-2-enoic acid (5). 3,5ꢀdihydroxy
methyl benzoate 18 is dialkylated by isopropyl iodide to obtain the ester 19, hydrolyzed to acid 20.
The acid function is transformed in acid chloride with oxalyl chloride, then condensed with 2ꢀaminoꢀ
pyridineꢀ5ꢀacrylate 22, obtained from 2ꢀaminoꢀ5ꢀiodoꢀpyridine 21 and methyl acrylate under Heck
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 16 of 35
1
2
3
4
5
conditions. Resulting condensation ester 23 is finally hydrolyzed to afford the desired acid 24
(compound 5 in the main text), in 99% purity (HPLC).
6
7
8
9
1H NMR (400 MHz, 300K, DMSOꢀd6): δ 12.5 (s, OH), 8.7 (s, 1H), 7.6 (d, J=15Hz, 1H), 6.6 (d, 1H),
4.7 (m, 2H); HRMS (ESI/FIA) for C21H24N2O5: [M+H]+(calculated) = 385.176348, [M+H]+(found) =
385.1771.
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
Protein expression and purification
Human glucokinase isoform 1 (pancreatic) was expressed and purified as previously described.17
Steady-state kinetics
GK activity was measured by monitoring the rate of G6P formation using the G6PDH/NADP
(glucoseꢀ6ꢀphosphate dehydrogenase/nicotinamide adenine dinucleotide phosphate)ꢀcoupled enzyme
assay. Initial rate measurements were carried out at 25 °C in buffer A (HEPES/Na 50 mM, NaCl 100
mM, MgCl2 5 mM, TCEP 2 mM, pH 7.1) on a PHERAstar microplate reader (BMG Labtech) by
following the increase of absorbance at 340 nm due to the reduction of NADP. Assays were
performed in a 384ꢀwell microplate with a final volume of 50 µl per well. To determine glucose or
MgꢀATPꢀrelated kinetic parameters, glucose concentration was varied (0ꢀ50 mM) while maintaining
a saturating (4 mM) MgꢀATP concentration. A lag phase observed in the early time course of the
reaction was omitted for the determination of the initial rate. This lag phase is likely due to the
coupledꢀenzyme assay, as already described.17 Initial rate data obtained with varying glucose
concentrations were fitted to the Hill equation. Nonlinear leastꢀsquares regression analysis was
performed with SigmaPlot 9.01.
Dissociation rate measurements by competitive displacement kinetics
ACS Paragon Plus Environment

Page 17 of 35
Journal of Medicinal Chemistry
1
2
3
4
5
6
7
8
9
Because of the relatively fast dissociation kinetics of the GKA studied, the stoppedꢀflow rapid mixing
technique was used, with the same setup described for the association kinetics measurements. For all
GKA except ligand 2, displacement experiments were performed by mixing an equal volume (75 µL)
of the preꢀequilibrated GK/GKA complex with an excess of competing ligand 2 in buffer A,
complemented with 5% DMSO and 100 mM Glc. GK, GKA, and ligand 2 concentrations after
mixing were 1 µM, 5 µM, and 50 µM respectively. Fluorescence traces were analyzed as previously
described.17
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
Protein crystallization and X-ray diffraction data collection
The crystallization samples were solutions of GK (at 6 or 10 mg ml/1) in 20mM HEPES, 50 mM KCl,
2 mM tris(2ꢀcarboxyethyl)phosphine (TCEP), 100 mM glucose, 5% (w/v) glycerol pH 7.5. When
present, activators (1 and 2) were diluted to a concentration of 80 mM and added to the protein
solution prior to concentration and AMPꢀPNP was added to a final concentration of 10 mM with the
addition of 1 mM MgCl2.
Initial crystals were obtained under conditions similar to those already reported14 and showed the
same unit cell parameters. However, in our hands, these crystals did not diffract well and grew very
slowly. After extensive screening in a matrix of various PEGs buffered using HEPES pH 7.5, a
different crystal form was obtained with space group P212121. This crystal was further used for
preliminary seeding experiments.
All crystals used in this work were grown under the same conditions, at room temperature using
hangingꢀdrop vapor diffusion in Linbro plates. 2 ml protein solution was mixed with 1 ml reservoir
solution (3232ro plates. 2 p 100 mM HEPES pH 7.5). After 3 d of equilibration, the drops were
seeded with crushed crystals diluted in reservoir solution.
Crystals appeared within 12 h in the presence of TAFMT and notably more slowly in its absence.
They grew to their final dimensions within 2 d (with TAFMT) or longer (without TAFMT). Crystals
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 18 of 35
1
2
3
4
5
were obtained over a range (32–38%) of PEG concentrations and cooled under liquid nitrogen
without the need for further cryoprotection.
6
7
8
9
Data sets were collected on the SLS PX6d beamline (Villigen, Switzerland) and were reduced using
XDS.20 To avoid any bias on the densities of the ligands, the following procedure was used for
molecular replacement and structure refinement. Firstly, molecular replacement was performed using
the program MOLREP.21 As a model, we used the structure 1V4S,14 from which we removed all
ligands or bound waters. This was followed by refinement omitting all ligands using REFMAC.22
After several refinement cycles, water molecules were built iteratively by ARP/wARP.23 At this stage,
the densities of the ligands were clearly visible and were then refined. Further refinement using TLS
was also performed.24 The structures have been deposited with the PDB ID 4NO718 (1) and 3F9M18
(2), the latter being a refinement of the structure 3D18.
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
MD and scaled MD simulations
Each compound was geometrically optimized via a quantum mechanical approach: electron density
calculations were performed in NWChem25 using the basis set 6ꢀ31G* at the HartreeꢀFock (HF) level
of theory. Partial charges were derived using the RESP methodology26 as implemented in
Antechamber, leading via a GAFF parameterization to a complete topological description of each
ligand to be used for classical simulations. The coordinates of the HF/6ꢀ31G* optimized ligands are
reported in the Supporting Information, section S7.
Proteinꢀligand experimental complexes were used as a starting point of the MD simulations. When
the experimental proteinꢀligand complex structure was not available the ligand was placed in the
binding site according to its best superimposition with the most corresponding moiety of the
available experimental structures (PDB ID: 4NO7 and 3F9M).
ACS Paragon Plus Environment

Page 19 of 35
Journal of Medicinal Chemistry
1
2
3
4
5
6
7
8
The proteinꢀligand complexes were then used as a starting point for molecular dynamics simulations
performed in a GROMACS 4.6.127 version customized to perform scaled molecular dynamics (scaled
MD)28 implemented by BiKi Life Sciences (BiKi Technologies, Genoa, Italy) as recently described.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
All the complexes were placed in the geometrical center of parallelepipedꢀshaped boxes of volume
equal to ~ 650 nm3. The simulation boxes were then solvated using tLeap, with ~70000 TIP3P water
molecules.29 Some water molecules were replaced with sodium ions in order to preserve the electroꢀ
neutrality of the system according to need, i.e. the charge of the protein plus the charge of the ligand,
which varied according to the case considered. The system was minimized with the steepest descent
method, followed by equilibration of the restrained protein (isotropic 1000ꢁkJ mol−1ꢁnm−1 force
applied to each heavy atom of the protein backbone) in NPT (up to 400ꢁps, pressureꢁ=ꢁ1ꢁatm) and
NVT (up to 400ꢁps) ensembles at 300ꢁK via a standard MD procedure. Electrostatics were treated
with the cutoff method for shortꢀrange interactions and with the Particle Mesh Ewald method for the
longꢀrange ones (rlistꢁ=ꢁ1ꢁnm, cutoff distanceꢁ=ꢁ0.9ꢁnm, vdW distanceꢁ=ꢁ0.9ꢁnm, PME orderꢁ=ꢁ4).30 The
constant temperature conditions were provided by using the Vꢀrescale thermostat,31 a modification
from Berendsen’s coupling algorithm.
A series of partially unrestrained (see below) scaled MD production runs were performed for each
complex until the occurrence of the unbinding event, defined as the situation where interactions
between the ligand and the binding site are no longer present, corresponding to a distance between
the ligand and the binding site centers of mass of 30ꢁÅ. This threshold has been chosen in order to
achieve complete solvation of the ligand (no remaining hydrogen bonds nor contacts between the
small molecule and the protein) mimicking the experimental conditions adopted for the experimental
measurement of the koff values. Twenty simulations were performed for each proteinꢀligand system
and for each scaling factor λ. We used two different scaling factors, λ = 0.4 and λ = 0.5, thus having
40 simulations for each complex and 280 simulations overall.
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 20 of 35
1
2
3
4
5
6
7
8
9
The restraints were adopted as previously described,8 i.e. weakly (50ꢁkJ mol−1 nm−1) restraining all
the backbone heavy atoms with the exception of those in the binding site, defined as those belonging
to the residues within 6ꢁÅ of the surface of the ligand computed on the starting crystal structure
(without hydrogen atoms). The residues left unrestrained were: 57 58 59 60 61 62 63 64 65 92 93 94
155 206 208 210 211 214 216 217 218 231 248 447 448 449 451 452 453 454. The unbound state
was determined using the same criterion of Mollica et al.8
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
All the simulations were set up using the BiKi software package and performed on a set of inꢀhouse
machines, equipped with two esacore Intel Xeon processors and 2 NVIDIA GTX 780ꢁGPUs, for a
total of 1300ꢁCPU days.
Abbreviations used
ADP, adenosine diphosphate; AMBER, assisted model building with energy refinement; ATP,
adenosine triphosphate; CD, circular dichroism; CPU, central processing unit; GPU, graphical
processing unit; DMSO, dimethyl sulphoxide; GAFF, generalized AMBER force field; GK,
glucokinase; GKA, glucokinase activator; GROMACS, Groningen machine for chemical
simulations; HEPES 4ꢀ(2ꢀhydroxyethyl)ꢀ1ꢀpiperazine ethanesulfonic acid; HF, HartreeꢀFock; HK,
hexokinase; HPLC, high performance liquid chromatography; MD, molecular dynamics; MS, mass
spectrometry; NADP, nicotinamide adenine dinucleotide phosphate; NMR, nuclear magnetic
resonance; NVT, canonical ensemble (amount of substance (N), volume (V), and temperature (T) are
conserved); NPT, isothermalꢀisobaric ensemble (amount of substance (N), temperature (T), and
pressure (P) are conserved); PES, potential energy surface; PME, particle mesh Ewald; RESP,
restrained electrostatic potential; SKR, structure kinetics relationship; SMD, scaled molecular
dynamics; SPR, surface plasmon resonance; T2DM, type 2 diabetes mellitus; TCEP, tris(2ꢀ
carboxyethyl)phosphine; vdW, van der Waals.
ACS Paragon Plus Environment

Page 21 of 35
Journal of Medicinal Chemistry
1
2
3
Supporting Information
4
5
6
7
8
Mathematical derivation of the koff and residence time scaling, starting poses for ligands 4, 5, 6a, and
7a, surface representation of the unbinding paths in the GK1 structure, ligand protein contact maps
during the unbinding process for the single ligands, ligand protein contact maps during the unbinding
process for the enantiomeric mixtures, scheme of reactions for the synthesis of compounds 7, 6 and 5,
ligand’s optimized structure coordinates with RESP charges in .mol2 format.
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
Corresponding Authors Information
Phone: +39 010 71781530
Fax: +39 010 7170187
Eꢀmail: andrea.cavalli@iit.it
Phone: +33 (0)1 55 72 28 15
Fax: +33 (0)1 55 72 52 15
Eꢀmail: pierre.ducrot@servier.com
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 22 of 35
1
2
3
4
5
6
7
8
9
References
1. Copeland, R. A.; Pompliano, D. L.; Meek, T. D. Drugꢀtarget residence time and its
implications for lead optimization. Nat. Rev. Drug Discov. 2006, 5, 730ꢀ739.
2. Swinney, D. C. Biochemical mechanisms of drug action: what does it take for success? Nat.
Rev. Drug Discov. 2004, 3, 801ꢀ808.
3. De Vivo, M.; Masetti, M.; Bottegoni, G.; Cavalli, A. Role of molecular dynamics and related
methods in drug discovery. J. Med. Chem. 2016, 59, 4035ꢀ4061.
4. Decherchi, S.; Berteotti, A.; Bottegoni, G.; Rocchia, W.; Cavalli, A. The ligand binding
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
mechanism to purine nucleoside phosphorylase elucidated via molecular dynamics and machine
learning. Nat. Commun. 2015, 6, 6155.
5.
Buch, I.; Giorgino, T.; De Fabritiis, G. Complete reconstruction of an enzymeꢀinhibitor
binding process by molecular dynamics simulations. Proc. Natl. Acad. Sci. U. S. A. 2011, 108,
10184ꢀ10189.
6.
drug molecule find its target binding site? J. Am. Chem. Soc. 2011, 133, 9181ꢀ9183.
7. Tiwary, P.; Limongelli, V.; Salvalaglio, M.; Parrinello, M. Kinetics of proteinꢀligand
Shan, Y.; Kim, E. T.; Eastwood, M. P.; Dror, R. O.; Seeliger, M. A.; Shaw, D. E. How does a
unbinding: Predicting pathways, rates, and rateꢀlimiting steps. Proc. Natl. Acad. Sci. U. S. A. 2015,
112, E386ꢀ391.
8.
Mollica, L.; Decherchi, S.; Zia, S. R.; Gaspari, R.; Cavalli, A.; Rocchia, W. Kinetics of
proteinꢀligand unbinding via smoothed potential molecular dynamics simulations. Sci. Rep. 2015, 5,
11539.
9.
pancreatic betaꢀcells and hepatocytes. Diabetes 1990, 39, 647ꢀ652.
10. Wild, S.; Roglic, G.; Green, A.; Sicree, R.; King, H. Global prevalence of diabetes: estimates
Matschinsky, F. M. Glucokinase as glucose sensor and metabolic signal generator in
for the year 2000 and projections for 2030. Diabetes Care 2004, 27, 1047ꢀ1053.
ACS Paragon Plus Environment

Page 23 of 35
Journal of Medicinal Chemistry
1
2
3
11.
Zimmet, P.; Alberti, K. G.; Shaw, J. Global and societal implications of the diabetes
4
5
6
7
8
epidemic. Nature 2001, 414, 782ꢀ787.
12. Agius, L.; Peak, M.; Newgard, C. B.; GomezꢀFoix, A. M.; Guinovart, J. J. Evidence for a role
9
of glucoseꢀinduced translocation of glucokinase in the control of hepatic glycogen synthesis. J. Biol.
Chem. 1996, 271, 30479ꢀ30486.
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
13.
Heredia, V. V.; Thomson, J.; Nettleton, D.; Sun, S. Glucoseꢀinduced conformational changes
in glucokinase mediate allosteric regulation: transient kinetic analysis. Biochemistry 2006, 45, 7553ꢀ
7562.
14.
regulation of the monomeric allosteric enzyme human glucokinase. Structure 2004, 12, 429ꢀ438.
15. Grimsby, J.; Berthel, S. J.; Sarabu, R. Glucokinase activators for the potential treatment of
type 2 diabetes. Curr. Top. Med. Chem. 2008, 8, 1524ꢀ1532.
16. Sarabu, R.; Berthel, S. J.; Kester, R. F.; Tilley, J. W. Novel glucokinase activators: a patent
review (2008 ꢀ 2010). Expert Opin. Ther. Pat. 2011, 21, 13ꢀ33.
17. Antoine, M.; Boutin, J. A.; Ferry, G. Binding kinetics of glucose and allosteric activators to
human glucokinase reveal multiple conformational states. Biochemistry 2009, 48, 5466ꢀ5482.
18. Petit, P.; Antoine, M.; Ferry, G.; Boutin, J. A.; Lagarde, A.; Gluais, L.; Vincentelli, R.;
Kamata, K.; Mitsuya, M.; Nishimura, T.; Eiki, J.; Nagata, Y. Structural basis for allosteric
Vuillard, L. The active conformation of human glucokinase is not altered by allosteric activators.
Acta Crystallogr. D Biol. Crystallogr. 2011, 67, 929ꢀ935.
19.
Rev. Drug Discov. 2009, 8, 399ꢀ416.
20. Kabsch, W. Integration, scaling, spaceꢀgroup assignment and postꢀrefinement. Acta
Crystallogr. D Biol. Crystallogr. 2010, 66, 133ꢀ144.
21. Vagin, A.; Teplyakov, A. Molecular replacement with MOLREP. Acta Crystallogr. D Biol.
Crystallogr. 2010, 66, 22ꢀ25.
Matschinsky, F. M. Assessing the potential of glucokinase activators in diabetes therapy. Nat.
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 24 of 35
1
2
3
22.
Murshudov, G. N.; Skubak, P.; Lebedev, A. A.; Pannu, N. S.; Steiner, R. A.; Nicholls, R. A.;
4
5
6
Winn, M. D.; Long, F.; Vagin, A. A. REFMAC5 for the refinement of macromolecular crystal
7
8
structures. Acta Crystallogr. D Biol. Crystallogr. 2011, 67, 355ꢀ367.
9
23.
Cohen, S. X.; Ben Jelloul, M.; Long, F.; Vagin, A.; Knipscheer, P.; Lebbink, J.; Sixma, T. K.;
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
Lamzin, V. S.; Murshudov, G. N.; Perrakis, A. ARP/wARP and molecular replacement: the next
generation. Acta Crystallogr. D Biol. Crystallogr. 2008, 64, 49ꢀ60.
24.
groups undergoing TLS motion. Acta Crystallogr. D Biol. Crystallogr. 2006, 62, 439ꢀ450.
25. Valiev M., B. E. J., Govind, N., Kowalski, K., Straatsma, T.P., van Dam, H.J.J., Wang, D.,
Painter, J.; Merritt, E. A. Optimal description of a protein structure in terms of multiple
Nieplocha, J., Apra, E., Windus, T.L., de Jong W.A. . NWChem: a comprehensive and scalable openꢀ
source solution for large scale molecular simulations. Comput. Phys. Commun. 2010, 181, 1477.
26.
Bailey, C. I., Cieplak, P., Cornell, W. D. & Kollman, P. A. A wellꢀbehaved electrostatic
potential based method using charge restraints for deriving atomic charges: the RESP model. J. Phys.
Chem. 1993, 97, 10269–10280
27.
Pronk, S.; Pall, S.; Schulz, R.; Larsson, P.; Bjelkmar, P.; Apostolov, R.; Shirts, M. R.; Smith,
J. C.; Kasson, P. M.; van der Spoel, D.; Hess, B.; Lindahl, E. GROMACS 4.5: a highꢀthroughput and
highly parallel open source molecular simulation toolkit. Bioinformatics 2013, 29, 845ꢀ854.
28.
scaled molecular dynamics. J. Phys. Chem. B 2013, 117, 12759ꢀ12768.
29. Jorgensen, W. L., Chandrasekhar, J., Madura, J. D., Impey, R. W. &Klein, M. L. .
Sinko, W.; Miao, Y.; de Oliveira, C. A.; McCammon, J. A. Population based reweighting of
Comparison of simple potential functions for simulating liquid water. J. Chem. Phys. 1983, 79, 926–
935
30.
Darden, T.; Perera, L.; Li, L.; Pedersen, L. New tricks for modelers from the crystallography
toolkit: the particle mesh Ewald algorithm and its use in nucleic acid simulations. Structure 1999, 7,
R55ꢀ60.
ACS Paragon Plus Environment

Page 25 of 35
Journal of Medicinal Chemistry
1
2
3
31.
Bussi, G.; Donadio, D.; Parrinello, M. Canonical sampling through velocity rescaling. J.
4
5
6
Chem. Phys. 2007, 126, 014101.
7
8
32.
Efron, B. Bootstrap methods: Another look at the jackknife. Ann. Statist. 1979, 7, 1ꢀ26.
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
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 26 of 35
1
2
3
Table 1
4
5
6
7
8
Experimental Kd, kon, koff, and residence time (tr,exp), and in silico residence time (tr,comp) for each
compound examined. Kds are expressed in ꢀM, kon in sꢀ1Mꢀ1, koff in sꢀ1, tr,exp = 1 / koff in s, and tr,comp
in ns. The activity of the enzyme is reported (in percentage) with reference to the yield of glucose
conversion in the presence of the ligand (details are provided in Materials and Methods). The tr,comp
values are reported for both λ values used: an estimation of the error is reported as standard error (se)
and via a bootstrapping procedure (bs).32 The experimental measurements for ligands from 1 to 5
have been performed on the enantiomers indicated in Fig.1, whereas the data for ligands 6 and 7 have
been collected for the racemic mixtures.
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
Act.
(%)
λ = 0.4
λ = 0.5
Ligand
Kd
kon
koff
tr,exp
tr,comp,se
tr,comp,bs
5.5 ± 0.9
3.7 ± 0.3
5.0 ± 0.4
tr,comp,se
tr,comp,bs
105.1 ± 10.1
29.3 ± 5.3
38.9 ± 7.1
92.9 ± 7.3
99.7 ± 6.7
25.9 ± 3.9
24.7 ± 3.0
0.83 0.1x106 0.12
0.36 1.1x106 0.43
0.40 0.9x106 0.37
8.3
2.3
2.7
1.6
6.3
0.7
0.2
90
170
130
20
5.5 ± 0.8
3.8 ± 0.2
5.0 ± 0.5
105.4 ± 9.2
29.4 ± 5.1
38.7 ± 7.2
92.5 ± 7.5
99.4 ± 6.6
25.7 ± 3.4
24.6 ± 2.6
1
2
3
4
5
6
7
3.3
0.2x106 0.61
12.9 ± 2.9 12.8 ± 2.8
14.1 ± 0.7 14.3 ± 0.6
0.72 0.2x106 0.16
60
2.0
4.4
0.7x106
1.3x106
1.4
5.6
n.a.
n.a.
6.5 ± 1.1
6.2 ± 1.2
6.5 ± 1.0
6.3 ± 1.1
ACS Paragon Plus Environment

Page 27 of 35
Journal of Medicinal Chemistry
1
2
3
Figure 1
4
5
6
7
8
Chemical structures of the ligands used in complex with GK1. The protonated structures are reported
for the acidic ligands, whereas the ionic form was used in the simulation.
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 2
Xꢀray structure of the three complexes of GK1 with 1, 2 , and reconstructed pose of 3 . a) Overall
structure of GK1 in complex with 1 (red), 2 (blue), and 3 (yellow) superimposed in the binding site.
bꢀd) Local structure of the binding pocket for 1 (b), 2 (c), and 3 (d).
Figure 3
Normalized experimental vs. computational residence times for λ = 0.4 (a) and λ = 0.5 (b) (errors for
the computational residence times have been normalized as well, according to error propagation).
The ratios between experimental or computational data and a reference value have been plotted: the
residence time values for ligand 3 have been used as a reference due to their average placement
within the overall data distribution.8 The experimental residence times (on the abscissa) were scaled,
after the aforementioned referencing scaling, using the same factor λ used for scaling the potentials
during the SMD simulations. Analytical details of the whole procedure are reported in the section S1
of the Supporting Information.
A linear fit was applied to the whole datasets (regression lines are not shown) as well with the
exclusion of the ligand 4, leading to a correlation factor R equal to 0.41 (0.72 excluding ligand 4) for
λ = 0.4 and to 0.81 (0.92 excluding ligand 4) for λ = 0.5.
Figure 4
Binding site of GK1 and its interactions with ligand 1. Hydrogen bonding between the amide proton
of the ligand and the R63 backbone carbonyl is reported as a dotted line.
ACS Paragon Plus Environment

Journal of Medicinal Chemistry
Page 28 of 35
1
2
3
4
5
6
Figure 5
7
8
9
Unbinding processes for key GK1 residues (in red), identified as explained in the text, are annotated
on the protein 3D structure. a) An overall representation. b) A detailed representation displaying all
the contacts in proximity of the binding site. (b) also reports a set of residues (in green using a stick
representation; the spheres correspond to the backbone alpha carbon atoms), These residues
correspond to pathological mutants in T2DM, which severely affect the enzymatic properties of GK1.
Compound 1 is reported in yellow sticks. The residues highlighted in red are 56ꢀ66, 89, 92ꢀ95, 154ꢀ
156, 196ꢀ199, 203, 206ꢀ218, 231ꢀ234, 241ꢀ248, 393ꢀ394, and 444ꢀ454.
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
The monoꢀdimensional contact map between ligand 5 ) and GK1 is reported for two different levels
of scaling adopted during scaled MD simulations, λ = 0.4 (black) and λ = 0.5 (red). The numbering
of residues is reported on the abscissa and the percentage of contacts computed for all simulation
replicas is reported on the ordinate. Details of the method adopted for generating these maps are
described within the text.
ACS Paragon Plus Environment

Page 29 of 35
Journal of Medicinal Chemistry
1
2
Figure 1
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
ACS Paragon Plus Environment